Неявный / Обратный Эйлер для DAE

Question

У меня вопрос по поводу неявного Эйлера. Я знаю, как вычислить неявный метод Эйлера, но моя проблема заключается в том, как использовать его в DAE (дифференциальное алгебраическое уравнение). Я получил правильное решение после того, как применил снижение индекса к своему исходному DAE, и, таким образом, я получил ODE, а затем применил свой неявный Эйлер. Однако задача заключалась в развертывании неявного Эйлера на DAE. Может кто-нибудь дать мне подсказку о том, как улучшить мой код, чтобы он работал и для 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

Answer 1

ode15s а также ode23t решатели могут решать линейно-неявные задачи с индексом 1 с матрицей сингулярных масс. MATLAB: решение дифференциальных алгебраических уравнений (DAE)