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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
  C3 \& l8 n) k/ d# w[\frac{dy}{dt} = f(t, y)]
8 D5 X5 `6 o- `+ B3 P& D9 Q- o: N1 n7 \这个方法的迭代公式如下：
; @1 ^' O# K( q" D4 n- A[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
0 P* L  x9 c  |[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
$ S; r% `( C) N# N5 h, D[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
  Y, F" {/ ?6 _+ o# @+ J# i$ V. ?其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
8 g3 M! W* x1 M$ o1 w这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   2 a: j4 O4 N- C. H7 K. B9 h
2. * s& Y  b) {' H) ?
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)8 o4 y6 O7 s, k1 K
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');5 B; ~( d9 E: y% X: t, o6 L
   
5.    disp('function z=f(x,y)');; F+ z. v5 F- b7 z+ T$ g' q' L5 x
   
6.    disp('z=y-2*x/y;');
   ) \, W6 Y9 u+ X- Q) V5 O6 ?
7.    disp('并将该文件保存在work文件夹下');
   $ a6 [: D$ w/ q& y
8. end 7 ^  H' B* ?1 _/ l' U. ^$ d8 u
   
9. 0 G, A+ d  G4 D3 D& L5 \8 I\" q) ^
   
10. X1=input('请输入求解区间的左端点X1=');
    ! k; b7 t& D5 J) e; I! @
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    , F! W$ p% S. @9 I0 P' s
12. Xn=input('请输入求解区间的右端点Xn=');
    - X4 B( w5 F( e
13. h=input('请输入求解步长h=');3 u; R( |( G* L% C# E: s5 `/ J
    
14. 
    ' p2 x% [( {/ T3 q, i# w\" N
15. X=X1;' L4 f% e4 p2 C% R
    
16. Y=Y1;                                                        %运算初始点
    ( W. o- ^/ O\" E2 n; c5 ?/ `
17. n=0;                                                         %节点序号变量置零
    1 k% }& }4 r& n% x0 N1 m
18. 
    # D0 z' G! A1 z$ ]( K% s
19. while X<=Xn-h
    1 N4 m: h. H, L1 n3 s* t' Z1 W: h
20.     K1=f(X,Y);
    3 v# F4 c5 v1 W' J! O5 L8 _0 G
21.     K2=f(X+h/2,Y+K1*h/2);
    ; t( V& R, y% T1 A- D6 ^' ~6 u
22.     K3=f(X+h/2,Y+K2*h/2);
    ! l) Z- `2 ]7 h/ v1 L4 b$ m
23.     K4=f(X+h,Y+K3*h);
    8 l7 v  f2 s4 z4 |* N3 [& `
24.     X=X+h;$ H( ~9 n! i0 Z9 d0 n. q# j
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
      z5 |- h) T& O' M2 ^( V
26.     n=n+1;                                                   %节点序号加1
    9 W/ `, A+ l* M3 b
27. 3 `# m! B2 W  u- p7 U) ?& @4 n
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    9 X5 E1 l6 ]! e
29.     plot(X,Y,'o')
      ]( N9 ?! E  _* z
30.     hold on
    7 t; X, r3 `+ ~5 D
31. end

复制代码

1. function z=f(x,y)  D8 u: O# R; }, R
   
2. z=y-2*x/y;

复制代码