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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
% G7 F0 j& Q# \1 Y, ^定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
4 @" X8 {# }. X4 b& O  I6 f0 d* F[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
& I4 E1 ]5 E: w[k4 = h \cdot f(tn + h, yn + k_3)]
0 G: P: r  d- L[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   & @1 W3 l+ G5 l- Q0 G
2. 
   ( |4 B7 M, u# x$ D
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   ' n- }6 k( _* B1 i
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');  u7 {\" W; y4 W; Y; X5 y' {
   
5.    disp('function z=f(x,y)');  }6 {- e9 Q; d2 H3 Y3 A
   
6.    disp('z=y-2*x/y;');/ Z( Z  e1 M/ a% k; L5 J1 _8 r3 P& |
   
7.    disp('并将该文件保存在work文件夹下');2 l( P$ o8 @% J9 W  F( {4 f# E
   
8. end % B+ J0 f% @$ Y+ d: O6 A+ ?4 D
   
9. 
   . t& d; \# r: |) M) K
10. X1=input('请输入求解区间的左端点X1=');% R. b2 C8 i1 @$ k\" T' G
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');+ L4 Q0 e, v. K& [- I
    
12. Xn=input('请输入求解区间的右端点Xn=');
    7 K! ^$ g3 t2 [8 Z3 w
13. h=input('请输入求解步长h=');! I) N1 e7 m& _! c: }# O/ h- N' o
    
14. 0 g/ e, @( L( \$ H6 q( a/ Q5 e
    
15. X=X1;0 J# w! U* F2 Z( t7 ]3 H
    
16. Y=Y1;                                                        %运算初始点
    3 E' e6 E8 l% z  Q3 k  M6 F
17. n=0;                                                         %节点序号变量置零/ l# X+ V1 O4 o# ]
    
18. 
    ; [8 ~; I& _$ Y4 c6 x
19. while X<=Xn-h
      h5 o2 g( J& s7 m1 I
20.     K1=f(X,Y);
    # ?: R4 x5 H$ o( q1 c+ ~5 O
21.     K2=f(X+h/2,Y+K1*h/2);
    # P- c. `& J8 Z
22.     K3=f(X+h/2,Y+K2*h/2);3 W' ?- ?7 T7 A9 |5 l0 q
    
23.     K4=f(X+h,Y+K3*h);
    : q- x0 Q5 U4 e6 y5 V4 k0 g& _
24.     X=X+h;& F0 K/ E4 I7 V$ m7 H8 d
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式, a7 m$ w6 T5 x; @* M5 A0 C
    
26.     n=n+1;                                                   %节点序号加18 r% i+ @! c' F( ]
    
27. 
    # r9 @7 ]8 h( ]  g! s
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    5 z5 W3 c2 ^; _. e0 D: p7 F) x2 q1 E
29.     plot(X,Y,'o'). n5 h) p6 {9 F/ ?. X
    
30.     hold on
    0 e6 M, c# K& k
31. end

复制代码

1. function z=f(x,y)
   $ B; Z( w6 G* \5 ]+ ^% r
2. z=y-2*x/y;

复制代码

$ R) {, o0 B& w2 E