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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
5 u! Y4 j( r7 j+ V[\frac{dy}{dt} = f(t, y)]
+ Z; r, |# d' v2 {: S  m: B# Z这个方法的迭代公式如下：
. U- \8 r/ p* T[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})]
% E, ?4 z$ t7 g1 C  g! ~2 d# }[k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
! s- u) O- p# _; Y. `: r5 P4 a0 F其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
) {4 Y+ X* K: X2 g这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   ; U& Q* ]* _$ T2 g( D* E. b9 D
2. 
   / u8 J\" R% ]5 \* Y
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   3 \9 n+ E# z* c; j6 O: E
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');7 `: C- D- K- B% H% v
   
5.    disp('function z=f(x,y)');
   9 |1 P3 J7 Z2 B. G
6.    disp('z=y-2*x/y;');6 Z) Q! p\" C; M
   
7.    disp('并将该文件保存在work文件夹下');
   ! H3 _9 B1 i, J+ ?$ P
8. end 7 T6 q* S8 v( ]8 Z
   
9. 
   0 F  @0 X! s- C8 o' {
10. X1=input('请输入求解区间的左端点X1=');- N( l! ?; B3 W9 Y- q/ ^, R( Y
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');) F/ v2 V9 u. T3 C  v  ~0 E
    
12. Xn=input('请输入求解区间的右端点Xn=');; K1 C) b7 d. F8 A
    
13. h=input('请输入求解步长h=');
    * Q\" I$ X- s0 W- H% u. m% @
14. 
    & N3 Y- t- E. _. z7 I8 E\" L+ e1 N
15. X=X1;
    . m7 l4 [* l7 c) G
16. Y=Y1;                                                        %运算初始点
    : D7 }' U' z( ?0 v1 @
17. n=0;                                                         %节点序号变量置零
    - _/ U* {8 A  i6 U3 c* m! F
18. 0 l  ]6 c: q& ~$ L  e; N( d
    
19. while X<=Xn-h
    & f% E5 I\" j1 i, J5 X% t
20.     K1=f(X,Y);4 P0 M( w- Y4 v- U6 l
    
21.     K2=f(X+h/2,Y+K1*h/2);
    ) H0 m  h, q8 G
22.     K3=f(X+h/2,Y+K2*h/2);
    \" t) t4 d# T\" _- w* J
23.     K4=f(X+h,Y+K3*h);
    & `! X# L- y1 g& W7 B\" g
24.     X=X+h;5 ^& X0 g# q! S
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式4 |: M6 \9 f/ N! L3 ~+ \* r4 J
    
26.     n=n+1;                                                   %节点序号加1
    9 k. E! i4 M+ J$ I* w0 |: I
27. 
    6 e$ Q; B% M% |8 L
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    ; \\" b4 z8 z% v5 a2 Z7 [# B: Y
29.     plot(X,Y,'o')  ~. u0 r0 A8 {8 I4 j
    
30.     hold on2 ]7 O2 k\" |/ r+ v  M* S1 i, K8 ^! }
    
31. end

复制代码

1. function z=f(x,y). u7 ^9 r/ Z& X4 Z' H$ p: d, ~3 ]
   
2. z=y-2*x/y;

复制代码