Re: bug in fsolve() affecting convergence depending on the returned args

[Christophe Prud'homme]

[ Vendredi 21 Avril 2006 21:15 ]
| Christophe Prud'homme wrote:
| > 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
|
| If ypu type "which fsolve" where does it tell you fsolve is coming from?
| Is it by chance the octave-forge version?
nope

octave2.9:11> which fsolve
fsolve is the dynamically-linked function from the file
/usr/lib/octave/2.9.5/oct/i486-pc-linux-gnu/fsolve.oct


octave-2.1.73:1> which fsolve
fsolve is the dynamically-linked function from the file
/usr/lib/octave/2.1.73/oct/i486-pc-linux-gnu/fsolve.oct

can you reproduce the behavior with your octave installation ?

cu
C.
-- 
Christophe Prud'homme
EPFL SB IACS CMCS
MA B2 534 (Bâtiment MA)
Station 8
CH-1015 Lausanne
Tel: +41 (0)21 693 25 47
Fax: +41 (0)21 693 43 03