欧拉法(Euler)求解常微分方程的Matlab程序及案例

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(n1),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(n1))

写成向量形式:

[ 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(n1))

令:

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(n1)

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(n1))=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(k1) 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(k1)y(k)y(k1)=hy(k)y(k1)

微分方程变为隐式差分方程:

{ 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(k1)+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(k1),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(k1)y(k+1)y(k1)=2hy(k+1)y(k1)

微分方程变为显式差分方程:

{ 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(k1)+2hf(t(k),y(k))y(0)=y0

上式是关于 y ( k − 1 ) {\bf{y}}\left( k-1 \right) y(k1) 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(k1)+2hf(t(k),y(k))=[y(k1)+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(k1)+hf(t(k),y(k)) 相当于隐式欧拉法由 y ( k − 1 ) {\bf{y}}\left( k-1 \right) y(k1)求解 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.51x(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时的数值计算结果:

  1. 步长 h = 0.001 s h=0.001s h=0.001s
    步长h=0.001s
  2. 步长 h = 0.01 s h=0.01s h=0.01s
    步长h=0.01s
  3. 步长 h = 0.1 s h=0.1s h=0.1s
    步长h=0.1s
  4. 步长 h = 0.5 s h=0.5s h=0.5s
    步长h=0.5s
  5. 步长 h = 0.6 s h=0.6s 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获取百度网盘下载链接。

在这里插入图片描述

  • 77
    点赞
  • 652
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

MatlabFans_Mfun

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值