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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
/ ~' P- Q  Y& F定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
) k' `! S2 R9 Y: C7 ?" H& D[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
8 b7 W  `3 O+ ~. s! w1 Z, ]* E[k4 = h \cdot f(tn + h, yn + k_3)]
! {6 X+ ?6 C' S% v, K! B[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   : Q4 _4 p8 @2 q/ p8 o
2. 7 i3 {  T1 W\" \: X) a# r0 Y
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   ( i( X2 F; E- g) m
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   3 ]5 n& k! `7 |# G2 D+ ^9 o
5.    disp('function z=f(x,y)');: O7 W+ W$ p5 q; m& c) W6 Q2 }
   
6.    disp('z=y-2*x/y;');
   ! p. N1 ]: O& S4 r# p- i$ F
7.    disp('并将该文件保存在work文件夹下');
   & G1 F2 ^# y\" B' ]5 t, ~- Q' C* I; J
8. end
   6 `5 y, s; g. w4 a$ m. N, B7 E1 j
9. ) j6 C) V+ W' E! W/ Z
   
10. X1=input('请输入求解区间的左端点X1=');
    . \' c2 |4 j/ X
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');' l; }, K8 A4 ], r* W' t
    
12. Xn=input('请输入求解区间的右端点Xn=');' y  [% I9 L; m/ L# T7 z
    
13. h=input('请输入求解步长h=');) Y! Q& d4 a3 f) e
    
14. 6 t$ @# ~% B7 K' S% I1 b
    
15. X=X1;0 V0 R5 b1 x4 ~\" [* y+ ]. t* g; c' D
    
16. Y=Y1;                                                        %运算初始点, d3 c3 E* _) u$ c7 I
    
17. n=0;                                                         %节点序号变量置零
    ; u4 F6 T8 o2 c\" K. ?, G* z  B$ V
18. 7 G2 d! w1 Z\" |8 K* D* n# q# H
    
19. while X<=Xn-h\" p* o2 f; ]1 o# @7 O
    
20.     K1=f(X,Y);
    / \9 x4 C# {, o  w  I; V5 w
21.     K2=f(X+h/2,Y+K1*h/2);
    # [9 d$ |6 C' \, @: ]: c3 H
22.     K3=f(X+h/2,Y+K2*h/2);  C9 X: I6 o1 u2 q2 s( `# ?% ~& D  Y
    
23.     K4=f(X+h,Y+K3*h);
    ; T+ m1 }. g+ q. }8 v0 K( v3 S
24.     X=X+h;
    $ }1 J- r5 l, z3 f
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    7 n3 ~, E' T% O* I5 H
26.     n=n+1;                                                   %节点序号加1
    5 c  E! W# G* i3 @) L% S
27.   O0 D* W- V% C/ t0 c% c/ j
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    \" ~$ e1 K4 u, q& I5 G
29.     plot(X,Y,'o')  }# f' h0 G' u/ F: o8 W9 r0 w3 N
    
30.     hold on: ^% X+ N5 d' [) j- x$ v) g6 T7 @
    
31. end

复制代码

1. function z=f(x,y)
   \" C/ s  L( {2 O3 p4 b) M
2. z=y-2*x/y;

复制代码

3 M4 D" k) ]$ [9 d$ B- D