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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
( L5 a( P+ @4 B2 r1 S5 Y7 P: q[\frac{dy}{dt} = f(t, y)]
% F  ^6 N  z; E; Y这个方法的迭代公式如下：
: Z0 ^- {3 H# N- z2 k[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})]
[k4 = h \cdot f(tn + h, yn + k_3)]
) y. v! Y$ x- P3 `$ u9 {5 {[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
     q; K# g* o9 d; d) i\" D
2. 
   & p1 |7 @( ^, d7 }\" e6 Y$ u
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   * B  k. e) A, W& ^8 k( ~
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');$ K\" B' y& ^/ c9 ~5 b* k
   
5.    disp('function z=f(x,y)');' f6 T0 U( s3 s, W) g# o
   
6.    disp('z=y-2*x/y;');0 u# I4 R' |8 ~) i! x+ T' O( v
   
7.    disp('并将该文件保存在work文件夹下');2 H: @$ F: V9 g% i
   
8. end 5 c: E/ V7 t  j! i! u3 q/ y
   
9. 
   3 t  E8 J+ T3 @8 G
10. X1=input('请输入求解区间的左端点X1=');
    5 \  c+ j% t; z& ~& C; K' [* Q
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');) L( H( i( z% L
    
12. Xn=input('请输入求解区间的右端点Xn=');\" F9 I% x, a' a# _
    
13. h=input('请输入求解步长h=');8 `6 o; B# Y4 q, q/ ^3 N
    
14. 0 I- {8 d* q! J. _( G2 L# ^6 H
    
15. X=X1;
    / F3 {; l. b\" L2 F4 @( ^$ V
16. Y=Y1;                                                        %运算初始点
    . m' I7 b0 B( n6 b
17. n=0;                                                         %节点序号变量置零
    7 x  l3 r: n; `\" }7 x& ~
18. 
    * t. B: I  I0 j6 O\" ~3 }# L
19. while X<=Xn-h
    2 j/ R\" r1 V% O6 Y  H. C( I0 V
20.     K1=f(X,Y);
    . u7 ]# N& L; F3 Y
21.     K2=f(X+h/2,Y+K1*h/2);' F8 ~: n9 K* t8 P4 U; {5 d+ D
    
22.     K3=f(X+h/2,Y+K2*h/2);# }( ~$ P1 r; O% c0 u6 y
    
23.     K4=f(X+h,Y+K3*h);2 q/ R2 u! \: i, D+ Z  O
    
24.     X=X+h;' F% _% y! R; U+ K+ H8 s% T
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    2 I) R3 Q0 L\" T# C0 C
26.     n=n+1;                                                   %节点序号加15 `) w7 Y, L. N  u
    
27. 
    3 A+ W. R( M/ o; o/ X
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    \" M/ x, f% o# @5 w
29.     plot(X,Y,'o')3 N1 ^$ g0 o7 R4 Y. t$ ~
    
30.     hold on
    * R! _! J0 T) `0 j
31. end

复制代码

1. function z=f(x,y)* B, D\" w. z/ ~/ k: y! T
   
2. z=y-2*x/y;

复制代码