Re: fsolve

[Ted Harding]

( Re Message From: Francesco Potorti` )
> 
> I wrote:
>    Hi, is there someone willing to help me to use fsolve?  
> 
> [skipped ... ]
> 
> The reason why it did not converge is a problem qith quad(), which I
> used inside the function to solve.  Indeed, I had to integrate from 0
> to Inf a gaussian with mean 100 and sigma 2, which is sensibly
> different from 0 only around 100.  quad missed that region and told
> that the integral was around 1e-14 (it should have been around 1).  I
> could not either give quad a singularity vector, because that one is
> not accepted with infinite intervals.  So I had to break the
> integration in two infinite intervals, (-Inf, 100) and (100, Inf).
> 
> 1st issue:
>   Shouldn't suggestions like this be given in the info manual?
> 
> 2nd issue:
>   I could not find a way to get the values of quad_options and
>   fsolve_options from inside a script.  Shoulnd't quad_options("abs")
>   suffice?  It does not work.
> 
> 3rd issue:
>   The info manual does not explicitely say if the fourth optional
>   argument to quad is an absolute or relative tolerance.
> 
> Thanks for listening.
> --
>        Francesco Potorti`     | address@hidden (Internet)
>        researcher at          | 39369::pot (DECnet)
>        CNUCE-CNR, Pisa, Italy | +39-50-593203 (voice) 904052 (fax)
> 

I don't have access to the detailed procedure of octave's "quad" function,
which appears to be built-in. However, I looked at the implementation in
an old version of Matlab (which is an m-file). This proceeds by "adaptive
recursive application of Simpson's rule", halving the interval each time
and comparing the new with the previous by relative error
("if abs(Q - Q0) > tol*abs(Q)"). 

With a function like the Gaussian distrbution, integrated over a very wide
range, "quad" could quite easily decide that this criterion was satisfied
after only a few iterations.  By splitting it as you did initially at the
peak of the function, you have ensured that this is less likely. All the
same, for this function I would not have adopted such a crude approach.

However, the real answer to such problems is to focus the numerical method
more precisely in terms of the detailed behaviour of the function. There
is no procedure which is capable of anticipating how a function might
possibly behave at points it has not looked at, unless the method is given
some information about the function or is adapted to a class of functions
which includes the one to which it is being applied. This is too wide an
issue (i.e. is really an issue of the user's appreciation of numerical
analysis) to be covered by your "1st issue: Shouldn't suggestions like
this be given in the info manual?".

However, I agree that in any computational package it is necessary for the
user to know how procedures work, since otherwise it becomes very
difficult for the user to "focus the numerical method more precisely in
terms of the detailed behaviour of the function". This, I agree, is a
documentation issue, and your "3rd issue" is relevant.

Ted.                                    (address@hidden)