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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
4 b9 ]0 m* G# R9 @, e# B3 j$ X4 ]6 K[\frac{dy}{dt} = f(t, y)]
6 I8 I: v8 K- M1 H这个方法的迭代公式如下：
$ K/ B: s$ o+ z[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{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m$ R6 [( p9 J/ _- r
   
2. 
   5 I9 {) j9 R4 T
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   + S% n/ ~! F* v9 F4 @
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');4 z4 j+ `5 f7 A7 W, ]# [
   
5.    disp('function z=f(x,y)');
   : l3 \% X5 K2 j\" |
6.    disp('z=y-2*x/y;');
   : c2 q' j1 k* L- Y
7.    disp('并将该文件保存在work文件夹下');
   * N2 ]$ n% {* O
8. end
   * T4 L, G3 `2 ~# E- t! b4 z
9. 
   0 \, U' U# D$ {( p' o
10. X1=input('请输入求解区间的左端点X1=');
    4 I9 ^6 b, k5 B, t8 m, w1 k: ?2 f, i
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    & q& M: Q, p! V
12. Xn=input('请输入求解区间的右端点Xn=');4 M0 s# {% L' M, o8 _0 A2 S* l
    
13. h=input('请输入求解步长h=');
    \" B& `8 a8 K' O
14. 
    5 X1 O9 J6 x1 L9 f  a9 V
15. X=X1;7 \! G3 f+ s( P) ~
    
16. Y=Y1;                                                        %运算初始点
    : h6 O* r8 F! \/ G0 N
17. n=0;                                                         %节点序号变量置零/ s/ }, [/ w( ^\" a' n1 Q& R7 o\" g
    
18. 
      Q, e1 r) H( x, z! |( q$ C: U
19. while X<=Xn-h. w5 B4 ]# k7 _( a0 s8 a- @7 K. S
    
20.     K1=f(X,Y);# P8 |8 G6 w3 B+ P1 F$ R
    
21.     K2=f(X+h/2,Y+K1*h/2);
    ! V3 |/ P6 U\" u- f5 |
22.     K3=f(X+h/2,Y+K2*h/2);! U- ~# H) K9 P0 i7 b3 r1 `5 L) v
    
23.     K4=f(X+h,Y+K3*h);5 @1 B\" G  Z5 t( m
    
24.     X=X+h;: P+ i, H& k( H# m
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    0 |4 G8 D! U+ C  Q0 s) Q$ p
26.     n=n+1;                                                   %节点序号加1
    9 F1 L( ~8 ~2 Q& I  l
27. 
    * Z, A; u% L2 F* \7 ]0 b
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    7 v- Q( z# P  J  q# t7 |- C
29.     plot(X,Y,'o')5 K- }1 K* }9 Q
    
30.     hold on
    ; R( ^4 U) F* N. v3 n
31. end

复制代码

1. function z=f(x,y)
   / j& a% q8 d5 I$ v* }
2. z=y-2*x/y;

复制代码