Provide Jacobian matrix to ode15i for solving DAE system
11 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Yanxin Liu
am 7 Okt. 2022
Kommentiert: Yanxin Liu
am 7 Okt. 2022
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 Kommentare
Akzeptierte Antwort
Torsten
am 7 Okt. 2022
Bearbeitet: Torsten
am 7 Okt. 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Linear Algebra finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!