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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
" O/ T6 a6 I! k" ^7 ]& w8 y5 s. n[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
7 Q5 e. h5 D( K& g/ F  D[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
8 v% H7 j) ^' i0 R, `: ^7 Z[k4 = h \cdot f(tn + h, yn + k_3)]
9 n# v: ^+ S# R# g" O[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
( y, s2 E  L: a6 i其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
; ~* i# z7 j4 b, E. _& z这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m, Q5 d0 y& h1 Y( D0 ~0 k$ ?7 q
   
2. : s\" j- R( y$ I/ Q$ u
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   5 C3 b+ K* d' @
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');/ V9 v6 `* Y( i& V
   
5.    disp('function z=f(x,y)');
   : P! p, `. a; g+ E, w
6.    disp('z=y-2*x/y;');; J% D# @# U; n  ]5 m
   
7.    disp('并将该文件保存在work文件夹下');+ f! o  `9 N  M' Q  k, g: z, v
   
8. end 9 e5 T& t1 H% T
   
9. 1 W4 G3 r7 P3 U! O
   
10. X1=input('请输入求解区间的左端点X1=');( U& d% h5 k0 @! G2 f, E
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    ' J& u7 s: s# l1 {& Z/ X' N3 A
12. Xn=input('请输入求解区间的右端点Xn=');
    : A# d7 k) X% a\" F' g
13. h=input('请输入求解步长h=');
    - M8 e& }# X- s. c/ J! f1 {
14. 
    ! A\" [- `6 y  A
15. X=X1;
    ' a6 }( y' W1 q5 _' @2 x
16. Y=Y1;                                                        %运算初始点
    7 T5 j# Z, P4 e
17. n=0;                                                         %节点序号变量置零
    ) q/ d* ^/ {3 r+ _
18. 8 N$ C; X  S# v/ e. J
    
19. while X<=Xn-h, Q2 ~1 ^; q7 e$ k
    
20.     K1=f(X,Y);
    ( N% y2 V+ O. x* v* O0 m$ F
21.     K2=f(X+h/2,Y+K1*h/2);$ w8 f* o! t) F0 d
    
22.     K3=f(X+h/2,Y+K2*h/2);
    8 [$ R' q& d6 |* b6 S
23.     K4=f(X+h,Y+K3*h);- u& B  g5 P6 n. [
    
24.     X=X+h;- a. }$ O# r) R5 B, M. t
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式4 n  N5 S2 q2 H* R/ `5 i1 Y
    
26.     n=n+1;                                                   %节点序号加1
    / |) A! w: N8 ?9 W0 D\" B5 i6 e
27. 
    ( b0 U1 B  y2 C1 ~2 z
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);( p9 g- U' n7 h7 u
    
29.     plot(X,Y,'o')2 r5 _, _1 F) b4 ~
    
30.     hold on
    % X% C4 B- W8 U) `& m2 R+ }: B
31. end

复制代码

1. function z=f(x,y)
   & {2 ?  R% N) P3 _# d; n7 d
2. z=y-2*x/y;

复制代码