help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## Supplying Jacobian to lsode

 From: John W. Eaton Subject: Supplying Jacobian to lsode Date: Sun, 6 Sep 2009 02:22:28 -0400

```On  5-Sep-2009, vHF wrote:

| I am relatively new to Octave and my problem is most likely rather trivial,
| so I apologise in advance.
| I am trying to create an .m file that would have a definition of a system of
| ODEs and call the lsode solver with the Jacobian supplied. What I have now
| is:
|
| ***begin script***
| u = linspace(0, 72, 72*60);
| init = [0; 10];
|
| function res = f(x, u)
|   res(1) = 2.0*x(2).*x(2) - 3.0*x(1);
|   res(2) = -2.0*x(2).*x(2) + 6.0*x(1);
| endfunction;
|
| function res = jac(x, u)
|   res(1) = -3.0;
| res(1,2) = 2.0*x(2);
|   res(2,1) = 6.0;
| res(2,2) = -2.0*x(2);
| endfunction;
|
| x = lsode(["f"; "jac"], init, u);
| ***end script***
|
| However, this fails with the following output:
|
| ***begin console output***
| error: `x' undefined near line 5 column 16
| error: evaluating binary operator `*' near line 5, column 15
| error: evaluating binary operator `.*' near line 5, column 20
| error: evaluating binary operator `-' near line 5, column 27
| error: evaluating assignment expression near line 5, column 10
| error: called from `f'
| error: evaluating assignment expression near line 16, column 40
| error: called from `__lsode_fcn__U__'
| error: lsode: evaluation of user-supplied function failed
| error: lsode: inconsistent sizes for state and derivative vectors
| error: evaluating assignment expression near line 16, column 3
| error: near line 16 of file `/home/marek/Desktop/octave/dimer.m'
| ***end console output***
|
|
| What am I doing wrong?

Try

x = lsode (address@hidden, @jac}, init, u);

instead.

I think there is an error in your Jacobian.  I think it should be

-3   4*x(2)
6  -4*x(2)

It would probably be slightly more efficient to write

function res = f (x, u)
res = zeros (2, 1);
tmp = 2*x(2)^2;
res(1) = tmp - 3.0*x(1);
res(2) = -tmp + 6.0*x(1);
endfunction;

function res = jac (x, u)
res = zeros (2, 2);
tmp = 4.0*x(2);
res(1,1) = -3.0;
res(1,2) = tmp;
res(2,1) = 6.0;
res(2,2) = -tmp;
endfunction;

or perhaps

function res = f (x, u)
tmp = 2*x(2)^2;
res = [tmp-3*x(1); -tmp+6*x(1)];
endfunction

function res = jac (x, u)
tmp = 4*x(2);
res = [-3, tmp; 6, -tmp];
endfunction

because your functions are repeatedly resizing the res vectors.

BTW, this system appears to be unstable, and all the interesting
behavior seems to happen before u = 1.

jwe

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]