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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
) ?' R4 P8 ^* r; D( N5 V, p[\frac{dy}{dt} = f(t, y)]
. V8 N% X8 }0 L/ T8 H) |- S" o; s1 J这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
* r; n. w0 _. G. J3 k6 M3 v[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)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m4 m0 |* [  e; P5 \/ O
   
2. 
   , M& a* Q- B; v- K( k4 z
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)$ U/ V  w# c- U
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');; i, {7 ]/ ^* U9 C
   
5.    disp('function z=f(x,y)');
   0 [$ P! n/ }+ ], @/ _$ s
6.    disp('z=y-2*x/y;');3 [% v5 r0 h# L
   
7.    disp('并将该文件保存在work文件夹下');
     ~4 |$ u- X* _7 p+ d) i+ X
8. end
   0 j; H8 b1 x* R& T$ ~2 V; [
9. 
   . q+ c2 }# T0 W0 s: X\" [5 s) P6 r
10. X1=input('请输入求解区间的左端点X1=');/ b9 j6 c( k0 }/ U  Z7 \- d' U
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    ' p. _- A& n8 N
12. Xn=input('请输入求解区间的右端点Xn=');, W! v) v; Q$ P
    
13. h=input('请输入求解步长h=');
    . _- e; P/ Z% [) E) P
14. 
      j2 \7 n% w6 p
15. X=X1;
    9 o4 x- S4 ^0 J
16. Y=Y1;                                                        %运算初始点; [' _( Z  J1 D+ @; V# K4 |) t
    
17. n=0;                                                         %节点序号变量置零5 b+ g9 F. \% ^6 a/ m( a1 w
    
18. 
    $ I3 H5 C( b% w/ o. |4 t1 W
19. while X<=Xn-h
    0 `1 a1 V( k/ q8 ]* @( k1 V9 [
20.     K1=f(X,Y);0 P) E- ]$ X4 w0 l- y& d* O) a5 `- P% }
    
21.     K2=f(X+h/2,Y+K1*h/2);
    ) }, N4 q7 d: n1 N1 j
22.     K3=f(X+h/2,Y+K2*h/2);+ ?\" x% w( J( q& r
    
23.     K4=f(X+h,Y+K3*h);
    8 C$ e( p4 ~$ M$ J% O/ ~5 r) t  _; F; c
24.     X=X+h;! |9 G; B4 z  i; h1 z: A
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    3 T% J% o% V1 [0 m3 `
26.     n=n+1;                                                   %节点序号加1( _+ O& s, W% h; ~
    
27. 4 d4 Z\" I% u! v\" M
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);; g2 S3 F/ e! z7 Y
    
29.     plot(X,Y,'o')6 h$ o. K1 I4 n3 ^! L: S
    
30.     hold on
    # J: Z/ i! {. i9 \
31. end

复制代码

1. function z=f(x,y)
   # ^& t$ G' b7 Y9 a  C9 n4 x
2. z=y-2*x/y;

复制代码