RE: fsolve not finding correct answer

[Thompson, David C]

It's not clear to me the relationship between the first
equation for p and the second. If the second should
in fact read
   f_B(y,\sigma)=(1/(pi.*sigma)).*besselk(0,(abs(y)./sigma));
then what is the value of p for which you hope to find
a corresponding k? Since besselk(0,y) is positive for all
positive y, its integral will always be positive and there
is no k that will solve p = 0 except k = 0. In fact, will
the integral of f_B even converge since besselk(0,y) diverges
as y approaches 0 (and you integrate from -k to k)?

    David

-----Original Message-----
From: address@hidden on behalf of Aivo Jürgenson
Sent: Sun 2/11/2007 12:17 PM
To: address@hidden
Subject: fsolve not finding correct answer

Hello,

I'm trying to use fsolve for finding unknown $k$ in a somewhat complex
equation

        $ p = \int_{-k}^{k} f_B(y,\sigma) dy $

where f_B is the probability density function of "Bessel's"
distribution and equals to

        p=(1/(pi.*sigma)).*besselk(0,(abs(y)./sigma));

where besselk is the standard Octave function.

The problem is that the fsolve function doesn't do a very good job and
sometimes gives "iteration is not making good progress" and
therefore completely wrong answer as well.

The appropriate code is given in the atree_bessel_find_range.m file in the
attachment. I'm using the 2.1.73-13 Debian package of Octave.

Anybody has some hints about how to persuade fsolve giving the correct
answer or perhaps somebody could suggest a different function for solving
the equation in the first place?

Aivo Jürgenson