ode15s using jacobian for a stiff problem
显示 更早的评论
Hello everyone!
I am trying to use ode15s to solve a stiff problem using a jacobian implementation, but I am getting some errors from matlab.
My code is the following:
script:
clear all
clc
tspan = [0,100];
options = odeset('jacobian','on','AbsTol',1e-13);
[T,X] = ode15s(@dx,tspan,[1 2 3],options);
function:
function dx = dx(t,x,flag)
dx = zeros(3,1);
c = [77.27;8.375*10^-6;1.161];
dx(1) = c(1).*(x(2)+x(1).*(1-c(2).*x(1)-x(2)));
dx(2) = (1./c(1)).*(x(3)-(1+x(1).*x(2)));
dx(3) = c(3).*(x(1)-x(3));
switch flag
case ''
dx = [ dx(1);dx(2);dx(3) ];
case 'jacobian'
dx = [ dx(1);dx(2);dx(3) ];
otherwise
error(['Unknown flag ''' flag '''.']);
end
Appreciate if somebody can help me find where the error is. Thanks!
回答（1 个）
Alan Stevens
2020-9-28
编辑：Alan Stevens
2020-9-28
The Jacobian option doesn't have an 'on' value; you have to supply a pointer to a Jacobian function (look at the documentation for odeset). However, ode15s seems to work perfectly well without this:(in fact I'm not sure what it does for you here)
tspan = [0,100];
options = odeset('AbsTol',1e-13);
[T,X] = ode15s(@dxfn,tspan,[1 2 3],options);
subplot(3,1,1)
plot(T,X(:,1)),grid
xlabel('t'),ylabel('x1')
subplot(3,1,2)
plot(T,X(:,2)),grid
xlabel('t'),ylabel('x2')
subplot(3,1,3)
plot(T,X(:,3)),grid
xlabel('t'),ylabel('x3')
function dx = dxfn(~,x)
dx = zeros(3,1);
c = [77.27;8.375*10^-6;1.161];
dx(1) = c(1).*(x(2)+x(1).*(1-c(2).*x(1)-x(2)));
dx(2) = (1./c(1)).*(x(3)-(1+x(1).*x(2)));
dx(3) = c(3).*(x(1)-x(3));
dx = [ dx(1);dx(2);dx(3) ];
end
This results in
4 个评论
Serhat Unal
2020-9-28
Alan Stevens
2020-9-29
编辑：Alan Stevens
2020-9-29
Unfortunately, I've never heard of, nor know anything about the Oregonator fitting. However, according to the odeset documentation you need to supply a pointer to the actual Jacobian; there is no 'on'/'off' switch.
If you have Simulink there is an Oregonator model in the File Exchange.
Serhat Unal
2020-9-29
Alan Stevens
2020-9-29
You can select the jacobian option in switch by calling the routine as follows:
[T,X] = ode15s(@dxfn,tspan,[1 2 3],[],'jacobian');
Note that I changed the name of the function to dxfn as you can't have the same function name as the variable
i.e. you cant have
function dx = dx(...etc)
so changed it to
function dx = dxfn(...etc)
However, there are two problems.
First, the jacobian can't be calculated from doubles, so you could define a separate function, say
function Jcb = Jcbfn()
c = [77.27;8.375*10^-6;1.161];
syms x y z
Jcb = jacobian([c(1).*(y+x.*(1-c(2).*x-y)), (1./c(1)).*(z-(1+x.*y)), c(3).*(x-z)], [x,y,z]);
end
then call it with
.....
case 'jacobian'
dx = Jcbfn();
.....
BUT, this isn't very useful as function dxfn must return a 3x1 column vector, but the Jacobian is a 3x3 matrix.
I confess, I don't know why the Jacobian is of use in this situation anyway!
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
2020-9-28
2020-9-29
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
