微分代數方程(DAE)的Matlab 解法
微分代數方程(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=[1 0 0;0 1 0;0 0 0];%品質矩陣
options=odeset('mass',M);%對以DAE問題，mass屬性必須設定
x0=[0.8;0.1;0.1];%初值
[t,x]=ode15s(odefun,[0 20],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]
%狀态變量一階微分初值，題