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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
3 m6 e: x9 H% K- _[\frac{dy}{dt} = f(t, y)]
/ ~. F$ b! g! c& m# w这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
2 j( A# S6 j# Z* {* j[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
3 T  _$ W; C, J- T[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
  A, }1 {% I, Q  ~' h/ S[k4 = h \cdot f(tn + h, yn + k_3)]
: e" v: s7 i' j[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m. I, Z& J/ B; [
   
2. 
   5 L, E9 D, b6 D8 n) m
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   \" H\" Q% \1 K\" ~4 l; P. x# f3 @
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');& Y- \1 ?) D3 U1 _% C, R  r
   
5.    disp('function z=f(x,y)');
   7 `1 j* Q\" |7 l9 p9 b
6.    disp('z=y-2*x/y;');
   / ^& I: g% v! {& Y. V% U) N( C2 L) z
7.    disp('并将该文件保存在work文件夹下');
   ; u\" |7 @8 |* f+ o/ j& Z
8. end
   ) o* i3 @! Z4 d+ Z* g$ z2 u
9. & ?. Y* j2 O$ ]9 C$ x* @/ b) ?
   
10. X1=input('请输入求解区间的左端点X1=');* \5 c) c( c) h6 m/ b3 O$ _
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');: E9 n2 U0 V\" |
    
12. Xn=input('请输入求解区间的右端点Xn=');6 ?3 m& n# w$ g\" X, S9 i
    
13. h=input('请输入求解步长h=');$ }; ^& h* U\" _0 Q, S3 \4 Q  |
    
14. 5 D) c* }, L$ ^/ {0 o
    
15. X=X1;6 U( V9 `$ e  |5 t* v
    
16. Y=Y1;                                                        %运算初始点
    & w9 ~- ^9 L$ V( f8 C  h
17. n=0;                                                         %节点序号变量置零
    , t\" H3 A. y% F% U% g' O
18. - h& y- Y, T# p$ E0 [3 ~/ [
    
19. while X<=Xn-h
    , g7 s7 f5 t$ v6 F; D$ L0 o
20.     K1=f(X,Y);1 o: e8 Y% d0 i  T' j7 Y. _
    
21.     K2=f(X+h/2,Y+K1*h/2);
    / ^, r; B3 j( H  n) y0 d
22.     K3=f(X+h/2,Y+K2*h/2);
    / a8 `4 N9 N( b9 }
23.     K4=f(X+h,Y+K3*h);8 `* v+ f$ x. ~5 H% ]
    
24.     X=X+h;
    ' A6 N0 j/ _\" H) q4 J. Q
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    ( T  S8 a$ W  I4 ]8 h' h) ]9 B
26.     n=n+1;                                                   %节点序号加1: d' _7 z) b1 j% _+ ?! J9 k
    
27. 8 C- f\" q( g$ V( m% B( W* Y8 s
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);* N\" l\" u% g/ t4 p
    
29.     plot(X,Y,'o')
    / z6 C& r% K. c8 H1 V9 h% L. t; ~
30.     hold on
    * y0 H# c; A\" v5 I7 p
31. end

复制代码

1. function z=f(x,y)( S: P; |( N* e\" r
   
2. z=y-2*x/y;

复制代码

4 z3 Z% l. F2 |1 ~; M/ g' I