fsolve() nasty bug ?

[John W. Eaton]

On 21-Apr-2006, Christophe Prud'homme wrote:

| while writing an octave/matlab compatible Crank Nicolson routine,
| I came up with, I believe, an fsolve issue.
| Please see the cranknic.m file attached and corresponding test 
| cranknic_test2.m that check the error for an ode system
| 
| If I write 
| 
| w = fsolve( 'myfun', the_y );
| 
| then fsolve does not converge
| 
| If I write 
| 
| [w info msg] = fsolve( 'myfun', the_y );
| 
| then fsolve _does_ converge. Actually returning just info is sufficient to 
| make fsolve() work.
| 
| I suspect there is a nasty bug somewhere in the implementation of fsolve(), 
| that is disabled when some extra arguments are returned by fsolve().
| 
| any comments ?
| 
| PS : I have tested this with 2.1.73 and 2.9.5

I don't see how this can happen.  In src/DLD-FUNCTIONS/fsolve.cc, the
following code is the only use of the nonlinear equation solving object:

      NLEqn nleqn (x, nleqn_fcn);
      nleqn.set_options (fsolve_opts);

      octave_idx_type info;
      ColumnVector soln = nleqn.solve (info);

      if (fcn_name.length())
        clear_function (fcn_name);
      if (jac_name.length())
        clear_function (jac_name);

      if (! error_state)
        {
          std::string msg = nleqn.error_message ();

          retval(2) = msg;
          retval(1) = static_cast<double> (hybrd_info_to_fsolve_info (info));

          retval(0) = soln;

          if (! nleqn.solution_ok () && nargout < 2)
            error ("fsolve: %s", msg.c_str ());
        }

I don't see how there could possibly be a difference among any of

  x = fsolve (...);
  [x, info] = fsolve (...);
  [x, info, msg] = fsolve (...);

since the code path is exactly the same and INFO is computed and returned
whether or not you ask for it.

Also, in the code you sent, the only calls to fsolve that I see are
 
| for the_t=tt(2:end)
|     if ( exist('OCTAVE_VERSION') )
|       [w, info] = fsolve("myfun",the_y);
|     else
|         w = fsolve('myfun',the_y,options);
|     end
|     u = [u w];
|     the_y = w;
| end

so when running your code with Octave, it seems that only one case is
tested.  Or am I missing something here?

If you think there is a bug, please submit a complete bug report to
the address@hidden list with an example that clearly demonstrates the
problem (it would help if you could eliminate the on-the-fly creation
of a function from another function if that is not central to the bug
report).

Thanks,

jwe