Supplying Jacobian to lsode

[vHF]

Hi all,

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?

Thanks in advance,

vHF

[vHF]

vHF wrote:

  res(1) = -3.0;

This obviously reads "res(1,1) = -3.0" in my script file.

vHF

[John W. Eaton-3]

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 ({@f, @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
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[vHF]

Dear John,

John W. Eaton-3 wrote:

Try

  x = lsode ({@f, @jac}, init, u);

instead.

Brilliant -- this works.

John W. Eaton wrote:

I think there is an error in your Jacobian.  
...
BTW, this system appears to be unstable, and all the interesting
behavior seems to happen before u = 1.

I should have mentioned that I submitted the smallest example
I could think of, not the real thing. Now I am going back to 60+
equations, so your performance tips are also going to be very useful!

Thanks,

Marek (vHF)