线性常系数微分方程(组)
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=tcos(t),x(1)=1,x(6)=−0.5
关于微分方程的更多了解
- 拉普拉斯变换 拉氏变换求解线性常微分方程
- 偏微分方程
- Simulink的微分方程框图求解
- 更多,欢迎补充
说明
- 以上的程序均在MATLAB 2019a测试通过
- 上面留下的练习题或多或少都有点坑,有的我解决了,有的还没,欢迎讨论
- 主要参考教程:薛定宇《高等应用数学问题的MATLAB求解 (第四版)》 清华大学出版社
- 有的问题我也一直没有解决,欢迎讨论