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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
8 {; r3 o; i8 B4 S" y定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
[\frac{dy}{dt} = f(t, y)]
这个方法的迭代公式如下：
[k1 = h \cdot f(tn, yn)]
9 {5 o, k- |- \[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
% K$ U- K  q* v, N& N! u[k4 = h \cdot f(tn + h, yn + k_3)]
[y{n+1} = yn + \frac{1}{6}(k1 + 2k2 + 2k3 + k4)]
: e( f$ y% R) y  y' [: x7 I其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
9 h! e' u2 M1 s8 G这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   7 N3 |/ S. ?. w3 b
2. * X6 A2 T* M9 N# I% V0 M# m\" d3 `
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   + i' \% B1 {$ p6 G( t: P
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');$ f8 Z2 `\" @) e# d5 Z\" L
   
5.    disp('function z=f(x,y)');
   ! v* H: k: o, E\" v* T3 P
6.    disp('z=y-2*x/y;');! W1 S& N, Z4 R: u- G
   
7.    disp('并将该文件保存在work文件夹下');. e2 J# p# N6 d& w# j1 i( e' M
   
8. end 6 N; B) u) g  d+ f( d# _1 R
   
9. 
   * ~7 b5 N! y2 p% h) X; v- }' l
10. X1=input('请输入求解区间的左端点X1=');9 ^1 d* K, ?' r: q: T0 B- V
    
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    * b( W# T: O* {
12. Xn=input('请输入求解区间的右端点Xn=');- w% S) A, m0 h
    
13. h=input('请输入求解步长h=');
    + S; a\" |1 v: I
14. 
    / }) y5 N\" x9 q  p
15. X=X1;
    8 D5 d/ C1 S  @: j# x
16. Y=Y1;                                                        %运算初始点
    . M# J) ?, s) u5 ~) G+ ^
17. n=0;                                                         %节点序号变量置零
    - G& c; M. @2 `; X* B. }: z
18. ; X: c+ }% W- n0 |* e. b' t
    
19. while X<=Xn-h* t; v) [3 \! _% p. e, |  I6 ^
    
20.     K1=f(X,Y);* S+ q. d: g8 O% Z. t: d( T
    
21.     K2=f(X+h/2,Y+K1*h/2);
    7 Z* g* R- |' n: f9 ^/ G' o# e% A, S
22.     K3=f(X+h/2,Y+K2*h/2);
    % e9 P- [% t/ T
23.     K4=f(X+h,Y+K3*h);) D/ J; C6 g# |- H3 [7 V0 q  l( \
    
24.     X=X+h;
    , G8 P3 F/ I& e. G. x+ Y% }4 c% e
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    ) ^; e3 Y% z) w5 L6 F
26.     n=n+1;                                                   %节点序号加1  `' a' N. S- `2 q( V5 `( P* a
    
27. 5 ]& I7 \8 b% p2 h
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);0 i3 M; m2 `% |# H# v
    
29.     plot(X,Y,'o')
    / ^. K; Q' O9 {3 v6 Q
30.     hold on\" q& c* R: u0 y' Z) e+ E
    
31. end

复制代码

1. function z=f(x,y)
   # L2 W9 h\" g% g/ i3 [1 @: q8 ?7 t
2. z=y-2*x/y;

复制代码