微分代数方程

微分代数方程(DAE)的Matlab解法

所谓微分代数方程，是指在微分方程中，某些变量满足某些代数方程的约束。假设微分方程的更一般形式可以写成

前面所介绍的ODEs数值解法主要针对能够转换为一阶常微分方程组的类型，故DAE就无法使用前面介绍的常微分方程解法直接求解，必须借助DAE的特殊解法。

其实对于我们使用Matlab求解DAE时，却没有太大的改变只需增加一个Mass

参数即可。描述f(t,x)的方法和普通微分方程完全一致。

注意：ode15i没法设置Mass属性，换句话说除了ode15i外其他ode计算器都可以求解DAEs问题

1.当M(t,y)非奇异的时候，我们可以将微分方程等效的转换为

y'=inv(M)*f(t,y)，此时就是一个普通的ODE(当然我们可以将它当成DAEs处理)，对任意一个给定的初值条件都有唯一的解

2.当m(t,y)奇异时，我们叫它为DAEs(微分代数方程)，DAEs问题只有在同时提供状态变量初值y0和状态变量一阶导数初值py0，且满足M(t0,y0)*yp0=f(t0,y0)时才有唯一解，假如不满足上面的方程，DAEs解算器会将提供的y0和py0作为猜测初始值，并重新计算与提供初值最近的封闭初值

3.质量矩阵可是一个常数矩阵(稀疏矩阵)，也可以是一个自定义函数的输出。但是ode23s只能求解Mass是常数的DAEs

4.对于Mass奇异的DAEs问题，特别是设置MassSingular为yes时，只能ode15s 和ode23t解算器

5.对于DAE我们还有几个参数需要介绍

a.Mass：质量矩阵，不说了，这个是DAE的关键，后面看例子就明白了

b.MStateDependence：质量矩阵M(t,y)是否是y的函数，可以选择

none|{weak}|strong，none表示M与y无关，weak和strong都表示与y相关c.MvPattern：注意这个必须是稀疏矩阵，S(i,j)=1表示M(t,y)的第i行中任意元素都与第j个状态变量yi有关，否则为0

d.MassSingular：设置Mass矩阵是否奇异，当设置为yes时，只能使用ode15s 和ode23t

e.InitialSlope：状态变量的一阶导数初值yp0，和y0具有相同的size，当使用ode15s和ode23t时，该属性默认为0

下面我们以实例说明，看下面的例子，求解该方程的数值解

【解】

真是万幸，选取状态变量和求状态变量的一阶导数等，微分方程转换工作，题目已经帮我们完成。

可是细心的网友会发现，最后一个方程不是微分方程而是一个代数方程(这就是为什叫DAE的原因)，其实我们可以将它视为对三个状态变量的约束。

(1)用矩阵形式表示出该DAEs

(2)编写Matlab代码

复制内容到剪贴板

代码:

odefun=@(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)*x(2);

x(1)+x(2)+x(3)-1];%微分方程组

M=[100;010;000];%质量矩阵

options=odeset('mass',M);%对以DAE问题，mass属性必须设置

x0=[0.8;0.1;0.1];%初值

[t,x]=ode15s(odefun,[020],x0,options);%这里好像不能使用ode45

figure('numbertitle','off','name','DAE demo—by Matlabsky')

plot(t,x)

legend('x1(t)','x2(t)','x3(t)')

在介绍隐式ODEs时，我们说了ode15i也可以求解DAEs，原理就是将约束看成一个隐式，好，下面我们看看，验证下所说的：

复制内容到剪贴板

代码:

odefun=@(t,x,dx)[dx(1)+0.2*x(1)-x(2)*x(3)-0.3*x(1)*x(2);

dx(2)-2*x(1)*x(2)+5*x(2)*x(3)+2*x(2)^2;

x(1)+x(2)+x(3)-1];

%状态变量初值，题目中给出

x0=[0.8;0.1;0.1];

%该初值全部给定，按理说应该全部为1，但是记住方程中有一个约束，故其实只有两个独立初值

x0_fix=[1;0;1];%随意一个改为0都可以，比如[0;1;1]或者[1;1;0]

%状态变量一阶微分初值，题目没有给出，可以任意写

dx0=[1;1;1];

%该初值一个都没有给出，故全部为0

dx0_fix=[0;0;0];

%时间变量的初值

t0=0;

%计算相容初值

[x0,dx0]=decic(odefun,0,x0,x0_fix,dx0,dx0_fix);

[x0,dx0]%相容初始条件

res=ode15i(odefun,[0,20],x0,dx0);

figure('numbertitle','off','name','DAE demo—by Matlabsky')

plot(res.x,res.y)

其实，有些简单DAEs是可以转换为ODEs的，由于代数方程只是状态变量的一

个约束，假如约束充分的话，我们就可以使用消参的思想，将约束反代回ODEs 中，这样就将那个约束方程消去了，最后只剩下那些微分方程了。

对于本例，我们只要通过第三个方程表示出

再将x3反代回前两个微分方程中则有

这个就是一个正常的ODEs了，好，下面我们就重新使用Matlab求解下看：

复制内容到剪贴板