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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
" r. I2 A, g$ U0 m* e. j; |9 S) {定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
! Z1 h8 O1 }3 D1 w  N, t+ m0 u这个方法的迭代公式如下：
0 y2 J& ~1 {7 e: q0 Q[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})]
+ {  s2 v2 h+ e+ r7 n0 d) g[k4 = h \cdot f(tn + h, yn + k_3)]
4 l5 U+ q+ j6 E4 A9 ]  y" I1 i[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
' p. f+ M$ w  J- t9 P$ R& ?" `这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m3 ]- ?! c: _* t. Q& x
   
2. 
   0 Z1 U6 f\" F2 c$ \
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   ) P% ~9 M3 d% c- r5 `+ T: w! v. Q) B\" H
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');; r- G* L' F/ [- w0 Y2 s7 u$ E: e
   
5.    disp('function z=f(x,y)');) V' w8 ~# K% Q$ l7 m. k, S
   
6.    disp('z=y-2*x/y;');
   1 g1 f4 v) N4 y$ m! l$ L
7.    disp('并将该文件保存在work文件夹下');
   7 s, S* Q, l/ c- z, }
8. end \" n( y# K0 h& c6 F5 V( X$ p
   
9. 
   - T- M\" `6 }  D/ \( Q2 I5 X
10. X1=input('请输入求解区间的左端点X1=');; _1 u# E. t( n; [
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    $ G: u2 i, u6 W, x
12. Xn=input('请输入求解区间的右端点Xn=');
    6 r. F7 R\" j6 ^, p2 U
13. h=input('请输入求解步长h=');6 X+ b( T5 `\" g
    
14. 
    5 J: I- C4 Q' ]
15. X=X1;* q7 L2 p7 N. N( b' z
    
16. Y=Y1;                                                        %运算初始点
    , N% r. C1 o# j5 N* h* s! H9 |
17. n=0;                                                         %节点序号变量置零4 H7 D5 f2 `/ y* C
    
18. 
    1 y3 l8 C6 o0 k* T; Q
19. while X<=Xn-h
    - N9 ]' k: L% @- U4 w) b# L
20.     K1=f(X,Y);
    . W7 E5 ^: \/ O\" Z
21.     K2=f(X+h/2,Y+K1*h/2);
    * J6 X4 y& g* N8 k4 h3 R\" d
22.     K3=f(X+h/2,Y+K2*h/2);
    6 r* n5 u  a7 J* i5 I' m# V% G
23.     K4=f(X+h,Y+K3*h);  F3 a0 f1 y7 F. ?, W5 ]8 _
    
24.     X=X+h;
    . v0 C  T- g  ~# U
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式6 X* n* ]. U! F) a$ Q; M% b+ T# w
    
26.     n=n+1;                                                   %节点序号加19 K1 M. T7 a' _' K) m# |6 {\" ~
    
27. 
    , l0 ^0 v( O, y! z
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    8 X$ l! E  B. F( p! `
29.     plot(X,Y,'o')6 f- I+ ^+ \5 s\" g# L
    
30.     hold on* s/ w9 w  p, m$ r7 X
    
31. end

复制代码

1. function z=f(x,y)
   - ?6 q: Q& p2 o0 U7 H
2. z=y-2*x/y;

复制代码

; U# |1 p3 N3 |! H, e6 H0 }