ODE15i DAE Warning:Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (0.000000e+00) at time t.

I'm solving a DAE using ode15i. I use the same syntaxs as the example of Robertson in Matlab help

https://it.mathworks.com/help/matlab/ref/ode15i.html?s_tid=doc_ta

The script is

h = 2/2; % h = L/n, L = 2, n = 2

y0 = [0;-2/h^2;1;2/h^2;1/h;-1/h;1/h;0];
yp0 = [
   -2.0000
   -0.0545
    2.0000
    0.0545
         0
         0
         0
         0];


%The initial conditions must be "consistent", meaning 
%   that f(T0,Y0,YP0)=0
% With y0 and yp0 given, Node_2(0,y0,yp0) = [0 ... 0] so initial conditions
% are consistent

%opt = odeset('RelTol',10^-3,'AbsTol',10^-6,'Jacobian',@Jac_Node_2); %Optionnal

tspan = [0 4];
[t,y] = ode15i(@Node_2,tspan,y0,yp0);%,opt);

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (0.000000e+00) at time t.

And the odefunction is :

function res = Node_2(t,y,yp)
h = 2/2;
m = 25000*0.12^2/9.81;

res = [
(y(3) - y(1))/h - y(5);
(y(6)-y(5))/h - m*yp(2);
y(1);
yp(1) - y(2);
(y(8)-y(3))/h - y(6);
(y(7)-y(6))/h - m*yp(4);
yp(3) - y(4);
y(8)
];

end 

The problem is that i get the following error :

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (0.000000e+00) at time t.

I tried to change the relative and absolute error tolerance 'RelTol' and 'AbsTol' but nothing change for every possible values, i tried to provide the Jacobian matrix but again it gives the same error.

The initial conditions are consistent because Node_2(0,y0,yp0) = [0 ... 0]

I would like to understand why there is this error so that i can find a solution, or i you spot the error in the code i would appreciate.

Answer 1

[Torsten]

It can be deduced from your equations that the only possible solution for your DAE system is y(i) = 0 for i=1,...,8.

Thus you cannot prescribe values different from 0 as you do in your vectors y0 and yp0.

If you have questions about the derivation below, feel free to ask.

   
(1) (y(3) - y(1))/h - y(5) = 0 
(2) (y(6)-y(5))/h - m*yp(2) = 0
(3) y(1) = 0
(4) yp(1) - y(2) = 0
(5) (y(8)-y(3))/h - y(6) = 0
(6) (y(7)-y(6))/h - m*yp(4) = 0
(7) yp(3) - y(4) = 0
(8) y(8) = 0


//1 (3) -> y1 = 0      
//2 (8) -> y8 = 0      
//3 (1) and //1 -> y3 = h*y5   
//4 (5) and //2 -> y3 = -h*y6  
//5 (//3 and //4) -> y5 = -y6
//6 (4) and //1 -> y2 = 0  
//7 (2) and //6  -> y5 = y6    
//8 (//7 and //5) -> y5 = y6 = 0
//9 (//1 and //8) and (1) -> y3 = 0 
//10 (7) and //9 -> y4 = 0
//11 (//8 and //10) and (6) -> y7 = 0