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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
" @( I7 J& R) T4 ~# Q& t[k1 = h \cdot f(tn, yn)]
2 q+ @4 }: o9 N[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
2 s( A$ J. Q7 D+ G4 V  d[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& H* ]; h! }, U. f' @/ o其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m! V( G1 `$ e' P8 d8 H1 R
   
2.   L6 I6 S6 w+ }$ F, m9 U
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   , O) a# _\" a* l\" N# e
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   ; Y( c2 `$ M8 C' E\" u/ x! J
5.    disp('function z=f(x,y)');
   # h1 p  j; e( e! w9 A/ L/ a. z# B
6.    disp('z=y-2*x/y;');
   / H+ P3 r/ d7 B. Z* O9 R
7.    disp('并将该文件保存在work文件夹下');# w! V0 B( L3 {7 _0 {
   
8. end \" `3 _0 S* N- t& K3 C; x
   
9. , U+ {0 F+ Y( `! o
   
10. X1=input('请输入求解区间的左端点X1=');9 h7 [2 c* u; i, k
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');2 L2 @& v% `* v# ?5 X7 d
    
12. Xn=input('请输入求解区间的右端点Xn=');
    6 [+ F/ s- J. y# _4 t/ {
13. h=input('请输入求解步长h=');
    ; o% b1 C' I( f: y& E9 Z
14. 
    9 ^0 ~% U* [9 B9 ]4 E3 K8 s
15. X=X1;
    % y7 a' M* Y, q0 ^
16. Y=Y1;                                                        %运算初始点
    / g9 {6 q( D6 v) `/ R$ e
17. n=0;                                                         %节点序号变量置零
    5 T; T: F; g$ i
18. 
    . y% E: D1 [* I' H
19. while X<=Xn-h
    4 U8 }  c- k# Q  N6 J\" V, Y
20.     K1=f(X,Y);
    ) C8 @1 D: G5 ~' U  F7 X+ C
21.     K2=f(X+h/2,Y+K1*h/2);, \2 m* }6 B  O
    
22.     K3=f(X+h/2,Y+K2*h/2);: [/ Y5 I& x- k/ P
    
23.     K4=f(X+h,Y+K3*h);$ q* @# I, s) _; J
    
24.     X=X+h;! p2 w7 ?7 `( Y% }
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    3 \, [% F1 M# l. h, n
26.     n=n+1;                                                   %节点序号加10 e: G! M5 ?9 [
    
27. 
    # X0 o, E- Y% h+ @+ u1 \+ y
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);% Y/ k7 {, l& K+ ~4 V  _2 G
    
29.     plot(X,Y,'o')  b1 g/ p! g3 e/ g\" ]
    
30.     hold on
    , E& E5 D. ~# }1 J2 W/ V
31. end

复制代码

1. function z=f(x,y)4 n; C! p$ R0 t) H' N
   
2. z=y-2*x/y;

复制代码

+ H# p) u" [7 m0 q