Skip site navigation (1)Skip section navigation (2)
Date:      Thu, 6 Sep 2018 04:09:05 +1000 (EST)
From:      Bruce Evans <brde@optusnet.com.au>
To:        Steve Kargl <sgk@troutmask.apl.washington.edu>
Cc:        freebsd-numerics@freebsd.org
Subject:   Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2)
Message-ID:  <20180906034525.A2959@besplex.bde.org>
In-Reply-To: <20180905152104.GA26453@troutmask.apl.washington.edu>
References:  <20180903235724.GA95333@troutmask.apl.washington.edu> <20180905201540.D1142@besplex.bde.org> <20180905152104.GA26453@troutmask.apl.washington.edu>

next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, 5 Sep 2018, Steve Kargl wrote:

> On Wed, Sep 05, 2018 at 10:06:29PM +1000, Bruce Evans wrote:
>> On Mon, 3 Sep 2018, Steve Kargl wrote:
>>
>>> Anyone know where the approximations for j0 (and y0) come from?
>>
>> I think they are ordinary minimax rational approximations for related
>> functions.  As you noticed, the asymptotic expansion doesn't work below
>> about x = 8 (it is off by about 10% for j0(2).  But we want to use the
>> single formula given by the asymptotic expansion for all the subintervals:
>
> I've scoured the literature and web for methods of computing
> Bessel functions.  These functions are important to my real
> work.  I have not found any paper, webpage, documentation, etc.
> that describes what "the related functions" are.

They are just the functions in the asymptotic expansion with errors corrected
as I discussed.

>> ...
>> However, we use the asympotic formula for all cases above x = 2.  It is
>> off by at least 10% unless pzero and qzero are adjusted.  But 10% is not
>> much, and adjusting only pzero and qzero apparently works to compensate
>> for this.  To use the same formula for all cases, the complicated factors
>> cc, ss and sqrt(x) must not be folded into pzero and qzero; instead these
>> factors supply some of the errors which must be compensated for.
>>
>> To find the adjustments, I would try applying the same scale factor to
>> the nominal pzero and qzero.  E.g., for j0, the error factor is
>> invsqrtpi*(u*cc-v*ss)/sqrt(x) / j0(x).  Divide the nominal pzero and qzero
>> by this.  Apply normal minimax methods to the modfided pzero and qzero.
>
> Obviously, I had not thought about scaling P0(x) and Q0(x) with a
> common factor to force agreement between j0(x) and its approximation.

I think I understand why we use the asymptotic formula as a base.  If
we choose a minimax function for each interval, this would need to
have a very high order.  Much the same order as the rational function
part of the asymptotic expansion.  But the latter has a clever grouping
of terms using cos and sin and in the internals of P0 and Q0.  This
makes it faster to calculate and also probably has better accuracy.

Bruce



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180906034525.A2959>