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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
% K" {- J' v% R; s4 u' h" b[\frac{dy}{dt} = f(t, y)]
. B: u7 j, K, z- c这个方法的迭代公式如下：
: y; b1 z8 \( d. C[k1 = h \cdot f(tn, yn)]
[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
  g1 p  j7 R6 z4 A/ N3 B[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
$ q& D9 K; V3 X. l[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
" k! c# L" f- x" P2 y其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
! z7 u) T& t+ |; j& Z这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m; K( ~3 ^+ @3 N0 {8 ^
   
2. ' w; B: r9 [* F( ^% F( b7 F
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)2 Q% `2 F/ Q2 p' D5 W3 K) Y5 w
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');( `0 U7 `6 K  R4 O. W
   
5.    disp('function z=f(x,y)');
   & i6 V9 O3 h& p9 T
6.    disp('z=y-2*x/y;');
   4 W1 N# h) s& [! q. B5 o; m
7.    disp('并将该文件保存在work文件夹下');
   6 j, B3 r3 B% T# i: I: {* W5 ^3 h
8. end $ [, d9 g( a2 I8 K
   
9. ) K+ c9 U; X3 ?/ Z
   
10. X1=input('请输入求解区间的左端点X1=');
    . o& ~6 `4 ]/ X! Z6 C
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');- ?5 q7 n3 G& q- C; e4 I& X
    
12. Xn=input('请输入求解区间的右端点Xn=');
    * `/ P3 \\" C9 S8 k  n
13. h=input('请输入求解步长h=');
    / u' d' \: _1 q8 l- ^, y0 _7 T2 Y
14. / {- |& E\" w6 W5 e* f+ T
    
15. X=X1;! Z( f, Q* K8 f# O. t9 x
    
16. Y=Y1;                                                        %运算初始点0 t  D* Y) W# |. `% @  U$ B
    
17. n=0;                                                         %节点序号变量置零% x. O% x: ?# [2 \4 w3 m
    
18. 
    ! r! Y) |( o8 E6 b, a
19. while X<=Xn-h
    9 k( k; ]% A9 F* k9 r9 V5 w9 o' z
20.     K1=f(X,Y);  e  C0 w% y( g6 J
    
21.     K2=f(X+h/2,Y+K1*h/2);) {& s8 P( n* w
    
22.     K3=f(X+h/2,Y+K2*h/2);- f$ f' L  g2 l& g2 S
    
23.     K4=f(X+h,Y+K3*h);
    % q# c: ~! Z- B$ K  \% p
24.     X=X+h;! s0 }: ^& ^9 y( l) w* J% Z
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式0 l0 {' _( y1 O* m
    
26.     n=n+1;                                                   %节点序号加1# u0 E+ }% l% G% @6 N6 N
    
27. 
    , n) e8 D  ~- n9 V
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);\" E3 t/ J' [; }& e
    
29.     plot(X,Y,'o')
    % M' I5 b4 ]! f4 X* j7 p; a
30.     hold on0 f& ~& R2 b8 `# N3 r
    
31. end

复制代码

1. function z=f(x,y)
   5 @6 g! z1 u+ @3 j
2. z=y-2*x/y;

复制代码

& ~# f3 R' H1 L) U0 y+ k