线性常系数微分方程（组）

y = d s o l v e ( f 1 , f 2 , … , f m ) % 默认自变量为 t y=dsolve(f_1,f_2,\dots,f_m)\quad\%默认自变量为t y=dsolve(f1,f2,…,fm)%默认自变量为t
y = d s o l v e ( f 1 , f 2 , … , f m , ′ x ′ ) % 指定自变量为 x y=dsolve(f_1,f_2,\dots,f_m,'x')\quad\%指定自变量为x y=dsolve(f1,f2,…,fm,′x′)%指定自变量为x
其中， f i 可以描述微分方程，初始条件或边界条件，可以用字符串（推荐）或符号表达式描述其中，f_i可以描述微分方程，初始条件或边界条件，可以用字符串（推荐）或符号表达式描述其中，fi可以描述微分方程，初始条件或边界条件，可以用字符串（推荐）或符号表达式描述
可以求出方程的解析解可以求出方程的解析解可以求出方程的解析解

例1

u ( t ) = e − 5 t c o s ( 2 t + 1 ) + 5 , 解方程 u(t)=e^{-5t}cos(2t+1)+5,解方程 u(t)=e−5tcos(2t+1)+5,解方程
y ( 4 ) ( t ) + 10 y ′ ′ ′ ( t ) + 35 y ′ ′ ( t ) + 50 y ′ ( t ) + 24 y ( t ) = 5 u ′ ′ ( t ) + 4 u ′ ( t ) + 2 u ( t ) y^{(4)}(t)+10y'''(t)+35y''(t)+50y'(t)+24y(t)=5u''(t)+4u'(t)+2u(t) y(4)(t)+10y′′′(t)+35y′′(t)+50y′(t)+24y(t)=5u′′(t)+4u′(t)+2u(t)

syms t;
u=exp(-5*t)*cos(2*t+1)+5;
uu=diff(u,t,2)+diff(u,t)+2*u;
%用字符串描述
y1=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',char(uu)]);%用字符串描述的另一种方法
y2=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=5*diff(u,t,2)+4*diff(u,t)+2*u(t)']);

假设已知初值， y ( 0 ) = 3 , y ′ ( 0 ) = 2 , y ′ ′ ( 0 ) = y ′ ′ ′ ( 0 ) = 0 假设已知初值，y(0)=3,y'(0)=2,y''(0)=y'''(0)=0 假设已知初值，y(0)=3,y′(0)=2,y′′(0)=y′′′(0)=0

syms t;
u=exp(-5*t)*cos(2*t+1)+5;
uu=diff(u,t,2)+diff(u,t)+2*u;
%用字符串描述
y3=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',char(uu)],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0');
y4=vpa(simplify(y3),3);
ezplot(y4,[0 10]);

方程组同理，线性状态空间方程也可以使用 d s o l v e 求解，某些极特殊的非线性微分方程也可以使用 d s o l v e 求解方程组同理，线性状态空间方程也可以使用dsolve求解，某些极特殊的非线性微分方程也可以使用dsolve求解方程组同理，线性状态空间方程也可以使用dsolve求解，某些极特殊的非线性微分方程也可以使用dsolve求解

微分方程的数值解

一阶显式微分方程（组）

一阶显式微分方程 : x ′ ( t ) = f ( t , x ( t ) ) 一阶显式微分方程:\quad \textbf{x}'(t)=\textbf{f}(t,\textbf{x}(t)) 一阶显式微分方程:x′(t)=f(t,x(t))
[ t , x ] = o d e 45 ( F u n , [ t 0 , t n ] , x 0 ) % 直接求解 [t,x]=ode45(Fun,[t_0,t_n],x_0)\quad\%直接求解 [t,x]=ode45(Fun,[t0,tn],x0)%直接求解
[ t , x ] = o d e 45 ( F u n , [ t 0 , t n ] , x 0 , o p t i o n s ) 带有控制选项 [t,x]=ode45(Fun,[t_0,t_n],x_0,options)\quad带有控制选项 [t,x]=ode45(Fun,[t0,tn],x0,options)带有控制选项
[ t , x ] = o d e 45 ( F u n , [ t 0 , t n ] , x 0 , o p t i o n s , p 1 , p 2 , … , p m ) % 带有附加参数 [t,x]=ode45(Fun,[t_0,t_n],x_0,options,p_1,p_2,\dots,p_m)\quad\%带有附加参数 [t,x]=ode45(Fun,[t0,tn],x0,options,p1,p2,…,pm)%带有附加参数
F u n 使用函数或匿名函数（推荐）方式描述， [ t 0 , t n ] 是求解区间， x 0 是初始状态变量 Fun使用函数或匿名函数（推荐）方式描述，[t_0,t_n]是求解区间，x_0是初始状态变量 Fun使用函数或匿名函数（推荐）方式描述，[t0,tn]是求解区间，x0是初始状态变量
关于控制选项的设定关于控制选项的设定关于控制选项的设定

参数名	参数说明
RelTol	相对误差容许上限，默认0.001
AbsTol	一个向量，分量表示每个状态变量允许的绝对误差，默认1e-6
MaxStep	最大允许步长
Mass	质量矩阵，用于微分代数方程

opt=odeset;
opt.RelTol=1e-7;

例2

{ x 1 ′ ( t ) = − β x 1 ( t ) + x 2 ( t ) x 3 ( t ) x 2 ′ ( t ) = − ρ x 2 ( t ) + ρ x 3 ( t ) x 3 ′ ( t ) = − x 1 ( t ) x 2 ( t ) + σ x 2 ( t ) − x 3 ( t ) \begin{cases} x'_1(t)=-\beta x_1(t)+x_2(t)x_3(t)\\ x'_2(t)=-\rho x_2(t)+\rho x_3(t) \\ x'_3(t)=-x_1(t)x_2(t)+\sigma x_2(t)-x_3(t) \end{cases} ⎩⎪⎨⎪⎧x1′(t)=−βx1(t)+x2(t)x3(t)x2′(t)=−ρx2(t)+ρx3(t)x3′(t)=−x1(t)x2(t)+σx2(t)−x3(t)
β = 2 , ρ = 5 , σ = 20 , x 1 ( 0 ) = 0 , x 2 ( 0 ) = 0 , x 3 ( 0 ) = 1 0 − 10 ，求解方程组 \beta=2,\rho=5,\sigma=20,x_1(0)=0,x_2(0)=0,x_3(0)=10^{-10}，求解方程组 β=2,ρ=5,σ=20,x1(0)=0,x2(0)=0,x3(0)=10−10，求解方程组

%使用匿名函数描述
f=@(t,x,beta,rho,sigma)[-beta*x(1)+x(2)*x(3) ; -rho*x(2)+rho*x(3) ; -x(1)*x(2)+sigma*x(2)-x(3)];
beta=2;rho=5;sigma=20;x0=[0 0 1e-10];
[t,x]=ode45(f,[0 20],x0,[],beta,rho,sigma);
subplot(221);plot(t,x);grid on;
subplot(222);plot3(x(:,1),x(:,2),x(:,3));grid on;%相空间轨迹
subplot(223);comet3(x(:,1),x(:,2),x(:,3));grid on;%动态相空间轨迹%使用M函数描述
function dx=lorenz1(t,x,beta,rho,sigma)
dx=[-beta*x(1)+x(2)*x(3) ; -rho*x(2)+rho*x(3) ; -x(1)*x(2)+sigma*x(2)-x(3)];
beta=2;rho=5;sigma=20;x0=[0 0 1e-10];
[t,x]=ode45(@lorenz1,[0 20],x0,[],beta,rho,sigma);
subplot(221);plot(t,x);grid on;
subplot(222);plot3(x(:,1),x(:,2),x(:,3));grid on;%相空间轨迹
subplot(223);comet3(x(:,1),x(:,2),x(:,3));grid on;%动态相空间轨迹

数值解的验证：将 R e l T o l 或 A b s T o l 设置为更小的值，观察解是否一致数值解的验证：将RelTol或AbsTol设置为更小的值，观察解是否一致数值解的验证：将RelTol或AbsTol设置为更小的值，观察解是否一致

高阶常微分方程（组）

高阶常微分方程： y ( n ) = f ( t , y , y ′ , … , y ( n − 1 ) ) 高阶常微分方程：y^{(n)}=f(t,y,y',\dots,y^{(n-1)}) 高阶常微分方程：y(n)=f(t,y,y′,…,y(n−1))
转换方法：选择一组状态变量 x 1 = y , x 2 = y ′ , … , x n = y ( n − 1 ) , 得到转换方法：选择一组状态变量x_1=y,x_2=y',\dots,x_n=y^{(n-1)},得到转换方法：选择一组状态变量x1=y,x2=y′,…,xn=y(n−1),得到
{ x 1 ′ = x 2 x 2 ′ = x 3 ⋮ x n ′ = f ( t , x 1 , x 2 , … , x n ) \begin{cases} x'_1=x_2\\ x'_2=x_3 \\ \quad \vdots \\ x'_n=f(t,x_1,x_2,\dots,x_n)\end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1′=x2x2′=x3⋮xn′=f(t,x1,x2,…,xn)
x 1 ( 0 ) = y ( 0 ) , x 2 ( 0 ) = y ′ ( 0 ) , … , x n ( 0 ) = y ( n − 1 ) ( 0 ) x_1(0)=y(0),x_2(0)=y'(0),\dots,x_n(0)=y^{(n-1)}(0) x1(0)=y(0),x2(0)=y′(0),…,xn(0)=y(n−1)(0)

例3

{ x ′ ′ + 2 y ′ x = 2 y ′ ′ x ′ ′ y ′ + 2 x ′ y ′ ′ + x y ′ − y = 5 \begin{cases} x''+2y'x=2y'' \\ x''y'+2x'y''+xy'-y=5 \end{cases} {x′′+2y′x=2y′′x′′y′+2x′y′′+xy′−y=5
x ( 0 ) = x ′ ( 0 ) = y ( 0 ) = y ′ ( 0 ) = 1 x(0)=x'(0)=y(0)=y'(0)=1 x(0)=x′(0)=y(0)=y′(0)=1
求解过程：求解过程：求解过程：
x 1 = x , x 2 = x ′ , x 3 = y , x 4 = y ′ x_1=x,x_2=x',x_3=y,x_4=y' x1=x,x2=x′,x3=y,x4=y′

syms x(t) dx dy
[dx,dy]=solve(dx+2*x(4)*x(1)-2*dy,dx*x(4)+3*x(2)*dy+x(1)*x(4)-x(3)-5,[dx dy]);
%有的方程可以直接转换，有的需要算一下

{ x 1 ′ = x 2 x 2 ′ = d x x 3 ′ = x 4 x 4 ′ = d y \begin{cases} x_1'=x_2 \\ x_2'= dx \\ x_3'=x_4 \\ x_4'=dy \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧x1′=x2x2′=dxx3′=x4x4′=dy

f=@(t,x)[x(2) ; dx ; x(4) ; dy];%这里可能需要将dx,dy的值复制放到相应位置，
[t x]=ode45(f,[0 10],[1 1 1 1]);

刚性微分方程

[ t , x ] = o d e 15 s ( F u n , [ t 0 , t n ] , x 0 , o p t i o n s , p 1 , p 2 , … , p m ) [t,x]=ode15s(Fun,[t_0,t_n],x_0,options,p_1,p_2,\dots,p_m) [t,x]=ode15s(Fun,[t0,tn],x0,options,p1,p2,…,pm)
用法，解法同 o d e 45 用法，解法同ode45 用法，解法同ode45

隐式微分方程

隐式微分方程： F [ t , x ( t ) , x ′ ( t ) = 0 , x ( 0 ) = x 0 , x ′ ( 0 ) = x 0 ′ 隐式微分方程：\textbf{F}[t,x(t),x'(t)=0,x(0)=x_0,x'(0)=x_0' 隐式微分方程：F[t,x(t),x′(t)=0,x(0)=x0,x′(0)=x0′
[ x 0 ∗ , x 0 ′ ∗ ] = d e c i c ( f u n , t 0 , x 0 , x 0 F , x 0 ′ , x 0 ′ F ) % 获得相容初始条件 [x_0^*,x_0^{'*}]=decic(fun,t_0,x_0,x_0^F,x_0',x_0'^F)\quad\%获得相容初始条件 [x0∗,x0′∗]=decic(fun,t0,x0,x0F,x0′,x0′F)%获得相容初始条件
[ t , x ] = o d e 15 i ( f u n , t s p a n , x 0 ∗ , x 0 ′ ∗ ) [t,x]=ode15i(fun,tspan,x_0^*,x_0'^*) [t,x]=ode15i(fun,tspan,x0∗,x0′∗)
使用 d e c i c 函数有时可能有一些坑 , 可能是由于 M A T L A B 版本的问题使用decic函数有时可能有一些坑,可能是由于MATLAB版本的问题使用decic函数有时可能有一些坑,可能是由于MATLAB版本的问题
x 0 F , x 0 ′ F 为 n 维列向量，值为 1 表示需要保留的初值，值为 0 表示需要求解的初值 x_0^F,x_0'^F为n维列向量，值为1表示需要保留的初值，值为0表示需要求解的初值 x0F,x0′F为n维列向量，值为1表示需要保留的初值，值为0表示需要求解的初值
%尽量避免使用decic函数，避免出现某些意想不到的情况（大部分情况还是可以用的）

例4

{ x ′ ′ s i n y ′ + y ′ ′ 2 = − 2 x y x − x ′ + x x ′ ′ y ′ x x ′ ′ y ′ ′ + c o s y ′ ′ = 3 y x ′ e − x \begin{cases} x''siny'+y''^2=-2xyx^{-x'}+xx''y' \\ xx''y''+cosy''=3yx'e^{-x} \end{cases} {x′′siny′+y′′2=−2xyx−x′+xx′′y′xx′′y′′+cosy′′=3yx′e−x
x ( 0 ) = 1 , x ′ ( 0 ) = 0 , y ( 0 ) = 0 , y ′ ( 0 ) = 1 x(0)=1,x'(0)=0,y(0)=0,y'(0)=1 x(0)=1,x′(0)=0,y(0)=0,y′(0)=1
求解过程：求解过程：求解过程：
x 1 = x , x 2 = x ′ , x 3 = y , x 4 = y ′ ，转化为 x_1=x,x_2=x',x_3=y,x_4=y'，转化为 x1=x,x2=x′,x3=y,x4=y′，转化为
{ x 1 ′ − x 2 = 0 s 2 ′ s i n x 4 + x 4 ′ + 2 e − x 2 x 1 x 3 − x 1 x 2 ′ x 4 = 0 x 3 ′ − x 4 = 0 x 1 x 2 ′ x 4 ′ + c o s x 4 ′ − e e − x 1 x 3 x 2 = 0 \begin{cases} x_1'-x_2=0 \\ s_2'sinx_4+x_4'+2e^{-x_2}x_1x_3-x_1x_2'x_4=0 \\ x_3'-x_4=0 \\ x_1x_2'x_4'+cosx_4'-ee^{-x_1}x_3x_2=0 \end{cases} ⎩⎪⎪⎪⎨⎪⎪⎪⎧x1′−x2=0s2′sinx4+x4′+2e−x2x1x3−x1x2′x4=0x3′−x4=0x1x2′x4′+cosx4′−ee−x1x3x2=0

f=@(t,x,xd)[xd(1)-x(2);
xd(2)*sin(x(4))+xd(4)^2+2*exp(-x(2))*x(1)*x(3)-x(1)*xd(2)*x(4);
xd(3)-x(4);
x(1)*xd(2)*xd(4)+cos(xd(4))-3*exp(-x(1))*x(3)*x(2)];
%注意xd
x0=[1 0 0 1];xd0=[0 1 1 -1];x0F=[1 1 1 1];xd0F=[];
[X0,XD0]=decic(f,0,x0,x0F,xd0,xd0F);
[t x]=ode15i(f,[0 2],X0,XD0);
plot(t,x);

使用经验与困惑：使用经验与困惑：使用经验与困惑：

1 、注意 x d 0 初值的选取，即使 x d 0 F 相应分量为 0 1、注意xd0初值的选取，即使xd0F相应分量为0 1、注意xd0初值的选取，即使xd0F相应分量为0
x ( 0 ) = 1 , x ′ ( 0 ) = 0 , y ( 0 ) = 0 , y ′ ( 0 ) = 1 带入，得到 x(0)=1,x'(0)=0,y(0)=0,y'(0)=1带入，得到 x(0)=1,x′(0)=0,y(0)=0,y′(0)=1带入，得到
{ x ′ ′ s i n ( 1 ) + y ′ ′ = x ′ ′ x ′ ′ y ′ ′ + c o s y ′ ′ = 0 \begin{cases} x''sin(1)+y''=x'' \\ x''y''+cosy''=0 \end{cases} {x′′sin(1)+y′′=x′′x′′y′′+cosy′′=0
求解求解求解

%testfunc.m文件
function F=testfunc(x)
F(1)=x(1)*sin(1)+x(2)-x(1);
F(2)=x(1)*x(2)+cos(x(2));
end

X=fsolve(@testfunc,rand(2,1));
%求解得 X(-1.85 0.34) 或 X(1.85  -0.34)

注意下面这两种写法：注意下面这两种写法：注意下面这两种写法：

%写法一
f=@(t,x,xd)[xd(1)-x(2);
xd(2)*sin(x(4))+xd(4)^2+2*exp(-x(2))*x(1)*x(3)-x(1)*xd(2)*x(4);
xd(3)-x(4);
x(1)*xd(2)*xd(4)+cos(xd(4))-3*exp(-x(1))*x(3)*x(2)];
x0=[1 0 0 1];xd0=[0 -1.85 1 0.34];x0F=[1 1 1 1];xd0F=[0 0 0 0];
%注意这里的xd0
[X0,XD0]=decic(f,0,x0,x0F,xd0,xd0F);
[t x]=ode15i(f,[0 2],X0,XD0);
plot(t,x);%报错
错误使用 decic (line 108)
DECIC 中出现收敛错误。%写法二
f=@(t,x,xd)[xd(1)-x(2);
xd(2)*sin(x(4))+xd(4)^2+2*exp(-x(2))*x(1)*x(3)-x(1)*xd(2)*x(4);
xd(3)-x(4);
x(1)*xd(2)*xd(4)+cos(xd(4))-3*exp(-x(1))*x(3)*x(2)];
x0=[1 0 0 1];xd0=[0 1.85 1 -0.34];x0F=[1 1 1 1];xd0F=[0 0 0 0];
%注意这里的xd0
[X0,XD0]=decic(f,0,x0,x0F,xd0,xd0F);
[t x]=ode15i(f,[0 2],X0,XD0);
plot(t,x);
%正确求解~~

d e c i c 函数就是为了求解那些隐式中的状态变量的一阶导数 , 得到一致性条件 decic函数就是为了求解那些隐式中的状态变量的一阶导数,得到一致性条件 decic函数就是为了求解那些隐式中的状态变量的一阶导数,得到一致性条件

2 、解方程 x ( 0 ) = 0 , x ′ ( 0 ) = 1 , y ( 0 ) = 0 2、解方程 x(0)=0,x'(0)=1,y(0)=0 2、解方程x(0)=0,x′(0)=1,y(0)=0
{ s ′ ′ s i n y ′ + y ′ + x y = 0 x y ′ + x c o s y ′ − x ′ y = 0 \begin{cases} s''siny'+y'+xy=0 \\ xy'+xcosy'-x'y=0 \end{cases} {s′′siny′+y′+xy=0xy′+xcosy′−x′y=0

%写法一
f=@(t,x,xd)[xd(1)-x(2);
xd(2)*sin(xd(3))+xd(3)+x(1)*x(3);
x(1)*xd(3)+x(1)*cos(xd(3))-x(2)*x(3)];
x0=[0 1 0];x0F=[1 1 0];xd0F=[0 0 0];
xd0=[1;5*(rand(2,1)-0.5)];
[X0,XD0]=decic(f,0,x0,x0F,xd0,xd0F);
[t x]=ode15i(f,[0 2],X0,XD0);
plot(t,x);%报错
错误使用 decic>sls (line 128)
请尝试释放 1 个固定分量。出错 decic (line 77)[dy,dyp] =sls(res,dfdy,dfdyp,neq,free_y,free_yp); %写法二
f=@(t,x,xd)[xd(1)-x(2);
xd(2)*sin(xd(3))+xd(3)+x(1)*x(3);
x(1)*xd(3)+x(1)*cos(xd(3))-x(2)*x(3)];
x0=[0 1 0];x0F=[1 1 0];xd0F=[0 0 0];
%注意这里的x0F
xd0=[1;5*(rand(2,1)-0.5)];
[X0,XD0]=decic(f,0,x0,x0F,xd0,xd0F);
[t x]=ode15i(f,[0 2],X0,XD0);
plot(t,x);
警告: 在 t=1.160344e-01 处失败。在时
间 t 处，步长大小必须降至所允许的最小
值(4.122369e-16)以下，才能达到积分容差
要求。%写法三
f=@(t,x,xd)[xd(1)-x(2);
xd(2)*sin(xd(3))+xd(3)+x(1)*x(3);
x(1)*xd(3)+x(1)*cos(xd(3))-x(2)*x(3)];
x0=[0 1 0];
yy=@(y) y(1)*sin(y(2))+y(2);
y=fsolve(yy,rand(2,1));
xd0=[1;y];
[t x]=ode15i(f,[0 10],x0,xd0);
plot(t,x);
警告: 在 t=1.160346e-01 处失败。在时
间 t 处，步长大小必须降至所允许的最小
值(4.122377e-16)以下，才能达到积分容差
要求。

这个方程组解的我很迷这个方程组解的我很迷这个方程组解的我很迷
可以参考的链接：http://www.matlabsky.com/forum.php?mod=viewthread&tid=529&highlight=%D2%FE%CA%BD%CE%A2%B7%D6%B7%BD%B3%CC

3 、上面的警告该如何解决呢？顺便再附一个题目 3、上面的警告该如何解决呢？顺便再附一个题目 3、上面的警告该如何解决呢？顺便再附一个题目
{ x 1 ′ x 2 ′ ′ s i n ( x 1 x 2 ) + 5 x 1 ′ ′ x 2 ′ c o s ( x 1 2 ) + t 2 x 1 x 2 2 = e − x 2 2 x 1 ′ ′ x 2 + x 2 ′ ′ x 1 ′ s i n ( x 1 2 ) + c o s ( x 2 ′ ′ x 2 ) = s i n ( t ) \begin{cases} x_1'x_2''sin(x_1x_2)+5x_1''x_2'cos(x_1^2)+t^2x_1x_2^2=e^{-x_2^2} \\ x_1''x_2+x_2''x_1'sin(x_1^2)+cos(x_2''x_2)=sin(t) \end{cases} {x1′x2′′sin(x1x2)+5x1′′x2′cos(x12)+t2x1x22=e−x22x1′′x2+x2′′x1′sin(x12)+cos(x2′′x2)=sin(t)
x 1 ( 0 ) = 1 , x 1 ′ ( 0 ) = 1 , x 2 ( 0 ) = 2 , x 2 ′ ( 0 ) = 2 x_1(0)=1,x_1'(0)=1,x_2(0)=2,x_2'(0)=2 x1(0)=1,x1′(0)=1,x2(0)=2,x2′(0)=2

微分代数方程

微分代数方程： M ( t , x ) x ′ = f ( t , x ) % 即方程中某些变量间满足某些代数方程的约束微分代数方程：\quad M(t,x)x'=f(t,x) \quad \%即方程中某些变量间满足某些代数方程的约束微分代数方程：M(t,x)x′=f(t,x)%即方程中某些变量间满足某些代数方程的约束

例5

{ x 1 ′ = − 0.2 x 1 + x 2 x 3 + 0.3 x 1 x 2 x 2 ′ = 2 x 1 x 2 − 5 x 2 x 3 − 2 x 2 2 0 = x 1 + x 2 + x 3 − 1 \begin{cases} x_1'=-0.2x_1+x_2x_3+0.3x_1x_2 \\ x_2'=2x_1x_2-5x_2x_3-2x_2^2 \\ 0=x_1+x_2+x_3-1 \end{cases} ⎩⎪⎨⎪⎧x1′=−0.2x1+x2x3+0.3x1x2x2′=2x1x2−5x2x3−2x220=x1+x2+x3−1
x 1 ( 0 ) = 0.8 , x 2 ( 0 ) = x 3 ( 0 ) = 0.1 x_1(0)=0.8,x_2(0)=x_3(0)=0.1 x1(0)=0.8,x2(0)=x3(0)=0.1
求解过程：求解过程：求解过程：
[ 1 0 0 0 1 0 0 0 0 ] [ x 1 ′ x 2 ′ x 3 ′ ] = [ − 0.2 x 1 + x 2 x 3 + 0.3 x 1 x 2 2 x 1 x 2 − 5 x 2 x 3 − 2 x 2 2 x 1 + x 2 + x 3 − 1 ] \left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{matrix} \right] \left[ \begin{matrix} x_1' \\ x_2' \\ x_3' \end{matrix} \right] = \left[ \begin{matrix} -0.2x_1+x_2x_3+0.3x_1x_2 \\ 2x_1x_2-5x_2x_3-2x_2^2 \\ x_1+x_2+x_3-1 \end{matrix} \right] ⎣⎡100010000⎦⎤⎣⎡x1′x2′x3′⎦⎤=⎣⎡−0.2x1+x2x3+0.3x1x22x1x2−5x2x3−2x22x1+x2+x3−1⎦⎤

f=@(t,x)[-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2) ;
2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)^2 ;
x(1)+x(2)+x(3)-1];
M=diag([1 1 0]);
opt=odeset;opt.Mass=M;
x0=[0.8 0.1 0.1];
[t x]=ode15s(f,[0 10],x0,opt);
%重点在于对Mass矩阵的描述，求解函数可以选择不同的

例6

{ x 1 ′ s i n x 1 + x 2 ′ c o s x 2 + x 1 = 1 − x 1 ′ c o s x 2 + x 2 ′ s i n x 1 + x 2 = 0 \begin{cases} x_1'sinx_1+x_2'cosx_2+x_1=1 \\ -x_1'cosx_2+x_2'sinx_1+x_2=0 \end{cases} {x1′sinx1+x2′cosx2+x1=1−x1′cosx2+x2′sinx1+x2=0
x 1 ( 0 ) = x 2 ( 0 ) = 0 x_1(0)=x_2(0)=0 x1(0)=x2(0)=0
求解过程：求解过程：求解过程：
[ s i n x 1 c o s x 2 − c o s x 2 s i n x 1 ] [ x 1 ′ x 2 ′ ] = [ 1 − x 1 − x 2 ] \left[ \begin{matrix} sinx_1 & cosx_2 \\ -cosx_2 & sinx_1 \end{matrix} \right] \left[ \begin{matrix} x_1' \\ x_2' \end{matrix} \right] = \left[ \begin{matrix} 1-x_1 \\ -x_2 \end{matrix} \right] [sinx1−cosx2cosx2sinx1][x1′x2′]=[1−x1−x2]

f=@(t,x)[1-x(1);-x(2)];
M=@(t,x)[sin(x(1)),cos(x(2));-cos(x(2)),sin(x(1))];
opt=odeset;opt.Mass=M;
[t x]=ode45(f,[0 10],[0 0],opt);
plot(t,x);

延迟微分方程

典型延迟微分方程

典型延迟微分方程： x ′ ( t ) = f ( t , x ( t ) , x ( t − τ 1 ) , x ( t − τ 2 ) , … , x ( t − τ m ) ) 典型延迟微分方程：x'(t)=f(t,x(t),x(t-\tau_1),x(t-\tau_2),\dots,x(t-\tau_m)) 典型延迟微分方程：x′(t)=f(t,x(t),x(t−τ1),x(t−τ2),…,x(t−τm))
s o l = d d e 23 ( f 1 , τ , f 2 , [ t 0 , t n ] , o p t i o n s ) sol=dde23(f_1,\tau,f_2,[t_0,t_n],options) sol=dde23(f1,τ,f2,[t0,tn],options)

例6

历史函数 x 1 ( t ) = e 2.1 t , x 2 ( t ) = s i n ( t ) , x 3 ( t ) = c o s ( t ) 历史函数x_1(t)=e^{2.1t},x_2(t)=sin(t),x_3(t)=cos(t) 历史函数x1(t)=e2.1t,x2(t)=sin(t),x3(t)=cos(t)
{ x ′ ( t ) = 1 − 3 x ( t ) − y ( t − 1 ) − 0.2 x 3 ( t − 0.5 ) − x ( t − 0.5 ) y ′ ′ ( t ) + 3 y ′ ( t ) + 2 y ( t ) = 4 x ( t ) \begin{cases} x'(t)=1-3x(t)-y(t-1)-0.2x^3(t-0.5)-x(t-0.5) \\ y''(t)+3y'(t)+2y(t)=4x(t) \end{cases} {x′(t)=1−3x(t)−y(t−1)−0.2x3(t−0.5)−x(t−0.5)y′′(t)+3y′(t)+2y(t)=4x(t)

f=@(t,x,Z)[1-3*x(1)-Z(2,2)-0.2*Z(1,1)^3-Z(1,1) ; x(3) ; 4*x(1)-3*x(3)-2*x(2)];
f2=@(t,x)[exp(2.1*t);sin(t);cos(t)];
tau=[0.5 1];
sol=dde23(f,tau,f2,[0 100]);
plot(sol.x,sol.y);
t=linspace(1,100,1000);
X=deval(sol,t);%计算对应于t的状态变量x的取值，这个题目可能画第二个图没什么意义，为了防止忘记就写在这里
Y=deval(sol,t-1);
figure;plot(X,Y);

变时间延迟微分方程

s o l = d d e s d ( f , f τ , f 2 , [ t 0 , t n ] , o p t i o n s ) sol=ddesd(f,f_{\tau},f_2,[t_0,t_n],options) sol=ddesd(f,fτ,f2,[t0,tn],options)

例7

历史函数 : x 1 ( t ) = s i n ( t + 1 ) , x 2 ( t ) = c o s ( t ) , x 3 ( t ) = e 3 t , 当 t < 0 时， x = 0 ; 当 t = 0 时， x 满足历史函数历史函数:x_1(t)=sin(t+1),x_2(t)=cos(t),x_3(t)=e^{3t},当t<0时，\textbf{x}=0;当t=0时，\textbf{x}满足历史函数历史函数:x1(t)=sin(t+1),x2(t)=cos(t),x3(t)=e3t,当t<0时，x=0;当t=0时，x满足历史函数
{ x 1 ′ ( t ) = − 2 x 2 ( t ) − 3 x 1 ( t − 0.2 ∣ s i n ( t ) ∣ ) x 2 ′ ( t ) = − 0.05 x 1 ( t ) x 3 ( t ) − 2 x 2 ( t − 0.8 ) + 2 x 3 ′ ( t ) = 0.3 x 1 ( t ) x 2 ( t ) x 3 ( t ) + c o s ( x 1 ( t ) x 2 ( t ) ) + 2 s i n ( 0.1 t 2 ) \begin{cases} x_1'(t)=-2x_2(t)-3x_1(t-0.2|sin(t)|) \\ x_2'(t)=-0.05x_1(t)x_3(t)-2x_2(t-0.8)+2 \\ x_3'(t)=0.3x_1(t)x_2(t)x_3(t)+cos(x_1(t)x_2(t))+2sin(0.1t^2) \end{cases} ⎩⎪⎨⎪⎧x1′(t)=−2x2(t)−3x1(t−0.2∣sin(t)∣)x2′(t)=−0.05x1(t)x3(t)−2x2(t−0.8)+2x3′(t)=0.3x1(t)x2(t)x3(t)+cos(x1(t)x2(t))+2sin(0.1t2)

f=@(t,x,Z)[-2*x(2)-3*Z(1,1) ;
-0.05*x(1)*x(3)-2*Z(2,2)+2 ;
0.3*x(1)*x(2)*x(3)+cos(x(1)*x(2))+2*sin(0.1*t^2)];
f_tau=@(t,x)[t-0.2*abs(sin(t));t-0.8];
f2=@(t,x)[sin(t+1);cos(t);exp(3*t)]*(t==0);
sol=ddesd(f,f_tau,f2,[0 10]);
plot(sol.x,sol.y);

中立型延迟微分方程

中立型延迟微分方程： x ′ ( t ) = f ( t ， x ( t ) , x ( t − τ p 1 ) , … , x ( t − t τ m ) , x ′ ( t − τ q 1 ) , … , x ′ ( t − τ q k ) ) 中立型延迟微分方程：x'(t)=f(t，x(t),x(t-\tau_{p_1}),\dots,x(t-t_{\tau_m}),x'(t-\tau_{q_1}),\dots,x'(t-\tau_{q_k})) 中立型延迟微分方程：x′(t)=f(t，x(t),x(t−τp1),…,x(t−tτm),x′(t−τq1),…,x′(t−τqk))
s o l = d d e n s d ( f , τ 1 , τ 2 , f 2 , [ t 0 , t n ] , o p t i o n s ) sol=ddensd(f,\tau_1,\tau_2,f_2,[t_0,t_n],options) sol=ddensd(f,τ1,τ2,f2,[t0,tn],options)

再附一个练习：当 t ≤ 0 时， x ( t ) = t , y ( t ) = e t ，求解方程再附一个练习：当t\le0时，x(t)=t,y(t)=e^t，求解方程再附一个练习：当t≤0时，x(t)=t,y(t)=et，求解方程
{ x ′ ( t ) = x 2 ( t − 0.2 ) + y 2 ( t − 0.2 ) − 6 x ( t − 0.5 ) − 8 y ( t − 0.1 ) y ′ ( t ) = x ( t ) [ 2 y ( t − 0.2 ) − x ( t ) + 5 − 2 t 2 ( t − 0.1 ) ] \begin{cases} x'(t)=x^2(t-0.2)+y^2(t-0.2)-6x(t-0.5)-8y(t-0.1) \\ y'(t)=x(t)[2y(t-0.2)-x(t)+5-2t^2(t-0.1)] \end{cases} {x′(t)=x2(t−0.2)+y2(t−0.2)−6x(t−0.5)−8y(t−0.1)y′(t)=x(t)[2y(t−0.2)−x(t)+5−2t2(t−0.1)]
如果考虑方程最后一项 x 2 ( t − 1 ) 变成 x ′ ( t − 0.1 ) ，重新求解如果考虑方程最后一项x^2(t-1)变成x'(t-0.1)，重新求解如果考虑方程最后一项x2(t−1)变成x′(t−0.1)，重新求解
（这个题目的 t s p a n 最好设置为（这个题目的tspan最好设置为（这个题目的tspan最好设置为 [ 0 1 ] \left[ \begin{matrix} 0 & 1 \end{matrix} \right] [01] , 在 t = 1.5 时，已经达到了 1 0 25 数量级） ,在t=1.5时，已经达到了10^{25}数量级） ,在t=1.5时，已经达到了1025数量级）

边值问题的求解

s i n i t = b v p i n i t ( v , x 0 , θ 0 ) sinit=bvpinit(v,x_0,\theta_0) sinit=bvpinit(v,x0,θ0)
s o l = b v p 5 c ( @ f u n 1 , @ f u n 2 , s i n i t , o p t i o n s , p 1 , p 2 , … , p m ) sol=bvp5c(@fun1,@fun2,sinit,options,p_1,p_2,\dots,p_m) sol=bvp5c(@fun1,@fun2,sinit,options,p1,p2,…,pm)

例8

− u ′ ′ ( x ) + 6 u ( x ) = e 10 x c o s ( 12 x ) -u''(x)+6u(x)=e^{10x}cos(12x) −u′′(x)+6u(x)=e10xcos(12x)
u ( 0 ) = u ( 1 ) = 1 ，解方程 u(0)=u(1)=1，解方程 u(0)=u(1)=1，解方程

f=@(t,x)[x(2);6*x(1)-exp(10*t)*cos(12*t)];
f2=@(xa,xb)[xa(1)-1;xb(1)-1];
sinit=bvpinit(linspace(0,2,100),rand(2,1));
sol=bvp5c(f,f2,sinit);
plot(sol.x,sol.y);

例9

{ x ′ = 4 x − α x y y ′ = − 2 y + β x y \begin{cases} x'=4x-\alpha xy \\ y'=-2y+\beta xy \end{cases} {x′=4x−αxyy′=−2y+βxy
α , β 是参数， x ( 0 ) = 2 ， y ( 0 ) = 1 ， x ( 3 ) = 4 ， y ( 3 ) = 2 ，求解方程和 α , β \alpha,\beta是参数，x(0)=2，y(0)=1，x(3)=4，y(3)=2，求解方程和\alpha,\beta α,β是参数，x(0)=2，y(0)=1，x(3)=4，y(3)=2，求解方程和α,β

f=@(t,x,v)[4*x(1)-v(1)*x(1)*x(2);-2*x(2)+v(2)*x(1)*x(2)];
f2=@(xa,xb,v)[xa(1)-2;xa(2)-1;xb(1)-4;xb(2)-2];
sinit=bvpinit(linspace(0,5,100),5*(rand(2,1)-0.5),5*(rand(2,1)-0.5));
%求解过程中状态变量初值x0和待定参数theta0（如果有的话）的初始搜索点入锅选取不当，出现Jacobi矩阵奇异无法求解的情况，变换不同的值多试几次
sol=bvp5c(f,f2,sinit);
plot(sol.x,sol.y);
sol.parameters

再附一个练习题再附一个练习题再附一个练习题
x ′ ′ + 1 t x ′ + ( 1 − 1 4 t 2 ) x = t c o s ( t ) , x ( 1 ) = 1 , x ( 6 ) = − 0.5 x''+\frac{1}{t}x'+(1-\frac{1}{4t^2})x=\sqrt{t} cos(t),x(1)=1,x(6)=-0.5 x′′+t1x′+(1−4t21)x=t cos(t),x(1)=1,x(6)=−0.5

关于微分方程的更多了解

拉普拉斯变换拉氏变换求解线性常微分方程
偏微分方程
Simulink的微分方程框图求解
更多，欢迎补充

说明

以上的程序均在MATLAB 2019a测试通过
上面留下的练习题或多或少都有点坑，有的我解决了，有的还没，欢迎讨论
主要参考教程：薛定宇《高等应用数学问题的MATLAB求解 (第四版)》清华大学出版社
有的问题我也一直没有解决，欢迎讨论

MATLAB微分方程学习总结相关推荐

系统辨识理论及MATLAB仿真——学习笔记（1）
系统辨识理论及MATLAB仿真学习笔记(1) 前言目录第1章绪论 1.1 建立数学模型的基本方法 1.2 系统辨识的定义 1.3 系统辨识的研究目的 1.4 数学模型的分类 1.5 几种常见的数 ...
matlab强化学习算例理/菜鸟理解1——双足机器人行走算例
目录 matlab双足机器人强化学习算例介绍强化学习的一些基础理解菜鸟对一些名词的理解 matlab强化学习库介绍双足机器人算例逻辑盘点如何改写算例做自己的强化学习. %写在前面: 本人大四狗 ...
【MATLAB深度学习工具箱】学习笔记--体脂估计算例再分析：拟合神经网络fitnet里面的数据结构】
原文链接如下 [MATLAB深度学习工具箱]学习笔记--体脂估计Body Fat Estimation_bear_miao的博客-CSDN博客介绍本示例展示一个函数拟合神经网络如何根据解剖学测量结果估 ...
【MATLAB深度学习工具箱】学习笔记--体脂估计算例再分析：拟合神经网络fitnet里面的函数】
介绍上一篇 [MATLAB深度学习工具箱]学习笔记--体脂估计算例再分析:拟合神经网络fitnet里面的数据结构]_bear_miao的博客-CSDN博客原文链接如下[MATLAB深度学习工具箱]学 ...
【MATLAB强化学习工具箱】学习笔记--actor网络和critic网络的结果放在哪里？
原算例见 [MATLAB强化学习工具箱]学习笔记--在Simulink环境中训练智能体Create Simulink Environment and Train Agent_bear_miao的博客- ...
MATLAB强化学习-appdesigner使用
MATLAB强化学习-appdesigner使用好像是MATLAB2017a之后都可以使用appdesigner,来代替老的MATLAB的guide. 首先在命令窗口输入:appdesigner打开 ...
第 09 章基于特征匹配的英文印刷字符识别 MATLAB深度学习实战案例
基于特征匹配的英文印刷字符识别 MATLAB深度学习实战话不多讲,直接开撸代码 MainForm函数 function MainForm global bw; global bl; global b ...
台大郭彦甫-Matlab软件学习课堂exercise示例（第二讲）
台大郭彦甫-Matlab软件学习课堂exercise示例 (仅供参考) 第二讲基本操作与矩阵输入 (P6 exercise) >> cos(((1+2+3+4)^3/5)^(1/2))a ...
基于Matlab深度学习Yolov4-tiny的交通标志识别道路标志识别检测
交通标志检测是辅助驾驶.自动驾驶系统中的重要组成部分,针对交通标志检测任务中复杂环境下的小目标检测精度低的问题,提出一种基于YOLOv4-tiny的交通标志检测方法. 基于Matlab深度学习的道路标 ...

MATLAB微分方程学习总结