Provide Jacobian matrix to ode15i for solving DAE system

10 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!

Translated by