Supplying Jacobian to lsode

[John W. Eaton]

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