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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
9 W9 p/ H0 ?) D' y4 L- k[\frac{dy}{dt} = f(t, y)]
" L+ o3 G7 D$ v这个方法的迭代公式如下：
0 E) G8 z1 Z0 p& a[k1 = h \cdot f(tn, yn)]
* p/ K' G! F7 L& ^7 \[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)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m: v: g  e! r* c6 M4 j
   
2. 
   ) ^+ d3 S9 o2 F$ o7 G) A2 p\" R
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   \" y& J! Y; c6 ^# |( `
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');8 r3 r0 d& l: G1 c2 v
   
5.    disp('function z=f(x,y)');7 ]  q4 s9 x. z( h+ V( A
   
6.    disp('z=y-2*x/y;');1 _/ \$ q; N8 c, Z% V' @
   
7.    disp('并将该文件保存在work文件夹下');
   ) U- n\" H2 ~: W! j, T
8. end \" F: c, Q; ]: d/ t+ l! ]* v8 N
   
9. 
   : X+ l- M1 G/ Z' R! b\" _& R0 l
10. X1=input('请输入求解区间的左端点X1=');
    & q' l/ z6 f: h\" v
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');; ]& K. g: b  F1 g7 j9 ~
    
12. Xn=input('请输入求解区间的右端点Xn=');
    6 |, v9 p( B. \- o: x3 p4 }
13. h=input('请输入求解步长h=');
    ' I5 c5 L% V, `8 z2 d- k$ t
14. # G- L1 C, V& p2 {
    
15. X=X1;! l1 y  ~0 I- p
    
16. Y=Y1;                                                        %运算初始点
    $ T3 r! X4 o5 _4 o& S
17. n=0;                                                         %节点序号变量置零9 L8 F/ b8 U, \& h- `
    
18. 0 W9 ]( J8 A) T; h
    
19. while X<=Xn-h5 d% F$ ]' ^( N* f/ o7 E& t
    
20.     K1=f(X,Y);3 S  J/ o! [( X
    
21.     K2=f(X+h/2,Y+K1*h/2);/ Z\" u0 x: o# V  c! S4 ~% M
    
22.     K3=f(X+h/2,Y+K2*h/2);
    - ~5 ^* }% P2 Y! S$ D7 s8 {- X
23.     K4=f(X+h,Y+K3*h);
    ( V) F4 \% f# {. c  ^
24.     X=X+h;, v5 Y) V/ x6 U9 N: H. [! k; `
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式( u. K5 N6 O5 m# K4 |$ T
    
26.     n=n+1;                                                   %节点序号加1
    - Z0 j5 M2 m) @, B) D! p' k8 z
27. 
    ' n! M9 P3 r- d
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    ! e! K. [2 N. b& M+ K
29.     plot(X,Y,'o')1 P  l6 I/ k: c7 M
    
30.     hold on/ d) Y* a* B- R$ w. y; `! B* O
    
31. end

复制代码

1. function z=f(x,y), g& p0 @\" u3 k/ O- T
   
2. z=y-2*x/y;

复制代码