DAE Problem for implicit method. Discrete dynamic system.

%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

[Image Analyst]

Is that daeic12.m?

Where is line 63 that says

F = UM

[Cesar García Echeverry]

no, it is an internal function of matlab for solving ode15s

Answer 1

[Walter Roberson]

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.

[Cesar García Echeverry]

Thank you for your response. It is a discrete model, where a need to solve the X balance 100 times (Nest dimension). It means that i have, for the explicit variable, X 100 initial condition vector and for the implicit variable, Y 100 initial conditions vector. What you mean, is that i need Mass matrix of 200x2?

[Walter Roberson]

You can either loop over all of the X Y pairs doing one at a time and using the 2 x 2 mass matrix; OR you can do all the initial conditions at the same time, in which case the mass matrix would have to be 200 x 200 .

For the 200 x 200 mass matrix, you could

tM = repmat({[0 1; 0 0]}, 1, Nest);
M = blkdiag(tM{:});

[Cesar García Echeverry]

Thank you for your approximation. But now i am getting this error: Error using daeic12 (line 78)

This DAE appears to be of index greater than 1. It is for the deffinitions of the mass matrix. For one set [0 1;0 0], i am obtaining a correct answer, but when i increase the dimension, it presents an error.

[Walter Roberson]

Ah... DAE of index 1 implies that only 1 derivative is required to solve the equations algebraically. But you are doing 100 different differential equations simultaneously, so you would be needing 100 different derivatives. You know that they would not interact, but MATLAB does not know; as far as it can tell the index is 100 .

You will probably need to switch to looping instead of running all of them at the same time.

[Cesar García Echeverry]

yeaph, this could be the solution or try with ode15 as a fully implicit method.

Thank you

[Cesar García Echeverry]

Nope, it does not work looping because the solutions of the next stage depend of the before stage:

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

It should be solved at the same time.

[Walter Roberson]

"If wishes were horses, then beggars would ride."

If your calculation for the N'th X, Y depends upon the results of the calculation for a different X, Y pair, then your DAE is going to be of order greater than 1, and ode15s() cannot be used.

Perhaps you can use ode15i()

[Cesar García Echeverry]

If you see the Figure; it is explained my problem. I need to solve the Problem in Figure 1 as is presented in Figure 2 using the set of equations for x and y. Being a implicit system to solve (DAE). The proble is that is a continuous-discrete system.

[Walter Roberson]

    vard(i,1)=beta.*X(i-1)+Y(i+1)-beta.*X(i)-Y(i);

What do your Nest represent? For example does Nest = 100 indicate that you are subdividing an object into 100 intervals, in which case the (i-1) and (i+1) indexing are talking about adjacent intervals? If it does then your DAE is of order greater than 1 and you cannot use ode15s.

Or does Nest = 100 indicate 100 time steps, in which case the reference to Y(i+1) is indicated to implicitly indicate how to find Y(i+1) with vard(i,1) already being known at the time Y(i+1) is to be calculated?

[Cesar García Echeverry]

"Nest" represent that i'm dividing an object into 100 intervarls or N. For this reason i have a DAE dynamic-discrete system. Any advice to solve it? Thank you in advance.

[Walter Roberson]

I do not have any experience with those. Perhaps PDE Toolbox might help?

[Cesar García Echeverry]

Actually, nope. I have a set of differential equations on discrete domain. But thank you for your time

[Walter Roberson]

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]

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