ode15s using jacobian for a stiff problem

3 vues (au cours des 30 derniers jours)
Serhat Unal
Serhat Unal le 28 Sep 2020
Commenté : Alan Stevens le 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:
clear all
tspan = [0,100];
options = odeset('jacobian','on','AbsTol',1e-13);
[T,X] = ode15s(@dx,tspan,[1 2 3],options);
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) ];
error(['Unknown flag ''' flag '''.']);
Appreciate if somebody can help me find where the error is. Thanks!

Réponses (1)

Alan Stevens
Alan Stevens le 28 Sep 2020
Modifié(e) : Alan Stevens le 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);
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) ];
This results in
  4 commentaires
Alan Stevens
Alan Stevens le 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]);
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!

Connectez-vous pour commenter.

Community Treasure Hunt

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

Start Hunting!

Translated by