Provide Jacobian matrix to ode15i for solving DAE system

8 views (last 30 days)
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.

Accepted Answer

Torsten
Torsten on 7 Oct 2022
Edited: Torsten on 7 Oct 2022
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
  3 Comments

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!