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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
7 `) y6 \# r5 K8 w2 T5 e[\frac{dy}{dt} = f(t, y)]
+ K/ b. v) p- Y1 u. P9 H6 S* A这个方法的迭代公式如下：
+ B6 Q6 k3 ?, S2 z0 A, }[k1 = h \cdot f(tn, yn)]
1 ]: X; t1 ~; z5 T9 z0 L3 m& Y- R[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
6 }5 G9 d' H3 D/ V[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
# l- f1 p0 A& ~[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
8 ?2 G. d4 Y2 P; s& h+ ^9 A- K其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
* N  x: I+ v( C) g这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   ; ~  d; V) V/ j9 A) a+ u
2. 
   0 R1 y* ~* p& ^( j- o( ?6 z
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名). |, y5 y5 Q3 F4 ?& Q# c
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');: j; i! i9 A( x2 V1 T2 I
   
5.    disp('function z=f(x,y)');
   ! A- l) }\" \0 W) k5 I2 u' a# }
6.    disp('z=y-2*x/y;');) ~* ^  V0 u  N7 a3 g7 p& J
   
7.    disp('并将该文件保存在work文件夹下');/ ?) w% S$ V. D
   
8. end ! M& q* p) w5 P
   
9. 
   * ^$ Z/ y4 o4 o9 x' x. u
10. X1=input('请输入求解区间的左端点X1=');
    , E( k- w2 d& j
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');. }; O$ }* v: j- I1 ~* g: m2 v9 L; ]
    
12. Xn=input('请输入求解区间的右端点Xn=');+ C. Z9 Y) \' e4 K
    
13. h=input('请输入求解步长h=');# h7 |  b4 I/ z& Q4 |
    
14. 
    \" ?, b0 L/ P& B8 i2 Y' `% q# H
15. X=X1;
    9 @2 D0 E; |9 M+ o, A
16. Y=Y1;                                                        %运算初始点
    0 U& `  w) X$ z' z+ D4 s( P) G( z# {
17. n=0;                                                         %节点序号变量置零
    % ?2 D3 K: Q  w2 [4 x1 P
18. 
    $ D+ J7 F1 U7 Y/ z2 H1 b7 N
19. while X<=Xn-h) Z0 t$ B; L, C* `
    
20.     K1=f(X,Y);
    * G( S6 n* C/ e& B
21.     K2=f(X+h/2,Y+K1*h/2);- n4 N+ p: Y3 ?, `; x6 v: X5 _3 x
    
22.     K3=f(X+h/2,Y+K2*h/2);
    . b, {; G& S$ A' ?\" D) K\" y
23.     K4=f(X+h,Y+K3*h);
    4 O; g- X+ K; W/ x' _0 ]
24.     X=X+h;3 N0 Y$ I8 ]3 Q3 l9 V, Q) f; C1 O
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式4 F6 K4 k7 b7 y0 O$ D2 Y
    
26.     n=n+1;                                                   %节点序号加1
    ; d- s* C) a3 r3 x0 i1 d
27. * M, n2 ~1 m2 y1 {3 j$ I
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);+ h. T: g; l7 h1 K$ @0 m$ k3 t
    
29.     plot(X,Y,'o')4 o3 s3 x5 N3 N' K/ N\" ^6 V  y
    
30.     hold on
    7 h& j9 j9 i( v8 j
31. end

复制代码

1. function z=f(x,y)
   5 e- K7 k8 k) f/ E& Q, F0 F& p
2. z=y-2*x/y;

复制代码