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!

Answer 1

[Alan Stevens]

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

[Serhat Unal]

Okey, the picture looks fine, but in my case I want to

use the Belousov-Zhabotinsky model which is my dx matrix and I want to fit this with

something called Oregonator fitting. In this case I have to implement the jacobian which is included in

the Oregonator fittin, but I am not

sure if I was on the right track which is shown in my code when I write:

switch flag

case ''

dx = [ dx(1);dx(2);dx(3) ];

case 'jacobian'

dx = [ dx(1);dx(2);dx(3) ];

otherwise

error(['Unknown flag ''' flag '''.']);

Is there maybe something wrong in this code, and also as you mentioned what should be written in option instead. Hope for help.

[Alan Stevens]

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]

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]

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!