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

1. 概念

1) 常微分方程

自变量只有一个的微分方程，称为常微分方程；自变量数量2个或以上时，称为偏微分方程。
绝大多数实际工程问题，常微分方程的自变量都是时间t，通常表达为：

$F\left(t,y,\dot{y},\cdots ,y^{\left(n-1\right)},y^{\left(n\right)}\right)=0$

2) 一阶常微分方程

求解上述关于$y^{\left(n\right)}$的方程，得到$y^{\left(n\right)}$的表达式为：

$y^{\left(n\right)}=f\left(t,y,\dot{y},\cdots ,y^{\left(n-1\right)}\right)$

写成向量形式：

$\left[\begin{array}{c}\dot{y}\\ \ddot{y}\\ ⋮\\ y^{\left(n\right)}\end{array}\right]=\left[\begin{array}{c}\ddot{y}\\ y^{\left(3\right)}\\ ⋮\\ f\left(t,y,\dot{y},\cdots ,y^{\left(n-1\right)}\right)\end{array}\right]$

令：

$\mathbf{y}=\left[\begin{array}{c}y\\ \dot{y}\\ ⋮\\ y^{\left(n-1\right)}\end{array}\right]$

$f\left(t,\mathbf{y}\right)=\left[\begin{array}{c}\ddot{y}\\ y^{\left(3\right)}\\ ⋮\\ f\left(t,y,\dot{y},\cdots ,y^{\left(n-1\right)}\right)\end{array}\right]=\left[\begin{array}{c}\ddot{y}\\ y^{\left(3\right)}\\ ⋮\\ f\left(t,\mathbf{y}\right)\end{array}\right]$

得到向量形式的一阶常微分方程：
$\dot{\mathbf{y}}=f\left(t,\mathbf{y}\right)$

3) 微分方程的解析解法和数值解法

微分方程的求解方法有解析解法和数值解法，解析法是求出因变量$\mathbf{y}$
关于时间$t$的具体函数式，表达为$\mathbf{y}=\mathbf{y}\left(t\right)$；数值法是解出因变量$\mathbf{y}$关于时间$t$的离散序列，通常表达为 { t ( k ) , y ( k ) } \left\{ {t(k),{\bf{y}}(k)} \right\} {t(k),y(k)}。
绝大多数的非线性常微分方程，不存在或难以求出解析解，大多数情况下只能求取微分方程的数值解。

2. 算法

欧拉法(Euler)是求解一阶常微分方程初值问题的数值方法，分为显式欧拉法、隐式欧拉法、两步欧拉法和改进欧拉法。

2.1 显式欧拉法

一阶微分方程初始问题：

$\{\begin{array}{l}\dot{\mathbf{y}}=f\left(t,\mathbf{y}\right)\\ \mathbf{y}\left(t_0\right)={\mathbf{y}}_0\end{array}$

式中， $t_0$为初始时间（已知常数），${\mathbf{y}}_0$为初始状态（已知向量），$f\left(t,\mathbf{y}\right)$为关于时间$t$和状态$\mathbf{y}$的函数（已知函数）。
使用一阶向前差商代替微分，即：

$\dot{\mathbf{y}}\approx \frac{\mathbf{y}\left(k+1\right)-\mathbf{y}\left(k\right)}{t\left(k+1\right)-t\left(k\right)}=\frac{\mathbf{y}\left(k+1\right)-\mathbf{y}\left(k\right)}{h}$

式中，h为时间步长。
微分方程变为显式差分方程：

$\{\begin{array}{l}\mathbf{y}\left(k+1\right)=\mathbf{y}\left(k\right)+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)\\ \mathbf{y}\left(0\right)={\mathbf{y}}_0\end{array}$

上式是关于$\mathbf{y}\left(k-1\right)$向$\mathbf{y}\left(k\right)$的递推形式，可以根据初始条件按照递推关系依次求出$\mathbf{y}\left(1\right),\mathbf{y}\left(2\right),\mathbf{y}\left(3\right),\mathbf{y}\left(4\right)\cdots ,\mathbf{y}\left(N\right)\cdots$，此离散序列即为微分方程的数值解。

2.2 隐式欧拉法

使用一阶向后差商代替微分，即

$\dot{\mathbf{y}}\approx \frac{\mathbf{y}\left(k\right)-\mathbf{y}\left(k-1\right)}{t\left(k\right)-t\left(k-1\right)}=\frac{\mathbf{y}\left(k\right)-\mathbf{y}\left(k-1\right)}{h}$

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

$\{\begin{array}{l}\mathbf{y}\left(k\right)=\mathbf{y}\left(k-1\right)+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)\\ \mathbf{y}\left(0\right)={\mathbf{y}}_0\end{array}$

在第$k$步迭代时， $\mathbf{y}\left(k-1\right),t(k)$已知， $\mathbf{y}\left(k\right)$为待求未知量， $f\left(t\left(k\right),\mathbf{y}\left(k\right)\right)$是关于待求未知量$\mathbf{y}\left(k\right)$的函数，即方程为关于$\mathbf{y}\left(k\right)$的非线性方程。通过非线性方程的迭代解法（如牛顿迭代法）可求解出$\mathbf{y}\left(k\right)$。依次求出$\mathbf{y}\left(1\right),\mathbf{y}\left(2\right),\mathbf{y}\left(3\right),\mathbf{y}\left(4\right)\cdots ,\mathbf{y}\left(N\right)\cdots$，此离散序列即为微分方程的数值解。

2.3 两步欧拉法

使用中心差商代替微分：

$\dot{\mathbf{y}}\approx \frac{\mathbf{y}\left(k+1\right)-\mathbf{y}\left(k-1\right)}{t\left(k+1\right)-t\left(k-1\right)}=\frac{\mathbf{y}\left(k+1\right)-\mathbf{y}\left(k-1\right)}{2h}$

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

$\{\begin{array}{l}\mathbf{y}\left(k+1\right)=\mathbf{y}\left(k-1\right)+2hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)\\ \mathbf{y}\left(0\right)={\mathbf{y}}_0\end{array}$

上式是关于$\mathbf{y}\left(k-1\right)$，$\mathbf{y}\left(k\right)$向$\mathbf{y}\left(k+1\right)$的递推形式，可以根据初始条件按照递推关系依次求出$\mathbf{y}\left(1\right),\mathbf{y}\left(2\right),\mathbf{y}\left(3\right),\mathbf{y}\left(4\right)\cdots ,\mathbf{y}\left(N\right)\cdots$，此离散序列即为微分方程的数值解。

这里需要特别注意$\mathbf{y}\left(1\right)$的计算方法：按照公式计算$\mathbf{y}\left(1\right)$需要用到无意义的$\mathbf{y}\left(-1\right)$，此步需要转换成显式欧拉法计算：

$\{\begin{array}{l}\mathbf{y}\left(1\right)=\mathbf{y}\left(0\right)+hf\left(t\left(0\right),\mathbf{y}\left(0\right)\right)\\ \mathbf{y}\left(0\right)={\mathbf{y}}_0\end{array}$

进一步分析：

$\mathbf{y}\left(k+1\right)=\mathbf{y}\left(k-1\right)+2hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)=\left[\mathbf{y}\left(k-1\right)+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)\right]+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)$

式中， $\bar{\mathbf{y}}\left(k\right)=\mathbf{y}\left(k-1\right)+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)$ 相当于隐式欧拉法由$\mathbf{y}\left(k-1\right)$求解$\mathbf{y}\left(k\right)$； $\mathbf{y}\left(k+1\right)=\bar{\mathbf{y}}\left(k\right)+f\left(t\left(k\right),\mathbf{y}\left(k\right)\right)$相当于显式欧拉法由$\mathbf{y}\left(k\right)$求解$\mathbf{y}\left(k+1\right)$。两步欧拉法是隐式欧拉法与显式欧拉法相结合的一种算法。

2.4 改进欧拉法

先使用显式欧拉法求解$k+1$步的预测值：

$\bar{\mathbf{y}}\left(k+1\right)=\mathbf{y}\left(k\right)+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)$

再使用$k$步的微分值$f\left(t\left(k\right),\mathbf{y}\left(k\right)\right)$和$k+1$步预测的微分值$f\left(t\left(k+1\right),\bar{\mathbf{y}}\left(k+1\right)\right)$的平均值来近似微分：

$\dot{\mathbf{y}}\approx \frac{f\left(t\left(k\right),\mathbf{y}\left(k\right)\right)+f\left(t\left(k+1\right),\bar{\mathbf{y}}\left(k+1\right)\right)}{2}$

再利用显式欧拉法校正$k+1$步的计算值：

$\mathbf{y}\left(k+1\right)=\mathbf{y}\left(k\right)+h\frac{f\left(t\left(k\right),\mathbf{y}\left(k\right)\right)+f\left(t\left(k+1\right),\bar{\mathbf{y}}\left(k+1\right)\right)}{2}$

综上所述，改进欧拉法的算法为：

$\{\begin{array}{l}\bar{\mathbf{y}}\left(k+1\right)=\mathbf{y}\left(k\right)+hf\left(t\left(k\right),\mathbf{y}\left(k\right)\right)\\ \mathbf{y}\left(k+1\right)=\mathbf{y}\left(k\right)+h\frac{f\left(t\left(k\right),\mathbf{y}\left(k\right)\right)+f\left(t\left(k+1\right),\bar{\mathbf{y}}\left(k+1\right)\right)}{2}\\ \mathbf{y}\left(0\right)={\mathbf{y}}_0\end{array}$

上式是关于$\mathbf{y}\left(k\right)$向$\mathbf{y}\left(k+1\right)$的递推形式，可以根据初始条件按照递推关系依次求出$\mathbf{y}\left(1\right),\mathbf{y}\left(2\right),\mathbf{y}\left(3\right),\mathbf{y}\left(4\right)\cdots ,\mathbf{y}\left(N\right)\cdots$，此离散序列即为微分方程的数值解。

2.5 四种欧拉方法的对比

微分方程的数值解法，本质是使用数值积分来实现$\dot{\mathbf{y}}$向$\mathbf{y}$的转换。显式欧拉法和隐式欧拉法，使用的是矩形数值积分方法，矩形积分是一阶算法，其截断误差为$O\left(h^2\right)$；两步欧拉法和改进欧拉法，使用的是梯形数值积分方法，梯形积分是二阶算法，其截断误差为$O\left(h^3\right)$。

算法类型	计算量	计算精度	稳定性
显式欧拉法	小	低	差
隐式欧拉法	大	低	差
两步欧拉法	中	较高	中
改进欧拉法	中	高	好

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. 案例

求解一阶常微分方程（式中向量$\mathbf{x}$等价于前文中的向量$\mathbf{y}$）：
$\dot{\mathbf{x}}=f\left(t,\mathbf{x}\right)=\left[\begin{array}{c}\mathbf{x}(2)\ast \mathbf{x}(3)\\ -\mathbf{x}(1)\ast \mathbf{x}(3)\\ -0.51\ast \mathbf{x}(1)\ast \mathbf{x}(2)\end{array}\right]$
$\mathbf{x}\left(0\right)=\left[\begin{array}{c}0\\ 1\\ 1\end{array}\right]$

% 欧拉算法测试程序
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$时的数值计算结果：

1. 步长$h=0.001s$
   
2. 步长$h=0.01s$
   
3. 步长$h=0.1s$
   
4. 步长$h=0.5s$
   
5. 步长$h=0.6s$
   

	显式欧拉法	隐式欧拉法	两步欧拉法	改进欧拉法
$h=0.001s$	精度高	精度高	精度高	精度高
$h=0.01s$	精度不高	精度降低	精度高	精度高
$h=0.1s$	精度很低	精度低	精度高	精度高
$h=0.5s$	发散	发散	精度很低	精度不高
$h=0.6s$	发散	发散	发散	精度低
$h=0.8s$	发散	发散	发散	发散

步长很小时，四种欧拉算法都能求得精度较高的数值解；随着步长增大，显式欧拉法和隐式欧拉法这两种一阶算法精度下降很快，直至发散；四种欧拉算法中改进欧拉法精度最高、稳定性最好，推荐使用。特别注意，隐式欧拉法计算量较大且计算精度不高，不推荐使用。下一篇介绍计算精度更高的四阶龙格库塔法(Runge-Kutta)。

5. 联系作者

有Matlab/Simulink方面的技术问题，欢迎发送邮件至944077462@qq.com讨论。
更多Matlab/Simulink原创资料，欢迎关注微信公众号：Matlab Fans

源程序下载:
1. csdn资源: 欧拉法(Euler)求解常微分方程的Matlab程序及案例
2. 扫码关注微信公众号Matlab Fans，回复BK08获取百度网盘下载链接。