DAE Problem for implicit method. Discrete dynamic system.

10 次查看（过去 30 天）
Cesar García Echeverry
编辑： Cesar García Echeverry 2021-5-11
%Parameters
Nest=100;
alfa = 1.2;
beta=0.75;
nu=alfa-1;
h = 0.9;
f = 0;
%Initial values
X=zeros(Nest,1);
Y=zeros(Nest,1);
X(1)=f;
Y(Nest)=h;
%Initial conditions to the DAE Method
init1=[X;Y];
tspan=[0 30];
rpar=[alfa,beta,nu,f,h,Nest];
M=[0 1 ;0 0];
options=odeset('Refine',5,'MassSingular','yes','Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,var]=ode15s(@diffX,tspan,init1,options,rpar);
function [vard]=diffX(t,var,rpar)
alfa=rpar(1);
beta=rpar(2);
nu=rpar(3);
f=rpar(4);
h=rpar(5);
Nest=rpar(6);
% var=reshape(var,[Nest,2]);
X=var(1:Nest);
Y=var(Nest+1:end);
for i=2:Nest-1
vard(i,1)=beta.*X(i-1)+Y(i+1)-beta.*X(i)-Y(i);
vard(i,2)=(alfa.*X(i))./(1+nu.*X(i-1));%Y(i+1)+nu.*X(i-1).*Y(i)-alfa.*X(i);
end
vard=reshape(vard,[2*(Nest-1),1]);
vard=[0;vard;Y(Nest)];
% vard=vard';
end

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2021-5-8
X=zeros(Nest,1);
Y=zeros(Nest,1);
init1=[X;Y];
So your initial conditions are 2*Nest x 1 and your function must return 2*Nest x 1.
For a mass matrix M then M*y_prime = f(x, y) and f must return 2*Nest x 1 and y_prime must be that size. In order to have M*y_prime be the same size as f(x, y) then M must be a square matrix with rows (and columns) the same size as the number of rows in y_prime. Which is 2*Nest
But does M have that size? No. M is 2x2 and so can only be used if Nest is 1.
  14 个评论
显示 12更早的评论
Walter Roberson
Walter Roberson 2021-5-10
PDE are often dealt with by dividing an area into discrete domains.
One of your recent questions mentions distillation columns. Physically those are continuous columns -- but subdividing into discrete parts to model behaviour through the column would be a common approach. It is also the most common approach to PDE.
Cesar García Echeverry
Cesar García Echeverry 2021-5-11
No, actually i solved. I used ode15s and ode15i. Defining M or the Jacobian matrix as {[],[eye(np-2) zeros(np-2); zeros(np-2) zeros(np-2)]} and the algebraic systemas as:
for i=2:n+1
res(i-1)=difXcol(i-1)-(sum(Am([1:np]+(i-1)*np).*Ycol)...
+beta*sum(An([1:np]+(i-1)*np).*Xcol(i))...
-beta*Xcol(i)-Ycol(i));
res(i-1+np-2)=Ycol(i)+nu.*Xcol(i).*Ycol(i)-alfa.*Xcol(i);
end
Plus the boundary conditions

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by