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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
( b  b. _3 `, v, w2 B# D定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
" \' L: H7 H& c7 p# X7 ~这个方法的迭代公式如下：
. b7 {7 q! f- u% z- _3 O) W[k1 = h \cdot f(tn, yn)]
$ E5 Q; d) i/ J" R  e0 C, s9 J[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
9 _1 U7 l) `' t8 B; [( b% k[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
) _; g. \7 h" B# \3 c7 s2 T其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
" o3 e+ _/ c2 T  x8 z/ M# M这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m! U: B4 O( g2 w+ |
   
2. 
   ! V7 b- W' D6 f8 }' [6 u
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   * E$ w8 j) E: H' D' W. G9 W
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   - p- B. w! h5 D
5.    disp('function z=f(x,y)');: u! [; a& a$ O. W0 n# U
   
6.    disp('z=y-2*x/y;');
   9 Z  Y\" q$ t1 Z. A
7.    disp('并将该文件保存在work文件夹下');
   & H/ R% \7 I! l: ~2 z2 ]4 x
8. end 3 \4 _2 Q5 |' n/ p
   
9. * d5 [' G+ g: x& Z! Q$ V\" y
   
10. X1=input('请输入求解区间的左端点X1=');
    6 f7 O! m# t' V- A# s! k1 [
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
      K8 X8 G0 N8 W( f' I! D: w
12. Xn=input('请输入求解区间的右端点Xn=');
    5 r1 p\" N, f( Y\" K5 ]
13. h=input('请输入求解步长h=');5 i& U- p6 M# t) u  R+ _
    
14. 
    % s/ |6 t. H% v. q, r
15. X=X1;3 `2 b9 h7 v; S- t\" g5 X5 ]. F
    
16. Y=Y1;                                                        %运算初始点
    3 [& e. o8 M$ Z/ N# u2 L
17. n=0;                                                         %节点序号变量置零: z4 H4 m' H7 }1 M\" F
    
18. 
    ) _\" H3 P! i1 y+ m
19. while X<=Xn-h) h3 B9 S' n, U% R2 U
    
20.     K1=f(X,Y);' U4 K5 S- O9 l/ Y/ p: x
    
21.     K2=f(X+h/2,Y+K1*h/2);
    ' R* \3 x, Y' y  z
22.     K3=f(X+h/2,Y+K2*h/2);% c2 e* v; ?! Q% y
    
23.     K4=f(X+h,Y+K3*h);% K% V& f) b9 U' W; t
    
24.     X=X+h;! I. u+ t5 ~7 a( E6 T
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    ( s# d3 D/ |\" L, j
26.     n=n+1;                                                   %节点序号加1/ K, ]6 }( [; U/ b8 {
    
27. / U* }/ Q. f6 |  z
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    $ H) ?0 t0 D& ]% M- X
29.     plot(X,Y,'o')* C! D2 S& v7 U3 a. k2 l
    
30.     hold on
    $ b8 F\" p  ?5 R: d  Y! B% c1 X, N& s
31. end

复制代码

1. function z=f(x,y)8 ~: s/ Y, z' y
   
2. z=y-2*x/y;

复制代码