Provide Jacobian matrix to ode15i for solving DAE system
4 views (last 30 days)
Show older comments
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.
0 Comments
Accepted Answer
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
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!