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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
: Z3 Z) a. S# v6 A* ^[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
% }. Y' J$ c$ I- {, \! J[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
, I, |8 y( H  E/ Y4 n8 W+ J[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
8 z9 g! g$ O4 [4 f5 Z- |这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   6 l) i2 R5 i+ {
2. 
   , ]8 i0 t: W! @' g\" W& c
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名); M# _4 e( B5 m. Q+ @& }1 v* J  C
   
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');- [\" s3 n3 e' B: n5 x) h
   
5.    disp('function z=f(x,y)');
   / ?+ G1 a. K5 f2 m* [
6.    disp('z=y-2*x/y;');
   : P7 c$ D) D0 D$ ~4 X' Q, f
7.    disp('并将该文件保存在work文件夹下');
   ) {( X$ Q! o$ w% P  v/ g
8. end
   8 r6 r5 K; t\" U  P
9. 
   , [3 W  {6 P# v1 w' n, G+ e
10. X1=input('请输入求解区间的左端点X1=');
    5 t: m3 I3 C\" z
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');; z9 A0 r: ?+ c( |5 {
    
12. Xn=input('请输入求解区间的右端点Xn=');
    . T: W1 G% v1 {
13. h=input('请输入求解步长h=');
    0 b, w7 ^\" g9 F4 w* c% _5 S- G
14. 
    ( M% Y. h+ K4 _! ^  B& {' T\" O  M
15. X=X1;7 X( o$ S- ^\" e9 E1 {/ ?& h
    
16. Y=Y1;                                                        %运算初始点
    9 m  |8 l! p  `# u
17. n=0;                                                         %节点序号变量置零3 O; }$ s2 [& X3 \. O
    
18. 
    0 ^  X& G! J# G4 F
19. while X<=Xn-h3 r9 b6 B+ @\" V9 J6 L  p; \
    
20.     K1=f(X,Y);) Z! Z, z! x) Y! W5 K3 i0 R
    
21.     K2=f(X+h/2,Y+K1*h/2);. [: r- D. C\" O0 u9 w/ ~8 F
    
22.     K3=f(X+h/2,Y+K2*h/2);
    % Z- p3 W7 c. [6 H
23.     K4=f(X+h,Y+K3*h);* h4 O/ r4 o, ?0 @
    
24.     X=X+h;
    ; L4 f3 Y( U3 T- t
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式\" \& g( q& D- m$ c
    
26.     n=n+1;                                                   %节点序号加1
    0 ^/ }, a6 O5 [$ O4 }( U/ I) [
27. 
    1 h9 D# H3 J, r$ r8 `  M5 {, C0 v
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    7 J6 `( B: n5 L9 |7 J
29.     plot(X,Y,'o')
    4 O8 i5 u; |, @, S# _  _
30.     hold on/ t5 Q, K1 z4 V  ^8 ?. x; v
    
31. end

复制代码

1. function z=f(x,y)
   ! n7 C2 s+ M& A3 m9 q
2. z=y-2*x/y;

复制代码

* Q5 N5 x+ Y; n( j8 r6 C( l! |% B