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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
1 ]2 E1 B) p5 E; X& E* S! V6 m. `6 `定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
. ~. V* o& x; b# Z; x5 l1 W# W这个方法的迭代公式如下：
8 v: l: y( X. r! F[k1 = h \cdot f(tn, yn)]
8 \, D( D- C8 ^1 ~* V& L[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
3 U0 ^# K3 Y, h8 |7 G[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
2 }0 M/ w; o' b6 t7 f其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
, i$ k; x9 T& @$ y1 ]+ p" M这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m( [7 `\" J. C\" e9 u8 k
   
2. 
   \" S, r$ c  L) n. v. p\" ^! ^\" Z
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   : p: z3 ?. u: g+ N7 y+ B
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   ; @+ P# r8 |1 O; J6 m
5.    disp('function z=f(x,y)');4 e. t1 _3 w\" W) t( Q. b7 r  O5 G
   
6.    disp('z=y-2*x/y;');& j7 s- {9 h6 C0 j\" Z. H- B
   
7.    disp('并将该文件保存在work文件夹下');
   ) k5 }5 z% N: X# m
8. end
   $ r' j% y3 Y, z& _6 c* w
9. 
   5 o9 L\" ]8 I, t- |* k! A
10. X1=input('请输入求解区间的左端点X1=');
    . @! X* ]# T8 {/ t7 D
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    ; \( Z2 l4 q/ `% V7 I# q  H
12. Xn=input('请输入求解区间的右端点Xn=');8 _. x1 f% [9 Z$ d+ i
    
13. h=input('请输入求解步长h=');
    / _9 E5 t4 i% }) m7 O9 o
14. 4 U+ O) h$ N; G( o
    
15. X=X1;
    0 P/ B* f9 b% |
16. Y=Y1;                                                        %运算初始点
    6 E. y7 u4 _; g; b9 H! t' b6 o$ M  p9 c
17. n=0;                                                         %节点序号变量置零
    $ N7 V7 Y) f' u# [0 m
18. 
    \" L9 U  v* m( b3 b/ n; h) N
19. while X<=Xn-h* V$ t# W6 }( r0 s
    
20.     K1=f(X,Y);
    ' }; G- r& Y5 j, Z5 F
21.     K2=f(X+h/2,Y+K1*h/2);\" u8 V, U6 k9 H  V$ L2 Y& v5 @+ e
    
22.     K3=f(X+h/2,Y+K2*h/2);/ n, \2 e+ j$ @$ \
    
23.     K4=f(X+h,Y+K3*h);& p* D& @9 |( w/ u+ b; b3 p7 M
    
24.     X=X+h;
    ! P2 o8 J; s7 E
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    $ V; y( `) X6 d+ v. x\" P
26.     n=n+1;                                                   %节点序号加1
    ( B' l: p9 Y1 |' T
27. 
    ! r/ e- I. k# Y0 y* R/ L% @
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);! [( H% M+ _7 I- |. [
    
29.     plot(X,Y,'o')
    1 h' K6 j  I4 v
30.     hold on8 I7 I3 C1 H4 |/ {% R8 v' V9 |/ m
    
31. end

复制代码

1. function z=f(x,y)
   1 `9 C/ F: R; Z  W
2. z=y-2*x/y;

复制代码

! f; Z  |* e1 r+ U