ode15s using jacobian for a stiff problem

Question

Serhat Unal 2020-9-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/601288-ode15s-using-jacobian-for-a-stiff-problem

评论： Alan Stevens 2020-9-29

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2020-9-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/601288-ode15s-using-jacobian-for-a-stiff-problem#answer_501766

编辑：Alan Stevens 2020-9-28

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

Serhat Unal 2020-9-29

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 2020-9-29

在 MATLAB Online 中打开

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!

请先登录，再进行评论。

ode15s using jacobian for a stiff problem

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

ode15s using jacobian for a stiff problem

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论