# ode15s using jacobian for a stiff problem

6 vues (au cours des 30 derniers jours)
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:
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 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponses (1)

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);
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 commentairesAfficher 2 commentaires plus anciensMasquer 2 commentaires plus anciens
Serhat Unal le 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 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]);
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!

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by