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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
1 J+ s: N; Q* r% E2 X' M: l[\frac{dy}{dt} = f(t, y)]
2 h9 d' w1 ^$ b6 l3 {2 \- a这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
# `; c" N; a8 F  p[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)]
. ~  b' v0 m" t% {0 v# l其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m  ]  V0 v7 z6 V4 z. ~1 Q
   
2. 
   9 V( X- q/ \\" ?5 W\" S6 B1 A
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)8 @7 j6 l3 E) ^- M
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   \" r( ?' K$ |- F' B0 I  T
5.    disp('function z=f(x,y)');
   6 F7 P\" I6 A; H2 D& N3 o
6.    disp('z=y-2*x/y;');. u% {$ @: V/ E& V/ y1 N- X
   
7.    disp('并将该文件保存在work文件夹下');
   . e3 `  B( I' K+ n0 y
8. end 3 {. s; K- i' j1 C. q$ o) u
   
9. 
   9 C% ^0 Q/ {* _/ N; `
10. X1=input('请输入求解区间的左端点X1=');
    1 z9 K. p  ?$ e0 K
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    / _- a0 B\" n1 D4 u\" I
12. Xn=input('请输入求解区间的右端点Xn=');' _8 R, q' j+ |
    
13. h=input('请输入求解步长h=');
    . i1 s# X/ W) z3 V1 x* I! C
14. 
    8 B* \$ e0 b5 R0 {! L- c8 a9 J
15. X=X1;
    4 M8 n+ U8 V3 }) V: H- O
16. Y=Y1;                                                        %运算初始点
    / x( Y1 @0 r% z/ b- `' ]4 }
17. n=0;                                                         %节点序号变量置零
    , S\" C& F) U\" ^& [7 s+ Z9 f
18. 
    1 ?5 ?' _% j# r2 a
19. while X<=Xn-h& s* f/ H, J9 T. Y+ [
    
20.     K1=f(X,Y);: J* g4 a1 s1 }! Z\" D6 c
    
21.     K2=f(X+h/2,Y+K1*h/2);
    - g; F; ^  `1 ?' u$ u4 \
22.     K3=f(X+h/2,Y+K2*h/2);( f6 ?2 G  ?/ H) ]5 x! s8 x$ S9 _
    
23.     K4=f(X+h,Y+K3*h);
    / n7 B2 S+ L' c3 |3 a$ X3 s
24.     X=X+h;
    * C3 G+ K1 [% J
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式; O$ A- Y9 I7 O\" S4 j
    
26.     n=n+1;                                                   %节点序号加14 ^% p: E, j5 Z1 m7 B9 X: v- j
    
27. & p. ]: b\" T- H$ f# C
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);/ b4 |  d# n# Y
    
29.     plot(X,Y,'o')7 r1 y. k! w6 k7 A
    
30.     hold on
    . ~\" Q0 I% F7 n, |2 |. }
31. end

复制代码

1. function z=f(x,y)
   ! G7 j! C, z6 w0 E' @! n
2. z=y-2*x/y;

复制代码