Provide Jacobian matrix to ode15i for solving DAE system

5 次查看（过去 30 天）
Yanxin Liu
Yanxin Liu 2022-10-7
评论： Yanxin Liu 2022-10-7
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.

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-10-7
编辑：Torsten 2022-10-7
Ran in:
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 个评论
显示 1更早的评论
Torsten
Torsten 2022-10-7
编辑：Torsten 2022-10-7
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];

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by