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

19 篇文章 63 订阅

### 1. 概念

#### 1) 常微分方程

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

#### 2) 一阶常微分方程

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 ˙ 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 ˙ ⋮ y ( n − 1 ) ] \mathbf{y}=\left[ \begin{matrix} y \\ {\dot{y}} \\ \vdots \\ {{y}^{\left( n-1 \right)}} \\ \end{matrix} \right]

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]

y ˙ = f ( t , y ) \mathbf{\dot{y}}=f\left( t,\mathbf{y} \right)

### 2. 算法

#### 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 ˙ ≈ 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 ( 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.

#### 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 ( 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.

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

#### 2.4 改进欧拉法

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 ˙ ≈ 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 ( 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 ) + 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.

### 3. 程序

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

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 ˙ = 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 ( 0 ) = [ 0 1 1 ] \mathbf{x}\left( 0 \right)=\left[ \begin{matrix} 0 \\ 1 \\ 1 \\ \end{matrix} \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')


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

h = 0.001 s h=0.001s 精度高精度高精度高精度高
h = 0.01 s h=0.01s 精度不高精度降低精度高精度高
h = 0.1 s h=0.1s 精度很低精度低精度高精度高
h = 0.5 s h=0.5s 发散发散精度很低精度不高
h = 0.6 s h=0.6s 发散发散发散精度低
h = 0.8 s h=0.8s 发散发散发散发散

### 5. 联系作者

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

• 77
点赞
• 652
收藏
觉得还不错? 一键收藏
• 打赏
• 3
评论
01-16 1万+
10-23
06-11 276
11-29 4508
03-16 8318
12-17 1673
07-06

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

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

MatlabFans_Mfun

¥1 ¥2 ¥4 ¥6 ¥10 ¥20

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