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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
$ S1 F4 O! d6 \6 }定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
1 y6 s+ f: r* e5 x4 F0 w* g[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
+ {& S5 M0 `- e  [[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)]
$ Q: M) B  D2 Y其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
# D; {  h6 A7 [: T  f1 r5 L) |) d这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m2 \! X% `3 R) h
   
2. $ ]7 [+ b\" e) @3 M
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)# D8 [) T' ^. b# ]$ _
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');% h1 @/ r4 K( w+ W
   
5.    disp('function z=f(x,y)');
   2 A  Y8 q! Y' \: m1 \, g. b
6.    disp('z=y-2*x/y;');
   5 Q3 Z$ |( P% ?- n4 ]9 ~
7.    disp('并将该文件保存在work文件夹下');4 N$ v. }- Q( P  b; w0 J8 K$ ?
   
8. end 2 Q7 e0 I$ O9 _- k: w0 {0 v; u
   
9. 
   ' |5 h, ~6 A. {9 J& c; Y
10. X1=input('请输入求解区间的左端点X1=');
    ) X# X1 Z2 j! \4 D+ s* C- o' o
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    + b) w: U& F4 _5 J4 C) D
12. Xn=input('请输入求解区间的右端点Xn=');
      i. J& o9 i# a$ r5 o; S
13. h=input('请输入求解步长h=');- h2 x2 I5 j& @4 Q' |+ z
    
14. 9 V5 \1 x7 H1 h
    
15. X=X1;
    : h$ b  N- i; m) n* O( m0 T1 H
16. Y=Y1;                                                        %运算初始点
    * ^  f; w7 G) g\" o6 n  C* S
17. n=0;                                                         %节点序号变量置零) P: V, ^; }9 S  C
    
18. ! X: D\" t( D1 J1 Z/ {3 ?
    
19. while X<=Xn-h
    7 |7 h. S$ i  ?\" |
20.     K1=f(X,Y);
    9 P% y7 t/ C7 Q. k0 o# V. {
21.     K2=f(X+h/2,Y+K1*h/2);! }5 C9 ^5 w; W& ?\" R# n% A
    
22.     K3=f(X+h/2,Y+K2*h/2);
    / d' c3 E/ C8 @0 n  t% o, ?2 }
23.     K4=f(X+h,Y+K3*h);% a+ P! w7 D- n' F6 x
    
24.     X=X+h;
    / \' ^) r5 `' B' K; n
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    0 ?  [5 K6 \: d3 S
26.     n=n+1;                                                   %节点序号加1
    3 H$ J1 G+ `7 X( P- l
27. 1 k% Z% [3 ]& P8 K: R  }# ]0 P
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);) F2 j/ H/ s4 v# J; M
    
29.     plot(X,Y,'o')
    1 Z( a5 N8 u\" t. ~
30.     hold on5 t1 H- Q: y8 q9 b6 T( f
    
31. end

复制代码

1. function z=f(x,y)
   \" W# k& F; @; e\" R! F0 ]3 Y
2. z=y-2*x/y;

复制代码

; n, B% K, ]3 ?( P