定步长四阶经典公式 解决数值积分

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
* ]4 h7 s& [+ n( B7 G, E! r定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
% E+ Z) [7 _, ?! z4 i* K[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
) s5 {" e( Z& R" q1 {2 ~[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
9 ~( V5 o: Y6 @* p# ?5 F$ @[k4 = h \cdot f(tn + h, yn + k_3)]
' o( I- d" J/ `[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
2 V8 }' A1 ^# d# |4 N其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
& h$ z8 Z* I/ ]" ]这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m8 L% D* a9 d# `8 B* `
   
2. 
   4 w( O6 q1 f3 E+ ^: g' M
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)0 m\" R3 g& T8 T: G' G
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   : d/ X$ g# R, z) U
5.    disp('function z=f(x,y)');) d- j( m$ L, e0 \  A4 c1 ]. D
   
6.    disp('z=y-2*x/y;');1 O, S' t+ X5 i
   
7.    disp('并将该文件保存在work文件夹下');
   - [: ~2 y* ]3 i4 F4 w\" c
8. end
   3 l% N# _  F/ ?5 [! z/ k) v, M* S
9.   G. \$ g! Z9 m# V
   
10. X1=input('请输入求解区间的左端点X1=');: G. A& F2 }- i5 i3 C4 F% o6 W
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');; c8 e5 K# }5 \6 T
    
12. Xn=input('请输入求解区间的右端点Xn=');9 c0 q& t$ Y+ o. a
    
13. h=input('请输入求解步长h=');; G- N. G7 h, b0 V3 ^( R3 l% @
    
14. % o3 l9 A( G5 _. I1 b+ |: l4 S
    
15. X=X1;
    ' j# z* d- X, ?! f& A8 J. [
16. Y=Y1;                                                        %运算初始点
    ) b1 ?4 {2 E4 X# ]
17. n=0;                                                         %节点序号变量置零7 |% S8 Z( V) |& |7 V% G( u
    
18. 
    : J  a/ K* r' R. F
19. while X<=Xn-h
    8 k; x- z! x1 M0 i  ?1 n
20.     K1=f(X,Y);
    ! M9 [9 D\" I( R, h% w: a- C
21.     K2=f(X+h/2,Y+K1*h/2);! c5 P8 F' W9 h. O2 \5 T8 I9 H
    
22.     K3=f(X+h/2,Y+K2*h/2);
    1 J9 h; X7 ^) N# p- E
23.     K4=f(X+h,Y+K3*h);
    : l9 J$ B$ Y  X$ T, }
24.     X=X+h;- Q/ a0 X- ^9 a# J6 P
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式% \1 L\" u; Y1 e+ r
    
26.     n=n+1;                                                   %节点序号加1) c* G6 H3 U3 V
    
27. 
    + R( \9 O5 T, L, P& J
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    8 e; o: D6 _3 y+ }: g- G, g4 N
29.     plot(X,Y,'o'): c; Y* P$ T# @! V% D
    
30.     hold on
    + ?) ^% h5 M! j4 l( h
31. end

复制代码

1. function z=f(x,y)
   ( g. S- P* l2 C
2. z=y-2*x/y;

复制代码

" }: V! A- d! t$ J4 F8 {; X4 w) [