How can I solve a second degree DAE ?

Hi !

For a project, I am currently needing to solve a second degree (or of index 2, I am not too familiar with those) DAE.

The equations are :

d²x/dt² = 1/m * (- f2(x-x1, dx/dt - dx1/dt))

0 = f2(x-x1, dx/dt - dx1/dt) - f1(x1,dx1/dt)

Where m is a scalar and f1 and f2 are functions looking like this :

function force=f1(x-x1, dx/dt - dx1/dt)

%Coefficients

p1 = 1.4347e+07;

p2 = 1.9757e+06;

p3 = 3.1841e+05;

if vrel<=0 %COMPRESSION

force = p1*(x-x1)^2 + p2*(x-x1) + p3;

else %DETENTE

force = p1*(x-x1)^2 + p2*(x-x1) + p3 - 260000;

end

So far, I have coded this :

m_materiel= ...;

M=[1 0 0 0

0 1 0 0

0 0 0 0

0 0 0 1]; % mass matrix

V0 =...; %initial condition

y0=[0 0 V0 V0];

dt=0.01;

tf=2;

tspan=0:dt:tf;

options = odeset('Mass',M,'Vectorized','on');

[t,Y]=ode15s(@(t,Y) f(t,Y,m_materiel),tspan,y0,options);

%-----------------------------------------

function out=f(t,Y,m_materiel)

out =[Y(3,:);

Y(4,:);

f_MI7984(Y(1,:),Y(3,:),data1,data2)-fQS_continue_MI20(Y(2,:)-Y(1,:),Y(4,:)-Y(3,:));

1/m*(fQS_continue_MI20(Y(2,:)-Y(1,:),Y(4,:)-Y(3,:)))];

I have used a similar technique than with classical second order differential equations,

Here "out" is dY where Y = [x1; x; dx1/dt; dx/dt]

When I run the code, I get multiple errors :

Error using vertcat

Dimensions of matrices being concatenated are not consistent.

Error in interp1q (line 31)

[~, j] = sort([x;xxi]);

Error in f (line 2)

out =[Y(3,:);

Error in @(t,Y)f(t,Y,m_materiel)

Error in odenumjac (line 143)

Fdel = feval(F,Fargs_expanded{:});

Error in daeic12 (line 37)

[DfDy,Joptions.fac,nF] = odenumjac(fun, {t0,y,args{:}}, f, Joptions);

Error in ode15s (line 310)

[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in attelageeqVSmur (line 24)

[t,Y]=ode15s(@(t,Y) f(t,Y,m_materiel),tspan,y0,options);

Can you help me out ? Do you know how to solve this kind of equation ? Thank you :)

Answer 1

[Torsten]

Write your equations as

dx/dt = x2
dx2/dt = 1/m * (- f2(x-x1, x2 - dx1/dt))
0 = f2(x-x1, x2 - dx1/dt) - f1(x1,dx1/dt)

Are you able to solve the third equation for dx1/dt ?

If yes, insert the expression for dx1/dt in the call to f2 in equation 2 and use ode15s, if no, use ode15i.

[Zoé Cord'homme]

Thank you for this quick answer ! I will have to use ode15i then, but I don't really understand how the equations will look ?

Thank you

[Torsten]

Read the documentation of ode15i.

In principle, the equations will be

res(1) = dy(1) - y(2);
res(2) = dy(2) + 1/m * f2(y(1)-y(3),y(2)-dy(3))
res(3) = f2(y(1)-y(3),y(2)-dy(3)) - f1(y(3),dy(3))

where

y = [x x2 x1]

But better pass y and dy to f1 and f2 instead of the derived values above and form the necessary expressions in the functions themselves. E.g. function definitions as

function force = f1(x-x1, dx/dt - dx1/dt)

are not allowed.

Thus better use

function force = f1(y,dy)
  diff1 = y(1) - y(3);
  diff2 = dy(1) - dy(3);
  ...

[Zoé Cord'homme]

Thank you very much ! Have a great day :)

[Zoé Cord'homme]

Hello again ! the program runs but the solution makes no physical sense at all...

I was wondering if maybe I should be including dx1/dt in y ? and write it as y=[x dx x1 dx1]

Because I feel like I am lacking an equation ?

I have read everything I found about ode15i but it still is quite mysterious to me !

[Torsten]

Please include the equations you try to solve and the code you use.

[Zoé Cord'homme]

found the issue thx ! my initial conditions were inconsistent !