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.