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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
0 i7 c! I) g; Y这个方法的迭代公式如下：
6 j6 I6 L  ^( u* F4 |[k1 = h \cdot f(tn, yn)]
$ S+ C* b" ^; X5 b  v[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[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)]
5 u% x8 G: ]/ n9 `其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m- [, X# p6 x  z
   
2. - n9 M7 P\" A$ n
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名); t\" o* o2 N, ?- f6 b; H
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   * a/ e+ Y3 [5 \  ?
5.    disp('function z=f(x,y)');
   ( x\" U( }8 X* R) d
6.    disp('z=y-2*x/y;');
   - {9 M0 h. c0 d+ G% ^- f
7.    disp('并将该文件保存在work文件夹下');\" s+ [: m7 ?* A
   
8. end
   ) U, q$ S0 \/ U& n; E, i
9. 2 Z$ ?- C7 u5 B5 z\" c
   
10. X1=input('请输入求解区间的左端点X1=');# A& @4 Z( F, T( q7 q
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');( z( d3 A4 i  {0 s, U1 j
    
12. Xn=input('请输入求解区间的右端点Xn=');
    : |; f0 \2 d# w- [+ m% U/ h! b
13. h=input('请输入求解步长h=');8 b& N* X1 z6 x
    
14. 
    , T( o7 I/ ?- h4 v\" q5 ?  V
15. X=X1;) k( R) P2 P8 _! Y
    
16. Y=Y1;                                                        %运算初始点
    # d* U+ Z1 F7 A
17. n=0;                                                         %节点序号变量置零7 _! P) v/ R( C6 L  ~( m- n0 W1 q
    
18. 
    ; N/ h. O8 M2 ^
19. while X<=Xn-h
    ) t! @' M4 v5 z1 J4 }, y8 i9 c5 n
20.     K1=f(X,Y);  m, Q- s! _9 i8 H: Y\" a5 ^\" H
    
21.     K2=f(X+h/2,Y+K1*h/2);6 h. l' W( G$ X( U8 m4 e! v
    
22.     K3=f(X+h/2,Y+K2*h/2);4 }# k2 \2 e2 ?# c1 z
    
23.     K4=f(X+h,Y+K3*h);
    2 v\" w, @1 r8 L9 q( ]( H; N0 A
24.     X=X+h;, D. S% e\" c# [8 X
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式/ D\" t1 y- I3 d# G* ^
    
26.     n=n+1;                                                   %节点序号加10 G, g8 K: W4 D
    
27. $ }1 n& S\" @& o1 g( Q3 D: D
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    : ?% d: x0 ^0 E7 a1 _, a( j: P0 D
29.     plot(X,Y,'o')
    4 j+ A+ P+ `* ~\" ~
30.     hold on
    ; v- B2 v# W' j) ?' J
31. end

复制代码

1. function z=f(x,y)
   ; f& v/ c9 }' t4 Q: W6 f# U
2. z=y-2*x/y;

复制代码