lsode error: "CORRECTOR CONVERGENCE FAILED REPEATEDLY"

[John B. Thoo]

Hi, everyone.

I'm having trouble with lsode, speficially, I'm getting the following  
error:

 LSODE--  AT T (=R1) AND STEP SIZE H (=R2), THE
       CORRECTOR CONVERGENCE FAILED REPEATEDLY
       OR WITH ABS(H) = HMIN
      In above,  R1 =   .4901607966255E-02   R2 =   .4558107208897-152
error: lsode: repeated convergence failures (t = 0.00490161perhaps  
bad jacobian supplied or wrong choice of integration method or  
tolerances)


So, I have a few questions.  (I searched the Octave help forum for  
"CORRECTOR CONVERGENCE FAILED REPEATEDLY" and that brought up four  
threads; unfortunately, none of them helped me much.  One thread did  
suggest relaxing the tolerances, but, being new to numerics and all  
this, I'm not sure what would be the effect of this on my particular  
equation.)

1) In my code (given below), I have dt = 0.02, so why does lsode  
report "repeated convergence failures (t = 0.00490161...)" when  
0.00490161 < 0.02?

2) Page 310 of the manual lists several lsode_options.  Should I try  
any of these?  Which ones?

3) Do you notice anything wrong with my code or any possible  
streamlining?


%%----------begin my code-----------
% This script solves the pde
%
%    u_t + 0.5*(H[(H[u])^2])_xx + H[u]*u_xx = 0,
%
% where  H[ ]  is the Hilbert transform, using a "pseudospectral"  
method.


clear all;

% set time parameters
t0 = 0;   % initial time
tf = 2;   % final time
M = 100;   % number of time intervals
dt = (tf - t0)/M;   % time mesh
t = t0:dt:tf;   % time interval

% set space parameters
ord = 8;   % 2^ord grid points
x = 0:(2^ord-1);   % }
x = x/2^ord;       % } space interval
x = 2*pi*x;        % }
dx = x(2)-x(1);   % space mesh width


% define initial condition  u0 = u(x,0)
u0 = cos (x);


global lx;
lx = length (x);

% define the equation
function k = make_k (n)   % to multiply  k  by the fft  correctly
  k = 1:n;
  k = heaviside(n/2 + 1 - k).*(-1*(k - 1)) + heaviside(k - n/2).*(n  
- k + 1);
  k = -k;
endfunction

global k;
k = make_k (lx);


function dydx = deriv (y, k)
  dydx = ifft (1i*k.*fft (y).');
endfunction

function d2ydx2 = dderiv (y, k)
  d2ydx2 = ifft (-k.*k.*fft (y).');
endfunction

function out = truncate (in, l)   % to truncate n/2 high wave numbers
  in = fft (in);                  % fft --> truncate ---> ifft
  in(l/4+1:3*l/4) = 0;
  out = ifft (in);
endfunction

function u_t = F (u, t)
global k;
global lx;

  Hu = hilbert (u, lx);
  u_xx = dderiv (u, k);

  Hu2 = Hu.*Hu;
  HHu2 = hilbert (Hu2, lx);
  halfHHu2_xx = 0.5*dderiv (HHu2, k);

  u_t = - halfHHu2_xx.' - Hu.*u_xx.';
  u_t = truncate (u_t, lx);
endfunction


% solve the equation
u = lsode ("F", u0, t).';


% plot the solution
clf

%% static plot
%plot(x, u(:,1), 'r', x, u(:,2), 'g', x, u(:,3), 'b')

%% animation
wait = max (1/(tf - t0 + 1), 0.05);   % pause time between frames

N = 10;   % number of time intervals to plot
time = ones (1, N);
time(2:N+1) = floor (M/N)*(1:N);

for j = 1:N+1
  plot (x, u(:,time(j)));
  title (strcat ('time =', num2str ((tf-t0)*time(j)/M)));
  axis ([0 2*pi -1 1]);
%  axis([0 2*pi 2*min(real(u0)) 2*max(real(u0))]);
  pause(wait);
endfor
%%----------end my code-----------


Thank you very much for any help you can give.

---John.