Attached you will find the 3 scripts
o a test script cranknic_test2
o the _good script works at all times when executing cranknic_test2
o the _bad script fails starting at N=2^5 when executing cranknic_test2
the only difference between _good and _bad is the fact that in _good I store
one extra output arg from fsolve()
ie [w info] = fsolve()
instead of just w = fsolve
This is highly suspicious since there is no reason why it should affect the
outcome of fsolve.
Note that the result give by _good when _bad fails is OK
here is the log
octave-2.9.5:1> cranknic_test2
get w and info
get w and info OK
get only w
get only w OK
N = 4
get w and info
get w and info OK
get only w
get only w OK
N = 8
get w and info
get w and info OK
get only w
get only w OK
N = 16
get w and info
get w and info OK
get only w
get only w OK
N = 32
get w and info
get w and info OK
get only w
error: fsolve: iteration is not making good progress
error: evaluating assignment expression near line 24, column 7
error: evaluating for command near line 22, column 1
error: called from `cranknic_bad' in file `/tmp/cranknic_bad.m'
error: evaluating for command near line 8, column 1
error: near line 26 of file `/tmp/cranknic_test2.m'
this is 100% reproducible on different octave versions(2.1.x and 2.9.y) on
different Debian/GNU/Linux distros(stable and unstable)
Best regards
C.
------------------------------------------------------------------------
function [t,u]=cranknic_bad(odefun,tspan,y,Nh,varargin)
%CRANKNIC Solve differential equations using the
% Crank-Nicolson method.
% [T,Y]=CRANKNIC(ODEFUN,TSPAN,Y0,NH) with TSPAN=[T0,TF]
% integrates the system of differential equations
% y'=f(t,y) from time T0 to TF with initial condition
% Y0 using the Crank-Nicolson method on an equispaced
% grid of NH intervals.Function ODEFUN(T,Y) must return
% a column vector corresponding to f(t,y). Each row in
% the solution array Y corresponds to a time returned
% in the column vector T.
% [T,Y] = CRANKNIC(ODEFUN,TSPAN,Y0,NH,P1,P2,...) passes
% the additional parameters P1,P2,... to the function
% ODEFUN as ODEFUN(T,Y,P1,P2...).
tt=linspace(tspan(1),tspan(2),Nh+1);
u(:,1)=y;
global glob_h glob_t glob_y glob_odefun;
glob_h=(tspan(2)-tspan(1))/Nh;
glob_y=y;
glob_odefun=odefun;
for glob_t=tt(2:end)
% get only X from fsolve
w = fsolve('cranknicfun',glob_y);
u = [u w];
glob_y = w;
end
t=tt;
clear glob_h glob_t glob_y glob_odefun;
end
function z=cranknicfun(w)
global glob_h glob_t glob_y glob_odefun;
z=w-glob_y-0.5*glob_h*(feval(glob_odefun,glob_t,w)+feval(glob_odefun,glob_t,glob_y));
end
------------------------------------------------------------------------
1;
y0=[1; 0];
tspan=[0 10];
f=inline('[-y(2); y(1)]', 't', 'y');
y='[cos(t); sin(t)]';
N=2;
for k=1:5
% call good
disp('get w and info')
[tt,u]=cranknic_good(f,tspan,y0,N);
t=tt(end);
e_good(k)=norm(u(:,end)-eval(y));
disp('get w and info OK')
% call bad
disp('get only w')
[tt,u]=cranknic_bad(f,tspan,y0,N);
t=tt(end);
e_bad(k)=norm(u(:,end)-eval(y));
disp('get only w OK')
N=2*N
end
------------------------------------------------------------------------
function [t,u]=cranknic_good(odefun,tspan,y,Nh,varargin)
%CRANKNIC Solve differential equations using the
% Crank-Nicolson method.
% [T,Y]=CRANKNIC(ODEFUN,TSPAN,Y0,NH) with TSPAN=[T0,TF]
% integrates the system of differential equations
% y'=f(t,y) from time T0 to TF with initial condition
% Y0 using the Crank-Nicolson method on an equispaced
% grid of NH intervals.Function ODEFUN(T,Y) must return
% a column vector corresponding to f(t,y). Each row in
% the solution array Y corresponds to a time returned
% in the column vector T.
% [T,Y] = CRANKNIC(ODEFUN,TSPAN,Y0,NH,P1,P2,...) passes
% the additional parameters P1,P2,... to the function
% ODEFUN as ODEFUN(T,Y,P1,P2...).
tt=linspace(tspan(1),tspan(2),Nh+1);
u(:,1)=y;
global glob_h glob_t glob_y glob_odefun;
glob_h=(tspan(2)-tspan(1))/Nh;
glob_y=y;
glob_odefun=odefun;
for glob_t=tt(2:end)
% get only X INFO from fsolve (getting MSG also works)
[w info] = fsolve('cranknicfun',glob_y);
u = [u w];
glob_y = w;
end
t=tt;
clear glob_h glob_t glob_y glob_odefun;
end
function z=cranknicfun(w)
global glob_h glob_t glob_y glob_odefun;
z=w-glob_y-0.5*glob_h*(feval(glob_odefun,glob_t,w)+feval(glob_odefun,glob_t,glob_y));
end