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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
" o! f/ q1 Y& I& `8 P[\frac{dy}{dt} = f(t, y)]
/ ?1 |7 \2 D5 S/ P6 p这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
$ D3 z3 ~) L1 ^, O* U$ o; J" A$ |% {! d[k4 = h \cdot f(tn + h, yn + k_3)]
4 N- u3 n4 p) {9 i( G$ Z[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
' r1 d  I0 r' V* E其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   4 n: @) H/ @- Z: c. p1 f8 ~  z
2. / \; u1 q0 q! ^% h1 R- b7 T6 P
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)& I4 m6 b, C3 {- a  s! m4 D& f
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');! X4 H) w& G8 g7 l0 z& y6 U
   
5.    disp('function z=f(x,y)');6 W\" \7 c4 e5 L  L; U0 O
   
6.    disp('z=y-2*x/y;');# I4 }! j9 C6 b) z( ?# K
   
7.    disp('并将该文件保存在work文件夹下');0 o4 S8 v1 n3 i2 [
   
8. end / q0 k- N( k, p/ \0 F
   
9. 6 o1 x' W7 p& @5 C  [: }: k- _
   
10. X1=input('请输入求解区间的左端点X1=');
    : X/ x2 V- {& U3 S
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');0 i3 O1 ^+ ]1 E6 e% I8 W: h, q8 _% D
    
12. Xn=input('请输入求解区间的右端点Xn=');. M. c) W# h3 r5 ?  g& q3 A
    
13. h=input('请输入求解步长h=');
    * w5 @\" j6 N5 e% r( Y8 c
14. & M# |; F- K7 D9 ~9 i, u6 k
    
15. X=X1;
    1 g0 d8 F2 z% A
16. Y=Y1;                                                        %运算初始点
    + [\" d. j6 _9 y\" [( H$ T+ P8 V
17. n=0;                                                         %节点序号变量置零4 x3 ~3 d) H% G, w
    
18. / ^, q1 y/ o. Y3 g
    
19. while X<=Xn-h$ X% A% f- A$ f7 A
    
20.     K1=f(X,Y);2 e! @* Z0 o; [! R3 {+ e5 @. P; D
    
21.     K2=f(X+h/2,Y+K1*h/2);
    8 k5 W2 ?9 K, H1 }9 [0 k
22.     K3=f(X+h/2,Y+K2*h/2);/ B3 g3 X4 M% u1 K, P2 |. p
    
23.     K4=f(X+h,Y+K3*h);( Z$ \7 o9 G& g$ W
    
24.     X=X+h;
    4 T' H( h( g$ `7 V
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式( R1 i0 H7 h4 g( E1 F$ F
    
26.     n=n+1;                                                   %节点序号加1. F3 v1 S, c% Q+ [. ?3 E
    
27. 
    , a+ f+ ~9 G\" D. U
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);9 j! r7 q# [$ N% Q
    
29.     plot(X,Y,'o')
    / v7 k% e3 E) J
30.     hold on5 w% j% e: u0 x+ T/ C1 I
    
31. end

复制代码

1. function z=f(x,y)
   2 F1 `( R* [% o1 m3 w% i8 m
2. z=y-2*x/y;

复制代码

8 x( v. n' {, N/ f