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

2744557306

"定步长四阶经典公式"通常指的是数值积分中的四阶Runge-Kutta方法。这是一种常用的数值解常微分方程（ODE）的方法，其主要思想是通过逐步逼近来估计微分方程的解。
# `! G) k  s' t; g" \3 ~定步长四阶经典公式是Runge-Kutta方法的一种，其中最常见的是经典的四阶Runge-Kutta方法。对于一个一阶常微分方程
# _9 ?1 R; l3 _[\frac{dy}{dt} = f(t, y)]
. O4 I$ j0 u: s1 c( Q这个方法的迭代公式如下：
, I/ x( K; S9 q7 g[k1 = h \cdot f(tn, yn)]
5 Z' A& j. f7 ]6 y# y[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)]
5 ]8 a% v/ E0 Q4 Q8 L6 ]其中，(tn) 是当前时间步，(yn) 是当前的解，(h) 是步长，(f(t, y)) 是微分方程右侧的函数。
2 C: _5 s4 ]! t& B' \这个方法的精度相对较高，因为它使用了函数 (f(t, y)) 在一个步长内的多个点上的信息。四阶Runge-Kutta方法在许多情况下被广泛应用，因为它相对简单且相对高效。

1. %四阶经典公式,微分方程为f.m2 }0 x3 u( Y$ w% ^  ?% |- O
   
2. $ U& \\" ]7 Y, s9 f) I0 A
   
3. if exist('f.m')==0                                           %在星号处输入文件名(把星号改为文件名)
   ! J6 ?( r0 c4 v\" U5 L! c6 w- ?8 P
4.    disp('没有为方程创建名为f.m的函数文件,请参照下例建立它');, b# N. Y! o4 i( F
   
5.    disp('function z=f(x,y)');7 J8 n\" S7 J) ~! G2 K! u# ^
   
6.    disp('z=y-2*x/y;');
   5 ?8 V! p\" w( O, o* O. \* A$ s
7.    disp('并将该文件保存在work文件夹下');
   ' C# l0 v5 G5 t& O
8. end 3 p\" r  z6 _+ E9 c( n) ]9 y/ V
   
9. ) I4 u  R+ g4 ~9 i# R
   
10. X1=input('请输入求解区间的左端点X1=');
    ( D! _6 Q\" R' d( ?! S
11. Y1=input('请输入微分方程的初始条件Y1=（X=X1时Y的值）');
    \" A4 D\" b1 e6 D. \
12. Xn=input('请输入求解区间的右端点Xn=');
    3 s/ R: E: u3 J  H2 Q7 T1 Q/ H9 C
13. h=input('请输入求解步长h=');
    / q$ l9 J% f! {. z0 l
14. 
    1 q! U$ h$ L4 W0 l
15. X=X1;  d5 t1 T/ }6 ?' ^( P) h0 t7 E
    
16. Y=Y1;                                                        %运算初始点& l% o+ s1 w3 }# B8 C3 @8 Z6 j& I8 ?
    
17. n=0;                                                         %节点序号变量置零
    0 s/ E0 R4 s/ u) K  d: ^8 C9 c
18. 
    + K1 P  x0 k$ I/ J' [
19. while X<=Xn-h% W+ F7 O4 l& P' H4 g8 [
    
20.     K1=f(X,Y);
    ! ^5 F7 m% R& o! P
21.     K2=f(X+h/2,Y+K1*h/2);
    2 h% Z/ B# ]. P0 {; k- D3 u) K& v
22.     K3=f(X+h/2,Y+K2*h/2);
    & H0 e8 m' K$ w\" j& e0 `
23.     K4=f(X+h,Y+K3*h);1 l2 ^\" \; c$ @9 i
    
24.     X=X+h;1 @5 N* Q3 j+ I, y$ Z! {+ E
    
25.     Y=Y+h*(K1+2*K2+2*K3+K4)/6;                               %四阶标准的龙格-库塔公式+ F5 A- o) y  x* d\" o
    
26.     n=n+1;                                                   %节点序号加16 H; d# ]6 e5 u6 h3 Y
    
27. 2 V, N7 n+ W- X8 c4 L/ F
    
28.     fprintf('第%d个点的计算结果为X=%10.8f,Y=%10.8f\n',n,X,Y);
    : s9 c# g* ~2 Q2 \/ w$ V* H% [
29.     plot(X,Y,'o')2 w0 Q9 l* i& F. c* y% S1 Y5 R
    
30.     hold on$ `9 r  O: _, D8 M
    
31. end

复制代码

1. function z=f(x,y)) N4 K' G* ]* {* t* @! j( q/ }9 g, _% ]$ {
   
2. z=y-2*x/y;

复制代码