Effacer les filtres
Effacer les filtres

ode15s using jacobian for a stiff problem

6 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:
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!

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);
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 commentaires
Serhat Unal
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
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.

Community Treasure Hunt

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

Start Hunting!

Translated by