2 views (last 30 days)

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!

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

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.

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!

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

Start Hunting!
## 0 Comments

Sign in to comment.