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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
( y8 t! A+ C: }5 ]1 y, t' J定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
4 L  n0 a; H. J! h. d[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
1 k0 V& A* R% v% [/ t- F  x[k1 = h \cdot f(tn, yn)]
! z' |, [& w% Y5 _! p. m- m[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
  V% j$ J3 S# j$ H( D[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
[k4 = h \cdot f(tn + h, yn + k_3)]
# z% {: S; v: t4 B& H( k; A# H[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
" c. B+ S' y0 `( H其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   : H. ?4 B! p  u
2. 
   ) p- @: p( l. R# Q; }
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   0 w, f4 T\" r) B! g7 r6 U
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
     y, X. I( K- E3 ~# `
5.    disp('function z=f(x,y)');
   / `, s\" M0 }7 [0 |8 E3 D3 r  P6 J- v
6.    disp('z=y-2*x/y;');
   / i6 s5 m' ~0 R( e3 A
7.    disp('并将该文件保存在work文件夹下');9 z/ v% Q/ l8 t# `
   
8. end ! ?. ^, `0 c8 t7 U/ T6 e+ Z& u
   
9. 
   $ f7 O' |2 n+ E/ a: @
10. X1=input('请输入求解区间的左端点X1=');
      s! v6 U$ D6 X; f4 ]$ a
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    & P\" ^9 q( t2 g5 Q
12. Xn=input('请输入求解区间的右端点Xn=');9 e& i; C1 k8 E* [1 ^! ]: R; w; V' c
    
13. h=input('请输入求解步长h=');! u  o. v& P* P2 u5 {9 N
    
14. * \% b! N/ A; {( O2 N
    
15. X=X1;, y\" a( p9 n) C) z
    
16. Y=Y1;                                                        %运算初始点
    - ]5 f: w: F' @' p5 D# J2 F  B6 Q
17. n=0;                                                         %节点序号变量置零6 n8 ?\" U4 J\" z: e& U1 a7 x
    
18. 
    9 w1 n  O4 Q0 [7 x
19. while X<=Xn-h, X! J* A' x* {: O+ h4 P5 `4 t
    
20.     K1=f(X,Y);
    ( W, W\" y9 h( V
21.     K2=f(X+h/2,Y+K1*h/2);+ q1 ]( n1 |% x' e7 M
    
22.     K3=f(X+h/2,Y+K2*h/2);
    5 d. q# R# Y% A; I9 q
23.     K4=f(X+h,Y+K3*h);/ I) y. R& t. _, L2 b0 f
    
24.     X=X+h;
    2 V) s$ `% U- r8 R9 E3 v
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    8 U+ l/ o: L7 h2 e
26.     n=n+1;                                                   %节点序号加13 ]4 l; X7 i4 i0 \1 t3 C# d7 c
    
27. ) A1 ^( {: [. ^
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    ' [$ K9 B2 X4 `- a8 E& x0 S
29.     plot(X,Y,'o')# r6 O( E& H8 n, B7 J3 l+ ]
    
30.     hold on& w7 y9 J- x$ ^* B4 V$ c7 J
    
31. end

复制代码

1. function z=f(x,y); `0 ]# a) ^5 h1 ~
   
2. z=y-2*x/y;

复制代码