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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
" L+ m1 k/ M1 I. @& h$ N. b定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
% Z/ a( D1 T( B[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
6 O; [) e; [( l2 P" C" `& |[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})]
[k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
" ~  B$ d8 X4 E其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   1 i# u, ~7 m& I- S! |
2. 
   . W4 }( A! S# I6 [/ Y* r# U4 Q$ q
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   8 K2 Q3 D- M5 v+ \( {& m
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');5 y: D: F  b* A1 h7 E0 Y8 I
   
5.    disp('function z=f(x,y)');
   & Q( r+ R\" i' r  J% t
6.    disp('z=y-2*x/y;');( D3 q. @% [; l6 y  T+ }' k
   
7.    disp('并将该文件保存在work文件夹下');
   ) F8 J+ U, {  S; [- T, B
8. end + S* Q; ^* n+ X2 i
   
9. 7 f- ~5 p! L# V3 C2 n8 T* ~
   
10. X1=input('请输入求解区间的左端点X1=');
    - ?, \  V$ ~, v
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    2 b' i8 v\" i2 H2 v! E3 p9 Q' K# W% C
12. Xn=input('请输入求解区间的右端点Xn=');
    / v2 M) C0 w3 n\" |\" x0 C0 Z
13. h=input('请输入求解步长h=');6 ^6 n5 ~2 ~6 z9 N) Q, _  A\" M% \( e3 }
    
14. 
    \" H3 _! q7 G\" @2 H
15. X=X1;
    5 T' k) W. H4 U6 H/ T' B, ^
16. Y=Y1;                                                        %运算初始点  B7 U2 ~% P4 d5 u
    
17. n=0;                                                         %节点序号变量置零
    + c) o2 J9 o: T) D7 K9 Y
18. 
    % D$ X4 W4 ?% k. F& A1 T( x1 M( i
19. while X<=Xn-h\" Q5 m' f2 H5 w  Z0 [$ _
    
20.     K1=f(X,Y);& c& b' f3 y! x2 d8 N9 c6 r* Q
    
21.     K2=f(X+h/2,Y+K1*h/2);
    ) e! Q. c- {\" Y  Z0 r
22.     K3=f(X+h/2,Y+K2*h/2);& L( n5 D: `' x, i
    
23.     K4=f(X+h,Y+K3*h);
    0 O2 `\" E; c% T/ y8 b- n9 D
24.     X=X+h;! Y& G2 ~2 @9 u' a% x: ]
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    ; G4 x\" S: ]  ?% ~$ n) K+ U
26.     n=n+1;                                                   %节点序号加1
    $ x) h* M+ b# e3 @* T\" @& R1 U
27. 
      h9 x7 E, \- S9 a
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);! ~! T' v1 L. Q% h; l
    
29.     plot(X,Y,'o')
      M: q9 P( W. O, d; ~
30.     hold on
    5 s! c( _/ O3 t& B2 Q1 ~) l  v- ~
31. end

复制代码

1. function z=f(x,y)3 P4 C- G! p9 @0 q$ d! A. c
   
2. z=y-2*x/y;

复制代码

( |& y' h& z2 \6 V