文章目录
1. 概念
1) 常微分方程
自变量只有一个的微分方程,称为常微分方程;自变量数量2个或以上时,称为偏微分方程。
绝大多数实际工程问题,常微分方程的自变量都是时间t,通常表达为:
F ( t , y , y ˙ , ⋯ , y ( n − 1 ) , y ( n ) ) = 0 F\left( {t,y,\dot y, \cdots ,{y^{\left( {n - 1} \right)}},{y^{\left( n \right)}}} \right) = 0 F(t,y,y˙,⋯,y(n−1),y(n))=0
2) 一阶常微分方程
求解上述关于 y ( n ) y^{\left( n \right)} y(n)的方程,得到 y ( n ) y^{\left( n \right)} y(n)的表达式为:
y ( n ) = f ( t , y , y ˙ , ⋯ , y ( n − 1 ) ) {y^{\left( n \right)}} = f\left( {t,y,\dot y, \cdots ,{y^{\left( {n - 1} \right)}}} \right) y(n)=f(t,y,y˙,⋯,y(n−1))
写成向量形式:
[ y ˙ y ¨ ⋮ y ( n ) ] = [ y ¨ y ( 3 ) ⋮ f ( t , y , y ˙ , ⋯ , y ( n − 1 ) ) ] \left[ \begin{matrix} {\dot{y}} \\ {\ddot{y}} \\ \vdots \\ {{y}^{\left( n \right)}} \\ \end{matrix} \right]=\left[ \begin{matrix} {\ddot{y}} \\ {{y}^{\left( 3 \right)}} \\ \vdots \\ f\left( t,y,\dot{y},\cdots ,{{y}^{\left( n-1 \right)}} \right) \\ \end{matrix} \right] ⎣⎢⎢⎢⎡y˙y¨⋮y(n)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y¨y(3)⋮f(t,y,y˙,⋯,y(n−1))⎦⎥⎥⎥⎤
令:
y = [ y y ˙ ⋮ y ( n − 1 ) ] \mathbf{y}=\left[ \begin{matrix} y \\ {\dot{y}} \\ \vdots \\ {{y}^{\left( n-1 \right)}} \\ \end{matrix} \right] y=⎣⎢⎢⎢⎡yy˙⋮y(n−1)⎦⎥⎥⎥⎤
f ( t , y ) = [ y ¨ y ( 3 ) ⋮ f ( t , y , y ˙ , ⋯ , y ( n − 1 ) ) ] = [ y ¨ y ( 3 ) ⋮ f ( t , y ) ] f\left( t,\mathbf{y} \right)=\left[ \begin{matrix} {\ddot{y}} \\ {{y}^{\left( 3 \right)}} \\ \vdots \\ f\left( t,y,\dot{y},\cdots ,{{y}^{\left( n-1 \right)}} \right) \\ \end{matrix} \right]=\left[ \begin{matrix} {\ddot{y}} \\ {{y}^{\left( 3 \right)}} \\ \vdots \\ f\left( t,\mathbf{y} \right) \\ \end{matrix} \right] f(t,y)=⎣⎢⎢⎢⎡y¨y(3)⋮f(t,y,y˙,⋯,y(n−1))⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y¨y(3)⋮f(t,y)⎦⎥⎥⎥⎤
得到向量形式的一阶常微分方程:
y
˙
=
f
(
t
,
y
)
\mathbf{\dot{y}}=f\left( t,\mathbf{y} \right)
y˙=f(t,y)
3) 微分方程的解析解法和数值解法
微分方程的求解方法有解析解法和数值解法,解析法是求出因变量
y
{\bf{y}}
y
关于时间
t
t
t的具体函数式,表达为
y
=
y
(
t
)
{\bf{y}} = {\bf{y}}\left( t \right)
y=y(t);数值法是解出因变量
y
{\bf{y}}
y关于时间
t
t
t的离散序列,通常表达为
{
t
(
k
)
,
y
(
k
)
}
\left\{ {t(k),{\bf{y}}(k)} \right\}
{t(k),y(k)}。
绝大多数的非线性常微分方程,不存在或难以求出解析解,大多数情况下只能求取微分方程的数值解。
2. 算法
欧拉法(Euler)是求解一阶常微分方程初值问题的数值方法,分为显式欧拉法、隐式欧拉法、两步欧拉法和改进欧拉法。
2.1 显式欧拉法
一阶微分方程初始问题:
{ y ˙ = f ( t , y ) y ( t 0 ) = y 0 \left\{ \begin{array}{l} {\bf{\dot y}} = f\left( {t,{\bf{y}}} \right) \\ {\bf{y}}\left( {{t_0}} \right) = {{\bf{y}}_0} \\ \end{array} \right. {y˙=f(t,y)y(t0)=y0
式中,
t
0
{t_0}
t0为初始时间(已知常数),
y
0
{{\bf{y}}_0}
y0为初始状态(已知向量),
f
(
t
,
y
)
f\left( {t,{\bf{y}}} \right)
f(t,y)为关于时间
t
{t}
t和状态
y
{{\bf{y}}}
y的函数(已知函数)。
使用一阶向前差商代替微分,即:
y ˙ ≈ y ( k + 1 ) − y ( k ) t ( k + 1 ) − t ( k ) = y ( k + 1 ) − y ( k ) h {\bf{\dot y}} \approx \frac{{{\bf{y}}\left( {k + 1} \right) - {\bf{y}}\left( k \right)}}{{t\left( {k + 1} \right) - t\left( k \right)}} = \frac{{{\bf{y}}\left( {k + 1} \right) - {\bf{y}}\left( k \right)}}{h} y˙≈t(k+1)−t(k)y(k+1)−y(k)=hy(k+1)−y(k)
式中,h为时间步长。
微分方程变为显式差分方程:
{ y ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) y ( 0 ) = y 0 \left\{ \begin{array}{l} {\bf{y}}\left( {k + 1} \right) = {\bf{y}}\left( k \right) + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) \\ {\bf{y}}\left( 0 \right) = {{\bf{y}}_0} \\ \end{array} \right. {y(k+1)=y(k)+hf(t(k),y(k))y(0)=y0
上式是关于 y ( k − 1 ) {\bf{y}}\left( k-1 \right) y(k−1)向 y ( k ) {\bf{y}}\left( k \right) y(k)的递推形式,可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {\bf{y}}\left( 1 \right),{\bf{y}}\left( 2 \right),{\bf{y}}\left( 3 \right),{\bf{y}}\left( 4 \right) \cdots ,{\bf{y}}\left( N \right) \cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯,此离散序列即为微分方程的数值解。
2.2 隐式欧拉法
使用一阶向后差商代替微分,即
y ˙ ≈ y ( k ) − y ( k − 1 ) t ( k ) − t ( k − 1 ) = y ( k ) − y ( k − 1 ) h {\bf{\dot y}} \approx \frac{{{\bf{y}}\left( k \right) - {\bf{y}}\left( {k - 1} \right)}}{{t\left( k \right) - t\left( {k - 1} \right)}} = \frac{{{\bf{y}}\left( k \right) - {\bf{y}}\left( {k - 1} \right)}}{h} y˙≈t(k)−t(k−1)y(k)−y(k−1)=hy(k)−y(k−1)
微分方程变为隐式差分方程:
{ y ( k ) = y ( k − 1 ) + h f ( t ( k ) , y ( k ) ) y ( 0 ) = y 0 \left\{ \begin{array}{l} {\bf{y}}\left( k \right) = {\bf{y}}\left( {k - 1} \right) + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) \\ {\bf{y}}\left( 0 \right) = {{\bf{y}}_0} \\ \end{array} \right. {y(k)=y(k−1)+hf(t(k),y(k))y(0)=y0
在第 k k k步迭代时, y ( k − 1 ) , t ( k ) {\bf{y}}\left( {k - 1} \right),t(k) y(k−1),t(k)已知, y ( k ) {\bf{y}}\left( k \right) y(k)为待求未知量, f ( t ( k ) , y ( k ) ) f\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) f(t(k),y(k))是关于待求未知量 y ( k ) {\bf{y}}\left( k \right) y(k)的函数,即方程为关于 y ( k ) {\bf{y}}\left( k \right) y(k)的非线性方程。通过非线性方程的迭代解法(如牛顿迭代法)可求解出 y ( k ) {\bf{y}}\left( k \right) y(k)。依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {\bf{y}}\left( 1 \right),{\bf{y}}\left( 2 \right),{\bf{y}}\left( 3 \right),{\bf{y}}\left( 4 \right) \cdots ,{\bf{y}}\left( N \right) \cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯,此离散序列即为微分方程的数值解。
2.3 两步欧拉法
使用中心差商代替微分:
y ˙ ≈ y ( k + 1 ) − y ( k − 1 ) t ( k + 1 ) − t ( k − 1 ) = y ( k + 1 ) − y ( k − 1 ) 2 h {\bf{\dot y}} \approx \frac{{{\bf{y}}\left( {k + 1} \right) - {\bf{y}}\left( {k - 1} \right)}}{{t\left( {k + 1} \right) - t\left( {k - 1} \right)}} = \frac{{{\bf{y}}\left( {k + 1} \right) - {\bf{y}}\left( {k - 1} \right)}}{{2h}} y˙≈t(k+1)−t(k−1)y(k+1)−y(k−1)=2hy(k+1)−y(k−1)
微分方程变为显式差分方程:
{ y ( k + 1 ) = y ( k − 1 ) + 2 h f ( t ( k ) , y ( k ) ) y ( 0 ) = y 0 \left\{ \begin{array}{l} {\bf{y}}\left( {k + 1} \right) = {\bf{y}}\left( {k - 1} \right) + 2hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) \\ {\bf{y}}\left( 0 \right) = {{\bf{y}}_0} \\ \end{array} \right. {y(k+1)=y(k−1)+2hf(t(k),y(k))y(0)=y0
上式是关于 y ( k − 1 ) {\bf{y}}\left( k-1 \right) y(k−1), y ( k ) {\bf{y}}\left( k \right) y(k)向 y ( k + 1 ) {\bf{y}}\left( k+1 \right) y(k+1)的递推形式,可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {\bf{y}}\left( 1 \right),{\bf{y}}\left( 2 \right),{\bf{y}}\left( 3 \right),{\bf{y}}\left( 4 \right) \cdots ,{\bf{y}}\left( N \right) \cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯,此离散序列即为微分方程的数值解。
这里需要特别注意 y ( 1 ) {\bf{y}}\left( 1 \right) y(1)的计算方法:按照公式计算 y ( 1 ) {\bf{y}}\left( 1 \right) y(1)需要用到无意义的 y ( − 1 ) {\bf{y}}\left( -1 \right) y(−1),此步需要转换成显式欧拉法计算:
{ y ( 1 ) = y ( 0 ) + h f ( t ( 0 ) , y ( 0 ) ) y ( 0 ) = y 0 \left\{ \begin{array}{l} {\bf{y}}\left( 1 \right) = {\bf{y}}\left( 0 \right) + hf\left( {t\left( 0 \right),{\bf{y}}\left( 0 \right)} \right) \\ {\bf{y}}\left( 0 \right) = {{\bf{y}}_0} \\ \end{array} \right. {y(1)=y(0)+hf(t(0),y(0))y(0)=y0
进一步分析:
y ( k + 1 ) = y ( k − 1 ) + 2 h f ( t ( k ) , y ( k ) ) = [ y ( k − 1 ) + h f ( t ( k ) , y ( k ) ) ] + h f ( t ( k ) , y ( k ) ) {\bf{y}}\left( {k + 1} \right) = {\bf{y}}\left( {k - 1} \right) + 2hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) = \left[ {{\bf{y}}\left( {k - 1} \right) + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right)} \right] + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) y(k+1)=y(k−1)+2hf(t(k),y(k))=[y(k−1)+hf(t(k),y(k))]+hf(t(k),y(k))
式中, y ˉ ( k ) = y ( k − 1 ) + h f ( t ( k ) , y ( k ) ) {\bf{\bar y}}\left( k \right) = {\bf{y}}\left( {k - 1} \right) + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) yˉ(k)=y(k−1)+hf(t(k),y(k)) 相当于隐式欧拉法由 y ( k − 1 ) {\bf{y}}\left( k-1 \right) y(k−1)求解 y ( k ) {\bf{y}}\left( k \right) y(k); y ( k + 1 ) = y ˉ ( k ) + f ( t ( k ) , y ( k ) ) {\bf{y}}\left( {k + 1} \right) = {\bf{\bar y}}\left( k \right) + f\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) y(k+1)=yˉ(k)+f(t(k),y(k))相当于显式欧拉法由 y ( k ) {\bf{y}}\left( k \right) y(k)求解 y ( k + 1 ) {\bf{y}}\left( k+1 \right) y(k+1)。两步欧拉法是隐式欧拉法与显式欧拉法相结合的一种算法。
2.4 改进欧拉法
先使用显式欧拉法求解 k + 1 k+1 k+1步的预测值:
y ˉ ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) {\bf{\bar y}}\left( {k + 1} \right) = {\bf{y}}\left( k \right) + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) yˉ(k+1)=y(k)+hf(t(k),y(k))
再使用 k k k步的微分值 f ( t ( k ) , y ( k ) ) f\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) f(t(k),y(k))和 k + 1 k+1 k+1步预测的微分值 f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) f\left( {t\left( {k + 1} \right),{\bf{\bar y}}\left( {k + 1} \right)} \right) f(t(k+1),yˉ(k+1))的平均值来近似微分:
y ˙ ≈ f ( t ( k ) , y ( k ) ) + f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) 2 {\bf{\dot y}} \approx \frac{{f\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) + f\left( {t\left( {k + 1} \right),{\bf{\bar y}}\left( {k + 1} \right)} \right)}}{2} y˙≈2f(t(k),y(k))+f(t(k+1),yˉ(k+1))
再利用显式欧拉法校正 k + 1 k+1 k+1步的计算值:
y ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) + f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) 2 {\bf{y}}\left( {k + 1} \right) = {\bf{y}}\left( k \right) + h\frac{{f\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) + f\left( {t\left( {k + 1} \right),{\bf{\bar y}}\left( {k + 1} \right)} \right)}}{2} y(k+1)=y(k)+h2f(t(k),y(k))+f(t(k+1),yˉ(k+1))
综上所述,改进欧拉法的算法为:
{ y ˉ ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) y ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) + f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) 2 y ( 0 ) = y 0 \left\{ \begin{array}{l} {\bf{\bar y}}\left( {k + 1} \right) = {\bf{y}}\left( k \right) + hf\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) \\ {\bf{y}}\left( {k + 1} \right) = {\bf{y}}\left( k \right) + h\frac{{f\left( {t\left( k \right),{\bf{y}}\left( k \right)} \right) + f\left( {t\left( {k + 1} \right),{\bf{\bar y}}\left( {k + 1} \right)} \right)}}{2} \\ {\bf{y}}\left( 0 \right) = {{\bf{y}}_0} \\ \end{array} \right. ⎩⎨⎧yˉ(k+1)=y(k)+hf(t(k),y(k))y(k+1)=y(k)+h2f(t(k),y(k))+f(t(k+1),yˉ(k+1))y(0)=y0
上式是关于 y ( k ) {\bf{y}}\left( k \right) y(k)向 y ( k + 1 ) {\bf{y}}\left( k+1 \right) y(k+1)的递推形式,可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {\bf{y}}\left( 1 \right),{\bf{y}}\left( 2 \right),{\bf{y}}\left( 3 \right),{\bf{y}}\left( 4 \right) \cdots ,{\bf{y}}\left( N \right) \cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯,此离散序列即为微分方程的数值解。
2.5 四种欧拉方法的对比
微分方程的数值解法,本质是使用数值积分来实现
y
˙
{\bf{\dot y}}
y˙向
y
{\bf{y}}
y的转换。显式欧拉法和隐式欧拉法,使用的是矩形数值积分方法,矩形积分是一阶算法,其截断误差为
O
(
h
2
)
O\left( {{h^2}} \right)
O(h2);两步欧拉法和改进欧拉法,使用的是梯形数值积分方法,梯形积分是二阶算法,其截断误差为
O
(
h
3
)
O\left( {{h^3}} \right)
O(h3)。
算法类型 | 计算量 | 计算精度 | 稳定性 |
---|---|---|---|
显式欧拉法 | 小 | 低 | 差 |
隐式欧拉法 | 大 | 低 | 差 |
两步欧拉法 | 中 | 较高 | 中 |
改进欧拉法 | 中 | 高 | 好 |
3. 程序
作者使用Matlab开发了四种欧拉方法求解常微分方程的程序,能够方便快捷的求解一阶常微分方程的初值问题。
function [T,X,dX] = ODE_ExplicitEuler( Hfun,t,h,x0 )
% [T,X,dX] = ODE_ExplicitEuler( Hfun,t,h,x0 ) 显式欧拉法求解常微分方程
% Hfun为描述一阶微分方程导数的函数句柄,格式为 dX = Hfun( t,X )
% t为时间节点,可为标量,时间范围为 T = 0:h:t
% 长2向量 ,时间范围为 T = t(1):h:t(2)
% 向量 ,时间范围为 T = t
% h为时间步长,在t的前两种情况下,必须给出h具体值
% 直接给出时间节点t时,h可用[]或任意值占位
% x0为状态量初始值
% 算法:
% X(k) = X(k-1) + h*dX(k-1)
% By ZFS@wust 2021
% 获取更多Matlab/Simulink原创资料和程序,清关注微信公众号:Matlab Fans
function [T,X,dX] = ODE_ImplicitEuler( Hfun,t,h,x0 )
% [T,X,dX] = ODE_ImplicitEuler( Hfun,t,h,x0 ) 隐式欧拉法求解常微分方程
% Hfun为描述一阶微分方程导数的函数句柄,格式为 dX = Hfun( t,X )
% t为时间节点,可为标量,时间范围为 T = 0:h:t
% 长2向量 ,时间范围为 T = t(1):h:t(2)
% 向量 ,时间范围为 T = t
% h为时间步长,在t的前两种情况下,必须给出h具体值
% 直接给出时间节点t时,h可用[]或任意值占位
% x0为状态量初始值
% 算法:
% 使用fzero函数求解隐式微分方程 X(k) = X(k-1) + h*dX(k)
function [T,X,dX] = ODE_2OrderEuler( Hfun,t,h,x0 )
% [T,X,dX] = ODE_2OrderEuler( Hfun,t,h,x0 ) 两步欧拉法求解常微分方程
% Hfun为描述一阶微分方程导数的函数句柄,格式为 dX = Hfun( t,X )
% t为时间节点,可为标量,时间范围为 T = 0:h:t
% 长2向量 ,时间范围为 T = t(1):h:t(2)
% 向量 ,时间范围为 T = t
% h为时间步长,在t的前两种情况下,必须给出h具体值
% 直接给出时间节点t时,h可用[]或任意值占位
% x0为状态量初始值
% 算法:
% X(k) = X(k-2) + 2*h*dX(k-1)
function [T,X,dX] = ODE_ImprovedEuler( Hfun,t,h,x0 )
% [T,X,dX] = ODE_ImprovedEuler( Hfun,t,h,x0 ) 改进欧拉法求解常微分方程
% Hfun为描述一阶微分方程导数的函数句柄,格式为 dX = Hfun( t,X )
% t为时间节点,可为标量,时间范围为 T = 0:h:t
% 长2向量 ,时间范围为 T = t(1):h:t(2)
% 向量 ,时间范围为 T = t
% h为时间步长,在t的前两种情况下,必须给出h具体值
% 直接给出时间节点t时,h可用[]或任意值占位
% x0为状态量初始值
% 算法:
% Xp = X(k-1) + h*dX(k-1)
% dXp = Hfun( T(k),Xp )
% X(k) = X(k-1) + (h/2)*[dX(k-1)+dXp]
下面结合实例进行演示和分析。
4. 案例
求解一阶常微分方程(式中向量
x
{\bf{x}}
x等价于前文中的向量
y
{\bf{y}}
y):
x
˙
=
f
(
t
,
x
)
=
[
x
(
2
)
∗
x
(
3
)
−
x
(
1
)
∗
x
(
3
)
−
0.51
∗
x
(
1
)
∗
x
(
2
)
]
\mathbf{\dot{x}}=f\left( t,\mathbf{x} \right)=\left[ \begin{matrix} \mathbf{x}(2)*\mathbf{x}(3) \\ -\mathbf{x}(1)*\mathbf{x}(3) \\ -0.51*\mathbf{x}(1)*\mathbf{x}(2) \\ \end{matrix} \right]
x˙=f(t,x)=⎣⎡x(2)∗x(3)−x(1)∗x(3)−0.51∗x(1)∗x(2)⎦⎤
x
(
0
)
=
[
0
1
1
]
\mathbf{x}\left( 0 \right)=\left[ \begin{matrix} 0 \\ 1 \\ 1 \\ \end{matrix} \right]
x(0)=⎣⎡011⎦⎤
% 欧拉算法测试程序
clear
clc
close all
%% 仿真步长 h=0.01 时
Hfun = @(t,x) [ x(2)*x(3); -x(1)*x(3); -0.51*x(1)*x(2)]; % 一阶微分方程导数表达式
% 参数
t = [0 12]; % 时间范围
h = 0.01; % 时间步长
x0 = [0 1 1]; % 初始状态
% 显式欧拉法求解
[T1,X1] = ODE_ExplicitEuler( Hfun,t,h,x0 );
% 隐式欧拉法求解
[T2,X2] = ODE_ImplicitEuler( Hfun,t,h,x0 );
% 两步欧拉法求解
[T3,X3] = ODE_2OrderEuler( Hfun,t,h,x0 );
% 改进欧拉法求解
[T4,X4] = ODE_ImprovedEuler( Hfun,t,h,x0 );
% Matlab自带ode45求解
[T5,X5] = ode45( Hfun,t,x0 );
% 绘图对比
figure
subplot(311)
plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1),T4,X4(:,1),T5,X5(:,1))
xlabel('Time(s)')
ylabel('x_1')
legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
subplot(312)
plot(T1,X1(:,2),T2,X2(:,2),T3,X3(:,2),T4,X4(:,2),T5,X5(:,2))
xlabel('Time(s)')
ylabel('x_2')
legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
subplot(313)
plot(T1,X1(:,3),T2,X2(:,3),T3,X3(:,3),T4,X4(:,3),T5,X5(:,3))
xlabel('Time(s)')
ylabel('x_3')
legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
不同时间步长 h h h时的数值计算结果:
- 步长
h
=
0.001
s
h=0.001s
h=0.001s
- 步长
h
=
0.01
s
h=0.01s
h=0.01s
- 步长
h
=
0.1
s
h=0.1s
h=0.1s
- 步长
h
=
0.5
s
h=0.5s
h=0.5s
- 步长
h
=
0.6
s
h=0.6s
h=0.6s
显式欧拉法 | 隐式欧拉法 | 两步欧拉法 | 改进欧拉法 | |
---|---|---|---|---|
h = 0.001 s h=0.001s h=0.001s | 精度高 | 精度高 | 精度高 | 精度高 |
h = 0.01 s h=0.01s h=0.01s | 精度不高 | 精度降低 | 精度高 | 精度高 |
h = 0.1 s h=0.1s h=0.1s | 精度很低 | 精度低 | 精度高 | 精度高 |
h = 0.5 s h=0.5s h=0.5s | 发散 | 发散 | 精度很低 | 精度不高 |
h = 0.6 s h=0.6s h=0.6s | 发散 | 发散 | 发散 | 精度低 |
h = 0.8 s h=0.8s h=0.8s | 发散 | 发散 | 发散 | 发散 |
步长很小时,四种欧拉算法都能求得精度较高的数值解;随着步长增大,显式欧拉法和隐式欧拉法这两种一阶算法精度下降很快,直至发散;四种欧拉算法中改进欧拉法精度最高、稳定性最好,推荐使用。特别注意,隐式欧拉法计算量较大且计算精度不高,不推荐使用。下一篇介绍计算精度更高的四阶龙格库塔法(Runge-Kutta)。
5. 联系作者
有Matlab/Simulink方面的技术问题,欢迎发送邮件至944077462@qq.com讨论。
更多Matlab/Simulink原创资料,欢迎关注微信公众号:Matlab Fans
源程序下载:
1. csdn资源: 欧拉法(Euler)求解常微分方程的Matlab程序及案例
2. 扫码关注微信公众号Matlab Fans,回复BK08获取百度网盘下载链接。