Failure when using ode15i : not meeting integration tolerances

Hello !

Thanks to your answers, I have been able to solve a DAE using ode15i. The code runs but only until a specific time, and after it it stops:

My code solves the following DAE :

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

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

and f1 and f2 are functions calculating the force corresponding to x1,x and dx/dt, dx1/dt.

f1 and f2 look like this :

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

if dx - dx1 <= 0 %COMPRESSION
    force = ... (a polynom P(x - x1))
 
else %DETENTE
    force = ... (another polynom P(x - x1))
        
end

I then use it in my DAE :

function res = f(t,Y,dY)

% Y = (x; x1; dx) and dY = (dx; dx1; d²x)

course1 = l0(1) - Y(2);
coursetot = l0(1) + l0(2) - Y(1);
course2 = coursetot - course1;

forcemi79 = f1(course1, dY(2))
forcemi20 = f_MI20ode15i(course2,Y(3)-dY(2));

res = [dY(1) - Y(3);
    - forcemi20 + forcemi79;
    dY(3) - 1/m_materiel*(forcemi20)];

And then I use ode15i :

sol = ode15i(@f,tspan,y0,yp0)

But when dx - dx1 goes from negative to positive, the code stops, with this error :

Warning: Failure at t=5.111363e-01.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.815921e-15) at time t. 
> In ode15i at 406
  In wagon79face20 at 31 

I don't understand how I can change f1 and f2 in order for it to run anyway ?

Thank you for your help !

[Torsten]

A function like your f1 will immediately throw an error:

a = 25;
b = 31;
force = f1(a,b)

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

Invalid use of operator.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

  if dx - dx1 <= 0 %COMPRESSION
    force = x-x1; 
  else %DETENTE
    force = (x-x1).^2;      
  end
end

The input arguments to a function cannot be arithmetic expressions of some variables.

It must read:

function force = f1(x ,x1, dx, dx1)
  if dx-dx1 <= 0 %COMPRESSION
    force = x-x1; 
  else %DETENTE
    force = (x-x1).^2;      
  end
end

[Zoé Cord'homme]

HI ! thanks for your answer. I have changed f1 and f2 but I still have the same error ! I did not have this issue when I used ode45 on a similar problem, so maybe it is due to my use of ode15i ?

[Walter Roberson]

@Zoé Cord'homme

What is your current code?

[Zoé Cord'homme]

I have my main :

clear
close all

l0 = [0.09 0.2];
             
V0 = -1;

y0 = [l0(1)+l0(2);l0(1);V0];

yp0 = [V0;V0;0];


dt = 0.01;
tf = 0.355;
tspan = 0:dt:tf;

sol = ode15i(@f,tspan,y0,yp0);

[y,dy]=deval(sol,tspan);

... some plots

And then my DAE :

function res = f(t,Y,dY)

l0 = [0.09 0.2];
m_materiel = 265959;

forcemi79 = f1(l0(1),Y(2),0, dY(2),0);

forcemi20 = f2(l0(2),Y(1),Y(2),Y(3),dY(2));

res = [dY(1) - Y(3);
    - forcemi20 + forcemi79;
    dY(3) - 1/m_materiel*(forcemi20)];

And my forces :

function forcemi79=f1(l0,y2,y1,v2,v1)

course = l0 - (y2 - y1);
vrel = v2 - v1;

if vrel <= 0 %COMPRESSION
    if course<=0
        forcemi79 = ...a polynom P(course)
    elseif course<=...
        forcemi79 = ...another polynom P(course)
    else
        forcemi79=0;
    end
    
else
    if course<=...
        forcemi79=...another polynom P(course)
    elseif course <=...
        forcemi79 = ...another polynom P(course)
    else
        forcemi79=0;
    end
end

And when I run it, I get the following error :

Warning: Failure at t=7.649144e-02.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.717522e-16) at time t. 
> In ode15i at 406
  In main at 16 

Thanks !

[Walter Roberson]

        forcemi79 = ...a polynom P(course)
    elseif course<=...

The three dots at in the forcemi79 assignment act as continuation markers. The rest of that line is ignored, and the next line is treated as if it was part of the same line, so you effectively have a line

        forcemi79 = elseif course<= forcemi79 else forcemi79=0

which is not valid syntax.

Reminder: if you are hoping that someone can debug your code, they need your actual code to test with.

[Zoé Cord'homme]

Sorry I will link the exact code :

clear
close all

l0 = [0.09 0.2];
             
V0 = -1;

y0 = [l0(1)+l0(2);l0(1);V0];

yp0 = [V0;V0;0];


dt = 0.01;
tf = 1.5;
tspan = 0:dt:tf;

sol = ode15i(@f,tspan,y0,yp0);%,options);

[y,dy]=deval(sol,tspan);

where f is :

%  Y = (x; x1; dx), dY = (dx; dx1; d²x)

function res = f(t,Y,dY)

l0 = [0.09 0.2];
m_materiel = 265959;

forcemi79 = f1(l0(1),Y(2),0, dY(2),0);

forcemi20 = f2(l0(2),Y(1),Y(2),Y(3),dY(2),1.3);

res = [dY(1) - Y(3);
    - forcemi20 + forcemi79;
    dY(3) - 1/m_materiel*(forcemi20)];

And f1 :

function forcemi79=f1(l0,y2,y1,v2,v1)

course = l0 - (y2 - y1);
vrel = v2 - v1;
if vrel <= 0 %COMPRESSION
    if course<=0
        forcemi79 = 0;
    elseif course<=0.01327
        forcemi79 = (6000/0.01327)*course;
    elseif course <=0.05073
        forcemi79 = ((310000-6000)/(0.05073-0.01327))*course + 310000 - ((310000-6000)/(0.05073-0.01327))*0.05073;
    elseif course <= 0.05074
        forcemi79 = ((750000-310000)/(0.05074-0.05073))*course + 750000 - ((750000-310000)/(0.05074-0.05073))*0.05074;
    elseif course <= 0.06073
        forcemi79 = ((823000-750000)/(0.06073-0.05074))*course + 823000 - ((823000-750000)/(0.06073-0.05074))*0.06073;
    elseif course <= 0.07823
        forcemi79 = ((748000-823000)/(0.07823-0.06073))*course + 748000 - ((748000-823000)/(0.07823-0.06073))*0.07823;
    elseif course <= 0.13353
        forcemi79 = ((999000-748000)/(0.13353-0.07823))*course + 999000 - ((999000-748000)/(0.13353-0.07823))*0.13353;
    elseif course <= 0.14073
        forcemi79 = ((930000-999000)/(0.14073-0.13353))*course + 930000 - ((930000-999000)/(0.14073-0.13353))*0.14073;
    else
        forcemi79=0;
    end
    
else
    if course<=0.01327
        forcemi79=0;
    elseif course <=0.05073
        forcemi79 = (106660/(0.05073-0.01327))*course + 106660 - (106660/(0.05073-0.01327))*0.05073;
    elseif course <= 0.05074
        forcemi79 = ((16000-106660)/(0.05074-0.05073))*course + 16000 - ((16000-106660)/(0.05074-0.05073))*0.05074;
    elseif course <= 0.11073
        forcemi79 = ((82000-16000)/(0.11073-0.05074))*course + 82000 -((82000-16000)/(0.11073-0.05074))*0.11073;
    elseif course1 <= 0.14073
        forcemi79 = ((290000-82000)/(0.14073-0.11073))*course + 290000 - ((290000-82000)/(0.14073-0.11073))*0.14073;
    else
        forcemi79=0;
    end
end

And f2 :

function forcemi20 = f_MI20ode15i(l0,y2,y1,v2,v1,alpha)

course = l0 - (y2-y1);
vrel = v2-v1;
mu = 0.15461;
sigma = 0.050566;
z = (course-mu)/sigma;

%Coefficients:
p1 = 364.9;
p2 = 4483.8;
p3 = 22046;
p4 = 54527;
p5 = 69312;
p6 = 43622;
p7 = 32992;
p8 = 75168;
p9 = 1.157e+05;
p10 = 1.4066e+05;
p11 = 1.7428e+05;

if vrel<= 0 %COMPRESSION
    if course<=0
        forcemi20 = 0;
    elseif course<=0.01
        forcemi20 = (150000/0.01)*course;
    elseif course<=0.2
        forcemi20 = p1*z^10 + p2*z^9 + p3*z^8 + p4*z^7 +p5*z^6 + p6*z^5 +p7*z^4 + p8*z^3 +p9*z.^2 + p10*z +p11+100000;
    else
        forcemi20 = 0; 
    end

else %DETENTE
    if course<=0
        forcemi20 =0;
    elseif course<=0.01
        forcemi20 =(50000/0.01)*course;
    elseif course<=0.199
        forcemi20 = p1*z^10 + p2*z^9 + p3*z^8 + p4*z^7 +p5*z^6 + p6*z^5 +p7*z^4 + p8*z^3 +p9*z.^2 + p10*z +p11;
    elseif course<=0.2
        forcemi20 = ((667950-543010)/(0.2-0.199))*course+667950-((667950-543010)/(0.2-0.199))*0.2;
    else
        forcemi20 =0; 
    end
end 

Thank you for your time !

Answer 1

[Walter Roberson]

I am a lot less familiar with ode15i, but I think your equations need to have continuous derivatives, at least one step. Which would mean that your different force functions need to be arranged to have at least first (and possibly second) derivatives continuous with each other across vrel changing between compression and relaxation . (The same problem occurs in friction system in effectively switching the direction of opposition to motion.)

If you were using ode45() or ode23s() or similar, then I would be certain that this is a problem; however, I see that ode15i uses a different calculation system that I know less about.

When you change state non-smoothly then you risk exactly the kind of error you see.

The work-around is to use an event function to detect each time that you would change state, and have it terminate the call there; then take the output of that, make any necessary adjustments to the boundary conditions (such as reflecting a bounce or an input of energy), and then call the ode* function again to continue from that time.

The rule for ode45() and ode23s() and similar is not exactly that you cannot use if statements: however if the if statements do not lead to smooth first and second derivatives at that location then you need to arrange so that the discontinuity of state does not occur within one call

[Zoé Cord'homme]

OK thank you very much for this thourough answer ! That's what I thought I might have to do, but I hoped I could avoid it !

Thank you :)