Euler implicite/arrière pour DAE

Question

J'ai une question concernant Euler implicite. Je sais comment calculer la méthode implicite d'Euler mais mon problème est de savoir comment l'utiliser sur DAE (équation algébrique différentielle). J'ai obtenu une solution correcte après avoir appliqué la réduction d'index sur mon DAE d'origine et ainsi j'ai obtenu l'ODE et ensuite j'ai appliqué mon Euler implicite. Cependant, la tâche consistait à déployer Euler implicite sur DAE. Quelqu'un peut-il me donner un indice sur la façon d'améliorer mon code afin qu'il fonctionne aussi pour DAE? Merci beaucoup et s'il vous plaît voir mon code ci-joint.

Voici ma solution pour le problème:

[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 

Answer 1

Les ode15s et ode23t solveurs peuvent résoudre index-1 problèmes de façon linéaire implicite avec une matrice de masse singulière. MATLAB: Solve Differential Algebraic Equations (DAEs)