Provide Jacobian matrix to ode15i for solving DAE system

As introduced by the help page of ode15i http://matlab.izmiran.ru/help/techdoc/ref/ode15i.html, providing Jacobian (df/dy, df/dyp) is important for the efficiency and reliability of ode15i.

For solving DAE system by ode15i, Matlab gives the example of Robertson https://uk.mathworks.com/help/matlab/ref/ode15i.html. In this example, only Jacobian of df/dyp was provided, while the Jacobian of df/dy was passed as anempy matrix [ ].

To practice, I tried to provide both df/dy and df/dyp to ode15i to solve Robertson problem.

The function to calculate these two Jacobian are:

function [dfdy, dfdyp] = Jac_robertsidae(t, y, yp)
 
dfdy = [-0.04, 1e4*y(3), 1e4*y(2); 0.04, -1e4*y(3)-3e7*2*y(2), -1e4*y(2); 1, 1, 1];
dfdyp = [-1, 0, 0; -1, 0, 0; 0, 0, 0];

end

The function for the Robertson DAE system is:

function res = robertsidae(t,y,yp)
   res = [- yp(1) - 0.04*y(1) + 1e4*y(2)*y(3);
   - yp(2) + 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];
end

Then I solve DAE by ode15i, proving the Jacobian for both df/dy and df/dyp:

y0 = [1; 0; 0];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[],yp0,[]);

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], 'Jacobian', @Jac_robertsidae);

tspan = [0 4*logspace(-6,6)];
[t,y] = ode15i(@robertsidae,tspan,y0,yp0, options);

However, a warning as below was given:

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 don't know if it's the problem of how the Jacobian was passed. Many thanks for any help in advance.

Answer 1

[Torsten]

y0 = [1; 0; 0];
yp0 = [0; 0; 0];
options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], 'Jacobian', @Jac_robertsidae);
[y0,yp0] = decic(@robertsidae,0,y0,[],yp0,[],options);

tspan = [0 4*logspace(-6,6)];
[t,y] = ode15i(@robertsidae,tspan,y0,yp0, options);

y(:,2) = 1e4*y(:,2);
semilogx(t,y)
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

function [dfdy,dfdyp] = Jac_robertsidae(t, y, yp) 
  dfdy = [-0.04, 1e4*y(3), 1e4*y(2); 0.04, -1e4*y(3)-3e7*2*y(2), -1e4*y(2); 1, 1, 1];
  dfdyp = [-1, 0, 0; 0, -1, 0; 0, 0, 0];
end

function res = robertsidae(t,y,yp)
  res = [-yp(1) - 0.04*y(1) + 1e4*y(2)*y(3);
         -yp(2) + 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2;
         y(1) + y(2) + y(3) - 1];
end

[Yanxin Liu]

It works, thank you so much, Torsten. I see the problem comes from my way of representing the DAE in the function 'robertsidae', which also leads to the sign flip of my Jacobian.

The ode15i page said the function for the DAE system needs to return the residual value of the equation. My understanding is as long as the terms adds up to 0, the signs can be flipped. And if I don't provide the Jacobian, solving the DAE in flipped sign by ode15i indeed gives the same result. Would you mind explaining why the signs matter when Jacobian are provided? Thank you so much.

(I have a DAE model consists of 36 equations which has been coded up. If I can change the model to a structure suitable for being solved by ode15i by just moving dy/dt to the right hand side (RHS) of ODEs by adding a minus sign, it would be much easier and avoid uneccesary errors when changing all the signs in the RHS of ODEs.)

[Torsten]

As you can see from "my" code, the only thing wrong with yours was

dfdyp = [-1, 0, 0; -1, 0, 0; 0, 0, 0];

It must be

dfdyp = [-1, 0, 0; 0, -1, 0; 0, 0, 0];

[Yanxin Liu]

Ah, I see! Thank you!