Mass Matrix SWITCH Problem in Matlab

I'm having a problem with my ODE in Matlab. The code to start the ODE solver is the following:

refine = 4; 
  options = odeset('Refine', refine,...
      'RelTol', 1e-9,'absTol',1e-9*ones(1,length(y_sym)),'Mass',Mss_bend,'MassSingular','no','MstateDependence','strong')

[t_bend_i,y_bend_i] = ode45(@(t,y) ...
      bending_DGL(t,y,time_r,dt,ug,vg,Mss_bend,Fss_bend,Sys,flag),time_r_i,IC_bend_i,options);

And the corresponding function file containing the ODE is:

function dy = bending_DGL(t,y,time_r,dt,ug,vg,Mss_bend,Fss_bend,Sys,flag)
r = Sys.r;
roh = Sys.roh;
A = Sys.A;
I = Sys.I;
ne = Sys.ne;
nnode = Sys.nnode;
ndof = Sys.ndof;
Le = Sys.Le;
E = Sys.E;
g = Sys.g;
EI = E*I;
EA = E*A;
%Interpolate Earthquake Excitation
i = floor(t/dt) + 1;
if i < length(time_r)
     u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
     v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));
else
      u = 0;
      v = 0;
end
%Load State Space Representation of Problem
switch flag
    case ''
        dy = Fss_bend;
    case 'mass'
        dy = Mss_bend;
end

Now I'm trying to solve the ODE where Mss_bend is the mass matrix which is symbolic, because it is statedependent so Mss_bend(y,t) and Fss_bend(y,t) is the right side of the Equation which is state dependent too. So now I wan't to solve the Problem Mss_bend(y,t)*dy = Fss_bend where dy is the derivative of the vector y with respect to time.

When I start the calculation it gives me the

Error: Switch expression must be a scalar or string constant.

Error in bending_DGL (line 31)
switch flag

Error in @(t,y)bending_DGL(t,y,time_r,dt,ug,vg,Mss_bend,Fss_bend,Sys,flag)

Error in odearguments (line 87)

What did i do wrong? Is there a problem because the matrices are symbolic?

Matlab tell's me in the function file that y is unused:

function dy = bending_DGL(t, y,time_r,dt,ug,vg,Mss_bend,Fss_bend,Sys,flag

How do I solve this Problem. y is as symbolics in the Mss_bend and Fss_bend where y(1), y(2) etc is declared as symbols.

Thank you for your answers!

[Marius]

I found the solution! I had to load the symbolic y vector

y_sym = [y(1), y(2), y(3) ...]

into my function file, end then make

Fss_bend = subs(Fss_bend,y_sym,y)
Fss_bend = subs(Fss_bend,{sym('u'), sym('v')}, {u, v}

dy = double(Fss_bend);

Answer 1

[Walter Roberson]

We would need to see how you initialized your flag variable to figure out why you are getting the error for switch.

You will need to change your code anyhow. Assuming that your Fss_bend and Mss_bend are expressions rather than functions:

%Load State Space Representation of Problem
switch flag
    case ''
        dy = double( subs(Fss_bend, {sym('t', sym('y')}, {t, y}) );
    case 'mass'
        dy = double( subs(Mss_bend, {sym('t', sym('y')}, {t, y}) );
end

However, it seems likely that your Mss_bend and Fss_bend also use some of the variables that you compute, such as u and v. We need to know whether Mss_bend and Fss_bend are functions or expressions, and if they are functions we need to know the order they expect their arguments in. If they are expressions then we need to know which symbolic variables need to be substituted.

[Walter Roberson]

In order for your line

 [t_bend_i,y_bend_i] = ode45(@(t,y) ...
      bending_DGL(t,y,time_r,dt,ug,vg,Mss_bend,Fss_bend,Sys,flag),time_r_i,IC_bend_i,options);

to execute at all, you need to have defined the variables time_r, dt, ug, vg, Mss_bend, Fss_bend, Sys, flag, time_r_i, and IC_bend_i .

I need to see how you initialized the value for flag before you reach that line.

But what I really suspect is that you are handling the mass matrix completely wrong.

MassFun = @(t,y) double(subs(Mss_bend, sym('y'), y));

refine = 4; 
options = odeset('Refine', refine,...
      'RelTol', 1e-9, 'absTol',1e-9*ones(1,length(y_sym)), 'Mass', MasFun, 'MassSingular', 'no', 'MstateDependence', 'strong');

[t_bend_i,y_bend_i] = ode45(@(t,y) bending_DGL(t, y, time_r, dt, ug, vg, Fss_bend, Sys), time_r_i, IC_bend_i ,options);

together with

function dy = bending_DGL(t, y, time_r, dt, ug, vg, Fss_bend, Sys)
r = Sys.r;
roh = Sys.roh;
A = Sys.A;
I = Sys.I;
ne = Sys.ne;
nnode = Sys.nnode;
ndof = Sys.ndof;
Le = Sys.Le;
E = Sys.E;
g = Sys.g;
EI = E*I;
EA = E*A;
%Interpolate Earthquake Excitation
i = floor(t/dt) + 1;
if i < length(time_r)
     u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
     v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));
else
      u = 0;
      v = 0;
end

dy = double( subs(Fss_bend, {sym('u'), sym('v'), sym('y')}, {u, v, y}) );

[Marius]

Thank you! I changed it as you said. But now I'm getting the error:

Error using odearguments (line 92)
@(T,Y)BENDING_DGL(T,Y,TIME_R,DT,UG,VG,FSS_BEND,SYS) returns a vector of length 256, but the
length of initial conditions vector is 16. The vector returned by
@(T,Y)BENDING_DGL(T,Y,TIME_R,DT,UG,VG,FSS_BEND,SYS) and the initial conditions vector must have
the same number of elements.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Why is it building an outputvector of 256 = 16^2 when it should be only 16?

[Torsten]

And Fss_bend really is a symbolic column vector of length 16 ?

Best wishes

Torsten.

[Marius]

Yes it's a symbolic 16x1 vector which looks like:

y(9)
y(10)
y(11)
y(12)
y(13)
y(14)
y(15)
y(16)
cos(y(1))*sin(y(2))*u+y(3)...
sin(y(1))*sin(y(2))*u+y(4)...
cos(y(1))*sin(y(2))*u+y(5)...
sin(y(1))*sin(y(2))*u+y(6)...
cos(y(1))*cos(y(2))*u+y(4)...
cos(y(1))*sin(y(2))*u+y(6)...
cos(y(1))*cos(y(2))*u+y(8)...
cos(y(1))*sin(y(2))*u+y(7)...

Where 'y(1)'...'y(16)' were defined as

%initialize statespace vector

j = 1;

while j<=2*nqi

yname = strcat('y(',num2str(j),')');

y(j) = sym(yname);

j = j+1;

end

[Torsten]

Can't be true.

If it were the case, the command

dy = double( subs(Fss_bend, {sym('u'), sym('v'), sym('y')}, {u, v, y}) );

would generate a 16x1 vector of doubles, not a 256x1 vector of doubles.

Check the sizes of Fss_bend (before dy is generated) and dy.

Best wishes

Torsten.

[Marius]

Hmm I tested it and you are right. If I just run my dy function it returns a 256x1 vector whith 16 different results each in 16 following rows. Like:

dy = 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,...

so it repeats each result 16 times?

[Torsten]

Maybe

dy = double( subs(Fss_bend, {sym('u'), sym('v'), sym('y')}, {u, v, y'}) );

?

Best wishes

Torsten.

Answer 2

[Marius]

Thank you for your answer!

1)We would need to see how you initialized your flag variable to figure out why you are getting the error for switch.

What do you mean with this? All I have done with the flag is already in the code. What did I miss? Sorry I have never worked with flags and the Mass representation of ODE's before.

2)You will need to change your code anyhow. Assuming that your Fss_bend and Mss_bend are expressions rather than functions.

Fss_bend is a 16x1 symbolic vector which contains 'u','v', and 'y(1)' - 'y(16)' as symbolic variables. Mss_bend is a strongly nonlinear 16x16 symbolic matrix which contains 'y(1)' - 'y(16)'. as symbolic variables.

What do I need to change?

Answer 3

[Marius]

Is there another possibility to do this? Is it true that now it is substituting the variables in each step of the ODE15s solver? This takes very long. Isn't there a way that doesn't use subs?

[Walter Roberson]

matlabFunction() the symbolic expressions into function handles outside of the objective function, and pass the function handles in to the objective.

[Marius]

Thank you for your answer! I'm trying to implement matlabFunction. Before I had for the Massfunction for ODE 15s:

MassFun = @(t,y) sparse(double(subs(Mss, y_sym, y)))

where y_sym is vector with the symbolic variables:

   y_sym = ['y(1)';'y(2)'...;'y(24)']

Now I changed it to:

MassFun = matlabFunction(Mss,'vars',y_sym)

When I start the script now it gives me the error message:

Error using sym/matlabFunction>checkVars (line 159)
Variable names must be valid MATLAB variable names.

Error in sym/matlabFunction (line 104)
vars = checkVars(funvars,opts);

What did I do wrong?

[Walter Roberson]

You could try

y_sym = [sym('y[1]'), sym('y[2]'), ..., sym('y[24]')]

I might be able to test this later today.

[Marius]

Mhh it didn't work.

[Walter Roberson]

Sorry, I had an install glitch and cannot test the symbolic toolbox at the moment.

[Marius]

I solved it! This works with:

MassFun = matlabFunction(Mss,'vars',{'t','y'})

But when I try to do the same with Fss, which is a vector with functions of y(1)...y(24), u and v:

FssFun = matlabFunction(Fss,'vars',{'t','y','u','v'})

It gives me the error:

Error using symengine>makeFhandle (line 109)
Error: ()-indexing must appear last in an index expression.

Error in symengine (line 60)

Error in sym/matlabFunction (line 125)
    g = symengine('makeFhandle',varnames,body);

Error in Processing (line 261)
FssFun = matlabFunction(Fss,'vars',{'t','y','u','v'})

[Walter Roberson]

Do not proceed by way of

y_sym = ['y(1)';'y(2)'...;'y(24)']

Instead

y = sym('y', [1, 24]);

which will show up as y1, y2, y3, etc, in your function. Then

matlabFunction(fss, 'vars', {'t', y, 'u', 'v'})

Notice this is y rather than 'y' : when you give a vector of names within {} then all the names are put into a single input argument that is indexed as needed.

Unfortunately the variable name that will be written into the anonymous function will be like 'in2', but since the arguments are positional rather than by name, it will still work.

[Marius]

Thanks for your answer. I tried it but it gives me the same error:

Warning: Function 'y1' is not verified to be a valid MATLAB function.

Warning: Function 'y1' is not verified to be a valid MATLAB function.

Error using symengine>makeFhandle (line 109)
Error: ()-indexing must appear last in an index expression.

Error in symengine (line 60)

Error in sym/matlabFunction (line 125)
    g = symengine('makeFhandle',varnames,body);

Error in Processing (line 263)
   FssFun = matlabFunction(Fss, 'vars', {'t', y, 'u', 'v'})

[Walter Roberson]

I need to see your Fss, at least part of it.

[Marius]

I found the Error. Now it computes MassFun and FssFun:

MassFun = 
    @(t,in2)reshape([1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0...

FssFun = 
    @(t,in2,u,v)reshape([1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0...

But how do I have to start in my ODE file?

I start the ODE solver with:

refine = 1; 
    options = odeset('Refine', refine,...
        'RelTol', 1e-5,'absTol',1e-5*ones(1,length(y_sym)),'Jacobian',@(t,y)JacobiFun(t,y,time_r,dt,ug,vg,J,y_sym,Sys),'JConstant','no','Mass',MassFun,'MassSingular','no','MstateDependence','strong','OutputFcn',@odetpbar,'Events',@events_overturn,'InitialStep',1e-5);

      [t_up,y_up,te,ye,ie] = ode15s(@(t,y) ...
          rocking_DGL(t,y,time_r,dt,ug,vg,Fss,Sys,y_sym),time_r,IC_rocking,options);

And my rocking_DGL file looks like:

function dy = rocking_DGL(t,y,time_r,dt,ug,vg,Fss,Sys,y_sym)

r = Sys.r;
roh = Sys.roh;
A = Sys.A;
I = Sys.I;
ne = Sys.ne;
nnode = Sys.nnode;
ndof = nnode*3;
Le = Sys.Le;
E = Sys.E;
g = Sys.g;
EI = E*I;
EA = E*A;

%Interpolate Earthquake Excitation

i = floor(t/dt) + 1;
if i < length(time_r)
     u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
     v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));

else
      u = 0;
      v = 0;

end

%Load State Space Representation of Problem

dy = FssFun(t,y,u,v);

end

But then it gives me the error:

Undefined function 'FssFun' for input arguments of type 'double'.
Error in rocking_DGL (line 38)
dy = FssFun(t,y,u,v);

[Walter Roberson]

In

function dy = rocking_DGL(t,y,time_r,dt,ug,vg,Fss,Sys,y_sym)

change the Fss to FssFun

[Marius]

Now i got the next error:

Error using odearguments (line 90)
@(T,Y)ROCKING_DGL(T,Y,TIME_R,DT,UG,VG,FSSFUN,SYS) must return a column vector.
Error in ode15s (line 148)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in Processing (line 270)
    [t_up,y_up,te,ye,ie] = ode15s(@(t,y) ...

[Walter Roberson]

What is your FssFun reshape()'ing to? It should be reshaping to a column vector.

If you are using matlabFunction(something) to generate FssFun then use matlabFunction(something(:))

Answer 4

[Marius]

Now it gives me the error:

Error using odearguments (line 92)
  @(T,Y)ROCKING_DGL(T,Y,TIME_R,DT,UG,VG,FSSFUN,SYS) returns a vector of length 576, but the length of initial conditions
  vector is 24. The vector returned by @(T,Y)ROCKING_DGL(T,Y,TIME_R,DT,UG,VG,FSSFUN,SYS) and the initial conditions
  vector must have the same number of elements.
  Error in ode15s (line 148)
  [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
  Error in Processing (line 270)
      [t_up,y_up,te,ye,ie] = ode15s(@(t,y) ...

[Marius]

I had to change it to:

   FssFun = matlabFunction(Fss(:,1), 'vars', {'t', y, 'u', 'v'})