Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Jacobian in stiff ode

Subject: Jacobian in stiff ode

From: Mateusz Gos

Date: 23 Mar, 2009 10:40:04

Message: 1 of 5

Hi,
I am running a code (pasted below) and it seems I cannot specify an analytical rectangular Jacobian matrix to support ode23s, it has to be square. But as far as I am concerned the Jacobian matrix may very well be rectangular, so why this limitation?

function model
coeff = [5 3/5];
ics = [4;7;-3];
options = odeset('RelTol',1e-6,'AbsTol',1e-6,'Jacobian',@(t,Y) jac(t,Y,coeff));
[t,Y] = ode23s(@(t, Y) integ(t,Y,coeff), [0 3000], ics, options);

display([t(end) Y(end,:)])
end
%--------------------------------------------------------------------------
function J = jac(t,Y,coeff)

a = coeff(1); b = coeff(2);
%{
rectangular Jacobian
J(1,:) = [ cos(a*Y(1))*Y(1), 0, cos(a*Y(1))*a, 0, 0];
J(2,:) = [ 0, -1/log(b*Y(2))^2*log(2)/b, 0, -1/log(b*Y(2))^2*log(2)/Y(2), 0];
J(3,:) = [ exp(a+b*Y(3)), Y(3)*exp(a+b*Y(3)), 0, 0, b*exp(a+b*Y(3))];
%}

J(1,:) = [ cos(a*Y(1))*a, 0, 0];
J(2,:) = [ 0, -1/log(b*Y(2))^2*log(2)/Y(2), 0];
J(3,:) = [ 0, 0, b*exp(a+b*Y(3))];

end
%--------------------------------------------------------------------------
function dY = integ(t,Y,coeff)

a = coeff(1); b = coeff(2);

dY(1) = sin(a*Y(1));
dY(2) = 1/log2(b*Y(2));
dY(3) = exp(a+b*Y(3));

dY = dY';
end


Cheers,
Mateusz

Subject: Jacobian in stiff ode

From: Torsten Hennig

Date: 23 Mar, 2009 11:00:03

Message: 2 of 5

> Hi,
> I am running a code (pasted below) and it seems I
> cannot specify an analytical rectangular Jacobian
> matrix to support ode23s, it has to be square. But as
> far as I am concerned the Jacobian matrix may very
> well be rectangular, so why this limitation?
>
> function model
> coeff = [5 3/5];
> ics = [4;7;-3];
> options =
> odeset('RelTol',1e-6,'AbsTol',1e-6,'Jacobian',@(t,Y)
> jac(t,Y,coeff));
> [t,Y] = ode23s(@(t, Y) integ(t,Y,coeff), [0 3000],
> ics, options);
>
> display([t(end) Y(end,:)])
> end
> %-----------------------------------------------------
> ---------------------
> function J = jac(t,Y,coeff)
>
> a = coeff(1); b = coeff(2);
> %{
> rectangular Jacobian
> J(1,:) = [ cos(a*Y(1))*Y(1),
> 0, cos(a*Y(1))*a,
> cos(a*Y(1))*a, 0,
> 0, 0];
> J(2,:) = [ 0,
> -1/log(b*Y(2))^2*log(2)/b, 0,
> , -1/log(b*Y(2))^2*log(2)/Y(2),
> 0];
> J(3,:) = [ exp(a+b*Y(3)),
> Y(3)*exp(a+b*Y(3)),
> a+b*Y(3)), 0,
> 0,
> 0, b*exp(a+b*Y(3))];
> %}
>
> J(1,:) = [ cos(a*Y(1))*a,
> 0, 0];
> J(2,:) = [ 0,
> -1/log(b*Y(2))^2*log(2)/Y(2),
> 0];
> J(3,:) = [ 0,
> 0,
> 0, b*exp(a+b*Y(3))];
>
> end
> %-----------------------------------------------------
> ---------------------
> function dY = integ(t,Y,coeff)
>
> a = coeff(1); b = coeff(2);
>
> dY(1) = sin(a*Y(1));
> dY(2) = 1/log2(b*Y(2));
> dY(3) = exp(a+b*Y(3));
>
> dY = dY';
> end
>
>
> Cheers,
> Mateusz

You need n differential equations to determine n
unknowns ; so why do you think square Jacobians
are a limitation ?

Best wishes
Torsten.

Subject: Jacobian in stiff ode

From: Mateusz Gos

Date: 23 Mar, 2009 12:06:01

Message: 3 of 5

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <20015.1237806033697.JavaMail.jakarta@nitrogen.mathforum.org>...
>
> You need n differential equations to determine n
> unknowns ; so why do you think square Jacobians
> are a limitation ?
>
> Best wishes
> Torsten.

Well, I admit this may be the problem of my not understanding what Jacobian matrix is, but when constructing Jacobian matrix we do not differentiate only over the variables we want to determine using the given set of differential equations (then indeed for a 'well-conditioned' problem we'll end up with a square matrix). In my case it would be Y(1:3), 3x3 matrix. I believe we differentiate over EVERY variable which varies with time. Actually in this case 'a' and 'b' do not differ, but they could just as well depend on Ys and some other constants - then we would have to differentiate over them and as a result and up with a rectangular matrix.

Unless there is a major error in my thinking… (which may be the case)

Regards,
Mateusz

Subject: Jacobian in stiff ode

From: Torsten Hennig

Date: 23 Mar, 2009 12:54:39

Message: 4 of 5

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in message
> <20015.1237806033697.JavaMail.jakarta@nitrogen.mathfor
> um.org>...
> >
> > You need n differential equations to determine n
> > unknowns ; so why do you think square Jacobians
> > are a limitation ?
> >
> > Best wishes
> > Torsten.
>
> Well, I admit this may be the problem of my not
> understanding what Jacobian matrix is, but when
> constructing Jacobian matrix we do not differentiate
> only over the variables we want to determine using
> the given set of differential equations (then indeed
> for a 'well-conditioned' problem we'll end up with a
> square matrix). In my case it would be Y(1:3), 3x3
> matrix. I believe we differentiate over EVERY
> variable which varies with time. Actually in this
> case 'a' and 'b' do not differ, but they could just
> as well depend on Ys and some other constants - then
> we would have to differentiate over them and as a
> result and up with a rectangular matrix.
>
> Unless there is a major error in my thinking…
> (which may be the case)
>
> Regards,
> Mateusz

For implicit solution methods in differential equations,
the Jacobian matrix is needed for Newton's method to
determine the solution variables for time t + deltat
given the solution at time t. For this purpose, the
derivatives of the right-hand sides of the system
of differential equations with respect to the unknowns
y1,...,yn are needed.
If your a and b were dependent on the y's,
you had to take this into account and explicitely
calculate da/dyi and db/dyi to supply the correct
derivatives.
(E.g.if your equation reads dy/dt = a*y and a = 5*y^2,
then you had to supply a + y*da/dy = 5*y^2 + 10*y^2 =
15*y^2 as Jacobian).

Best wishes
Torsten.

Subject: Jacobian in stiff ode

From: Mateusz Gos

Date: 23 Mar, 2009 13:50:18

Message: 5 of 5

Thanks a lot, that's the answer I needed.

Mateusz

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us