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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
2 _) ~: l- n# R, v5 w定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
. Y4 q( g1 |# C% D) z[\frac{dy}{dt} = f(t, y)]
* H3 w* ^# D  Z! G这个方法的迭代公式如下：
/ {6 \+ x$ j, ?+ l% o% K6 ^* ][k1 = h \cdot f(tn, yn)]
' P- _; g8 U% [7 i/ T/ \[k2 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k1}{2})]
[k3 = h \cdot f(tn + \frac{h}{2}, yn + \frac{k2}{2})]
" I* v; w5 e" }$ [[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)) 是微分方程右侧的函数。
6 U# ^  c9 `, l+ {3 `- F# t5 e这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m
   3 a' `  L1 `+ ]$ g+ V+ }! {0 B. U
2. 3 ]+ h4 ?9 j2 t+ D! T. Q
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
     P& q1 U/ c3 D' d9 g% H
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');
   : l3 M1 J. l' c; K\" H# }3 `
5.    disp('function z=f(x,y)');
   # u6 {. D& E1 v
6.    disp('z=y-2*x/y;');+ {% h6 N5 S* T# o# [% k
   
7.    disp('并将该文件保存在work文件夹下');
   ' ^\" h5 \1 G7 y
8. end # d: P$ f7 Y/ Z5 J
   
9. # `/ C0 D4 R, V6 k
   
10. X1=input('请输入求解区间的左端点X1=');
    - f+ ^! d8 ?3 V: n/ f& c
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');% P8 `; y/ T+ d+ u8 y( h\" ^/ G! p
    
12. Xn=input('请输入求解区间的右端点Xn=');& ^$ M9 [& z; m! q( z9 j
    
13. h=input('请输入求解步长h=');
    ( ]5 K7 F\" T\" e$ n$ B) |
14. 
    0 Z* N$ Y5 F' E, U9 c\" \
15. X=X1;6 T4 J6 H: ]1 V: B$ K: h
    
16. Y=Y1;                                                        %运算初始点
    5 n9 C: s: O2 p2 C4 i. Y
17. n=0;                                                         %节点序号变量置零1 x# Z\" X0 R! P: ?! [
    
18. ) I2 w: [- \' o# M( b
    
19. while X<=Xn-h
    0 p- n1 ]/ {. G/ y, E
20.     K1=f(X,Y);\" l2 y$ Q- ?( v& a8 U7 v0 Y+ e( ]
    
21.     K2=f(X+h/2,Y+K1*h/2);8 {8 Y! r' W( r) n
    
22.     K3=f(X+h/2,Y+K2*h/2);
    0 ^! N& }1 ^( x3 Y, @6 Q8 W  W
23.     K4=f(X+h,Y+K3*h);  E! K% X# V% F4 F
    
24.     X=X+h;
    , z) I9 y0 [* E7 O5 Y2 f4 K
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式
    ; d% X3 N) P+ {$ k' z+ g
26.     n=n+1;                                                   %节点序号加1) F* i' j* x7 F2 I. f0 ?1 n
    
27. 
    ) |5 m. ]) I* V# n5 K6 r9 I) E$ e
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    0 W+ ~9 k, p) ~% E0 e& d
29.     plot(X,Y,'o')  |* {- A5 y$ P  f8 A3 F
    
30.     hold on' u% E1 X+ L* B6 N& ~7 \' M( P8 H
    
31. end

复制代码

1. function z=f(x,y)5 L% l4 d4 x' W\" ]\" y; u. t
   
2. z=y-2*x/y;

复制代码

" Q9 ^# p5 w3 P6 _8 c( y