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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[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})]
* t7 \- T  o9 Y4 A7 U# T. {; q! ][k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
9 a" X7 T% d2 D& }$ w1 m& D这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m8 ?( \- C2 G! E, T
   
2. % _) ]: C& W2 ?; x
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名). a\" W+ V8 l1 x$ h2 Q7 ~
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   + R! k( h7 I- D  h7 G% t
5.    disp('function z=f(x,y)');
   ( V3 a0 M! R9 b, {
6.    disp('z=y-2*x/y;');$ A5 j$ L3 y1 U
   
7.    disp('并将该文件保存在work文件夹下');. @$ W/ d; t( I6 b! c3 }% ^# W$ e' V8 L
   
8. end
   ; a2 d( [- V9 m7 {- V5 ^
9. 
   0 s( w+ ~% ], A5 S1 R3 Q
10. X1=input('请输入求解区间的左端点X1=');. q# g$ h, _' ~( ^9 m; `
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    ' V' v  i4 r  x9 ^1 g( p) U
12. Xn=input('请输入求解区间的右端点Xn=');
    0 E( u4 V0 ^\" b
13. h=input('请输入求解步长h=');1 p4 I; J5 a  G
    
14. 
    : n* A) ]- s7 r
15. X=X1;+ s% o9 b- T; I
    
16. Y=Y1;                                                        %运算初始点
    2 J7 P8 E3 I9 {7 [
17. n=0;                                                         %节点序号变量置零$ H/ w, A; U* K) z0 Y+ f+ \
    
18. 
    6 F1 V% b  D- H  k
19. while X<=Xn-h
    , q( ^\" Z9 B7 ]: z\" h5 q! |& a- _
20.     K1=f(X,Y);5 l2 x0 c6 j4 _; [6 w+ ~
    
21.     K2=f(X+h/2,Y+K1*h/2);
    7 j) t3 B; [\" d6 j
22.     K3=f(X+h/2,Y+K2*h/2);2 H8 }. }+ H* t; `5 G
    
23.     K4=f(X+h,Y+K3*h);! x\" Z; U* H0 C  p; N
    
24.     X=X+h;( ^/ p2 i+ [$ M
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式0 w: K) d7 t6 h+ f7 d! I
    
26.     n=n+1;                                                   %节点序号加1& w. y# S1 r' q/ ?  _# a9 N& ?
    
27. 
    7 q7 o$ n* V( t6 n# d
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    - d3 w+ r+ d3 }6 g2 s) {% L
29.     plot(X,Y,'o')# p\" \: C# k, D/ ?$ x1 y- z\" Z
    
30.     hold on8 k5 b( e, D. l7 b1 ?4 K
    
31. end

复制代码

1. function z=f(x,y)
   0 Q: M# i  l$ J7 j% s7 X
2. z=y-2*x/y;

复制代码

* R+ b  v2 |( I( i) q' o7 Z