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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
! ~& w3 f8 ~5 G$ d  f这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
2 h- L- ~4 n" ~[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
5 t$ p" D9 E; b6 o. [$ R4 g1 z[k4 = h \cdot f(tn + h, yn + k_3)]
+ [; A% N- t0 R9 C) t6 v# m# }[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
   / R* G7 p\" X. \' j; X\" x7 f
2. 
   / ^) [, J7 m7 ], G
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   9 l' Z0 P# C: S9 y4 b( K$ R2 t
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   3 t5 M. r- F/ D7 p
5.    disp('function z=f(x,y)');- n1 ]6 f3 E3 ?: {3 K/ @
   
6.    disp('z=y-2*x/y;');
   ! w: P* B1 R4 l& m$ u( {( Q, h( K
7.    disp('并将该文件保存在work文件夹下');8 ~1 }& u0 J( M+ R
   
8. end
   ! e2 E+ S% F/ U
9. 
   ( R  S) j0 K5 i1 z& F/ |1 E9 X
10. X1=input('请输入求解区间的左端点X1=');
    & |0 J8 y$ o. e( r2 V5 D& j
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    ) X8 H: K. g- f: @: r& s) O
12. Xn=input('请输入求解区间的右端点Xn=');
    . b; N' o1 h. y- a
13. h=input('请输入求解步长h=');
    \" a+ w% c3 n+ v
14. 6 v4 M) k\" N: x6 Q
    
15. X=X1;6 d- I1 z* \5 A2 a$ m! _
    
16. Y=Y1;                                                        %运算初始点
    6 ?' t) ]! o\" E- ~! u
17. n=0;                                                         %节点序号变量置零4 N9 m( i$ S9 ]
    
18. * V! J, e& P! ?! `3 v6 q  P8 H
    
19. while X<=Xn-h: k) o! a$ Y7 D% H0 b6 x% R
    
20.     K1=f(X,Y);- S1 N! g$ V$ q) H6 Q: B! A+ _  [
    
21.     K2=f(X+h/2,Y+K1*h/2);6 I- D\" Y& R) Q8 I$ V8 ~! [, k
    
22.     K3=f(X+h/2,Y+K2*h/2);! R' o* P! _& G: H4 l
    
23.     K4=f(X+h,Y+K3*h);
    1 g2 h4 T  V8 |# L0 m$ g
24.     X=X+h;% K- O, F( V) Z3 g$ D. i$ ^
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式6 C6 {8 Q- {. J% k
    
26.     n=n+1;                                                   %节点序号加12 J+ r! V( e% m
    
27. 4 H( P- J$ h: A1 W
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    4 `/ ~: ~' h! ^: ~: |/ W
29.     plot(X,Y,'o')3 M$ r* I0 S8 v) Y7 Z
    
30.     hold on
    7 A, x8 K- T3 p  f6 ^6 \; M
31. end

复制代码

1. function z=f(x,y)8 I& g1 K* L+ j
   
2. z=y-2*x/y;

复制代码