DAE 的隐式/后向欧拉

Posted 2023-04-13

【中文标题】DAE 的隐式/后向欧拉

【英文标题】：Implicit/Backward Euler for DAE

【发布时间】：2017-10-18 17:55:01

【问题描述】：

我有一个关于隐式欧拉的问题。我知道如何计算隐式欧拉方法，但我的问题是如何在 DAE（微分代数方程）上使用它。在我对原始 DAE 应用索引缩减后，我获得了正确的解决方案，因此我获得了 ODE，然后我应用了我的隐式欧拉。但是，任务是在 DAE 上部署隐式 Euler。谁能给我一个关于如何改进我的代码以便它也适用于 DAE 的提示？非常感谢，请参阅我的附加代码。

这是我的解决方案：

[t,y]=beul('system','dsystem',[-1,1,-1],0,1,100);
plot(t,y);

function yp=system(t,y)
yp(2)=y(1);  % equations
yp(3)=y(2);
yp(1)=exp(-t);  % after applying index reduction we obtain this
end

function y=dsystem(t,x)
y(1,1)=-1;
y(1,2)=0;
y(1,3)=0;
y(2,1)=0;
y(2,2)=-1;
y(2,3)=0;
y(3,1)=0;
y(3,2)=0;
y(3,3)=-1;


end

function [t,y]=beul(f,df,y0,t0,tf,n)
h=(tf-t0)/n;
t=linspace(t0,tf,n+1);
y=zeros(n+1,length(y0));
y(1,:)=y0;
for i=1:n
k0=y(i,:)';
k1=k0-inv(eye(length(y0))-h*feval(df,t(i),k0))*(k0-h*feval(f,t(i),k0)'-y(i,:)');
while (norm(k1-k0)>0.0001) % Newton evaluation
k0=k1;
k1=k0-inv(eye(length(y0))-h*feval(df,t(i),k0))*(k0-h*feval(f,t(i),k0)'-y(i,:)');
end
y(i+1,:)=k1;
end
end

【问题讨论】：

求解器应该如何自动工作？您可以手动求解隐式欧拉方程，即x3[k+1]=-u3[k+1], x2[k+1]=(x3[k+2]-x3[k])/h, x1[k+1]=(x2[k+1]-x2[k])/h，还是从输入方程自动导出？

【参考方案1】：

ode15s 和 ode23t 求解器可以求解具有奇异质量矩阵的 index-1 线性隐式问题。 MATLAB: Solve Differential Algebraic Equations (DAEs)