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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
: S* N# W" i7 e5 s) @" M& p+ [定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
( \" N2 {& t6 k6 \2 ~+ Q5 N& U这个方法的迭代公式如下：
. t% z3 N5 X' r[k1 = h \cdot f(tn, yn)]
7 U) e7 ~) G/ y) I2 b* I[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
" x; h- R/ {" l[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
* U9 X" S& p9 F+ `" ?4 C' n) D[k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
8 [4 S0 }" H- U4 `+ _6 I3 C, D3 d其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
4 N/ ^, p, J. D+ f+ k1 k1 y这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   / W8 J8 s; _0 G
2. : q6 B! ^9 b% e9 V! X  D
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   8 E6 e3 Z2 D3 W1 _- t( a, ^
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   ( w8 v' W# e9 H6 e1 N  A
5.    disp('function z=f(x,y)');- R$ v3 D1 x2 g4 G; r
   
6.    disp('z=y-2*x/y;');) V: R8 }) r% E, w
   
7.    disp('并将该文件保存在work文件夹下');3 T) d9 x6 t* k: J; k6 [\" ]
   
8. end $ {* z/ Q; l2 V
   
9. 
   ( i: X  ]+ r+ `\" y
10. X1=input('请输入求解区间的左端点X1=');\" P9 y) R  D) N( X\" N( q6 I) D
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    . r/ O/ V, ]- m2 l
12. Xn=input('请输入求解区间的右端点Xn=');$ G* G9 D, B; c$ s* l5 a, r
    
13. h=input('请输入求解步长h=');
    $ q! c5 ^' q; f
14. ' a! e! t: R% o7 w8 C# K
    
15. X=X1;- L1 o/ K. ]' y; r1 L
    
16. Y=Y1;                                                        %运算初始点
    9 M0 A  d9 x# m4 H0 i) g& k
17. n=0;                                                         %节点序号变量置零' W+ I& J+ C  T/ I+ t% T2 ^) \8 b
    
18. # T, k; r4 v5 G4 X7 g5 H/ J
    
19. while X<=Xn-h% c3 j) u. `) T% z
    
20.     K1=f(X,Y);
    ! V\" H4 r. `  |+ H' @9 J  ]
21.     K2=f(X+h/2,Y+K1*h/2);
    ! j- Y- \, d, N% `  [7 v2 _
22.     K3=f(X+h/2,Y+K2*h/2);; H  n) C' q  x! u* L$ _
    
23.     K4=f(X+h,Y+K3*h);
    5 g6 ^' O/ Z  k' A3 E
24.     X=X+h;& T+ l, n/ P# A9 m- {
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    $ F# C- B2 c' O* l! G
26.     n=n+1;                                                   %节点序号加1# Z1 a* A& o, A. t+ C. Y  i9 p
    
27. 
    7 V3 z: |( Z+ \' v1 V' y
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);3 `+ T8 r0 K' B9 `2 c
    
29.     plot(X,Y,'o')
    ) }\" O* Z$ A. O8 j' b
30.     hold on
    \" y4 f% l7 w/ ^& \* B, {+ T
31. end

复制代码

1. function z=f(x,y)5 `# _: N! R# @$ {\" \  X
   
2. z=y-2*x/y;

复制代码

9 d6 k6 }9 z' D) k- I; d% d7 v# H