Vectorizing DAE with strong State Dependance Mass matrix

Hello,

Please help me to vectorize my DAE, I want Vectorized='on'

Thank you

Andrii

a=0.4; b=0.1; c=0.9;
A = @(y) c/y(1);
Ms=@(t,y) [A(y) 0 -a/y(2); 0 b 0; 0 0 0];
Ds=@(t,y) [y(1)*y(3)-y(2)+1;  y(1)-1; a*y(1)+y(2)+y(3)];
Y0=[1, 1, -2];

opts = odeset('MStateDependence','strong','MassSingular','yes','Mass',Ms,Vectorized='off');
%opts = odeset('MStateDependence','strong','MassSingular','yes','Mass',Ms,Vectorized='on');
[t,Y] = ode15s(Ds,[0 10],Y0,opts);

plot(t,Y,"-x")

Answer 1

[Chuguang Pan]

The Vectorized option in odeset means that your DAEs can be writen as , apparently your DAEs can not be expressed as separable form.

[a,b,c] = deal(.4,.1,.9);
Ms = @(t,y) [c./y(1),0,-a./y(2);0,b,0;zeros(1,3)];
Ds = @(t,y) [y(1).*y(3)-y(2)+1;y(1)-1;a.*y(1)+y(2)+y(3)];
Y0=[1;1;-2];
opts = odeset('MStateDependence','strong','MassSingular','yes','Mass',Ms,Vectorized='off');
[t,Y] = ode15s(Ds,[0 10],Y0,opts);

plot(t,Y)

[Andrii]

I was thinking if I vectorize Ms and Ds it will do the job. I there any way to vectorize Ms and Ds? I could try to solve it as ODE, I can do Ms not singular.

[Chuguang Pan]

@Andrii There is an example which use Vectorized option=="on" for solving stiff ODEs, it maybe helpful.

Solve ODEs with Vectorize

[Torsten]

@Andrii

You should supply consistent initial values Y0 for the unknowns. Since your Y0 vector does not satisfy the last algebraic equation in t = 0, you will get almost arbitrary solutions for y(1) - y(3). The reason is the following:

If Y0 does not satisfy the algebraic equation, MATLAB must modify one (or more) of the given initial values in Y0=[1;1;-2]. But which one ?

If MATLAB assumes that y(1) and y(2) are differential variables and y(3) is an algebraic variable, it will start with y(1) = y(2) = 1 and y(3) = -a*1 - 1 = -1.4.

If MATLAB assumes that y(1) and y(3) are differential variables and y(2) is an algebraic variable, it will start with y(1) = 1, y(3) = -2 and y(2) = -a*1 + 2 = 1.6.

If MATLAB assumes that y(2) and y(3) are differential variables and y(1) is an algebraic variable, it will start with y(2) = 1, y(3) = -2 and y(1) = (-1 + 2)/a = 2.5.

In all three cases, you will get different solutions.

Looking at what MATLAB returns for y at time t = 0 shows that the chosen Y0 is even more complicated:

[a,b,c] = deal(.4,.1,.9);
Ms = @(t,y) [c./y(1),0,-a./y(2);0,b,0;zeros(1,3)];
Ds = @(t,y) [y(1).*y(3)-y(2)+1;y(1)-1;a.*y(1)+y(2)+y(3)];
Y0=[1;1;-2];
opts = odeset('MStateDependence','strong','MassSingular','yes','Mass',Ms,Vectorized='off');
[t,Y] = ode15s(Ds,[0 10],Y0,opts);
Y(1,:)

ans = 1×3

1.2692 1.0054 -1.5131

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The vectorized option for a system of only three equations really doesn't matter.

[Andrii]

Thanks for the help, my problem is order of magnitude bigger (19x19) but can reach around (40x40) if it would make sence. Problem has several nonlinear functions (B-H curves) in Mass matrix. I was thinking that Vectorization will improve speed of calculation. My last equation is Kirchhoff's law about sum of currents. I can leave them as is in this case I will get zerous in Mass matrix, or can differentiate them and get ones in mass matrix and zero in right part. Both solution gives me the same result. What I see that there no ways to vectorize my functions - Mass or Mass^-1 in case of non- singular. Which is not obvios why not.

[Torsten]

Which is not obvios why not.

If MATLAB would accept 'Vectorize','on' in case of a mass matrix, you would have to return a matrix which has vectors as elements, thus a 3d- instead of a 2d - object.

If you want to try to vectorize the code, you could:

• Differentiate the algebraic equation
• Don't define a mass matrix, but just multiply f by M^(-1) in the function defining the ODEs

Then, without vectorization, your code would be like below. I'm not able to vectorize it - I had to use a for-loop in function "fun" to deal with matrix inputs for y.

Note that you get different results compared to the code with algebraic equation. The reason again is that your Y0 vector is not consistent.

[a,b,c] = deal(.4,.1,.9);
Y0=[1;1;-2];
options = odeset('Vectorized','on');
[t,Y] = ode15s(@(t,y)fun(t,y,a,b,c),[0 10],Y0,options);
plot(t,Y)

function dydt = fun(t,y,a,b,c)
  dydt = zeros(3,size(y,2));
  for i = 1:size(y,2)
    M = [c/y(1,i),0,-a/y(2,i);0,b,0;a,1,1];
    f = [y(1,i)*y(3,i)-y(2,i)+1;y(1,i)-1;0];
    dydt(:,i) = M\f;
  end
end

[Andrii]

Thanks, appreciate your help. Lets see if vectorization will save me some calculation time.

Regards

Andrii

[Torsten]

Lets see if vectorization will save me some calculation time.

But the code from above is not vectorized - you are still left with the need to vectorize the loop.