MATLAB Answers

ode15s using jacobian for a stiff problem

2 views (last 30 days)
Serhat Unal
Serhat Unal on 28 Sep 2020
Commented: Alan Stevens on 29 Sep 2020
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!

  0 Comments

Sign in to comment.

Answers (1)

Alan Stevens
Alan Stevens on 28 Sep 2020
Edited: Alan Stevens on 28 Sep 2020
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 Comments

Show 1 older comment
Alan Stevens
Alan Stevens on 29 Sep 2020
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
Serhat Unal on 29 Sep 2020
function dx = dx(t,flag)
syms x y z
c = [77.27;8.375*10^-6;1.161];
switch flag
case ''
dx = [ c(1).*(y+x.*(1-c(2).*x-y));(1./c(1)).*(z-(1+x.*y));c(3).*(x-z) ];
case 'jacobian'
dx = jacobian([c(1).*(y+x.*(1-c(2).*x-y)), (1./c(1)).*(z-(1+x.*y)), c(3).*(x-z)], [x,y,z]);
otherwise
error(['Unknown flag ''' flag '''.']);
end
I have tried so far and have come to here in my function, but I get a error message about that switch that it
is not a scalar, hope one can help to figure out the errors here.
Alan Stevens
Alan Stevens on 29 Sep 2020
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!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by