Calculation of Hankel Function

Michael_ 10 років тому оновлено Pavel Holoborodko 9 років тому 9

I tried to calculate the Hankel Function with mp objects. But as I saw it is not implemented yet. Therefore I used the definition to calculate it via Bessel functions. Following example:

First with double precision:
besselh(0,1,71.62771527592053131293674232738074i) =
0.000000000000000e+00 - 7.347258660573146e-33i

since besselh is not implemented yet, I used besselh(0,1,x)=besselj(0,x)+1i*bessely(0,x) to calculate it.

mp.Digits(34) %Quadruple precision

-6.103515625e-05 - 0.0003108494982263580776736011003369421i

which is quite different. Now with higher precision

mp.Digits(200) %Higher Precision
0 - 7.34725866057318428016340868623722...e-33i

seems better.

mp.Digits(500) %Higher Precisionbesselj(0,mp('71.62771527592053131293674232738074i'))+1i*bessely(0,mp('71.62771527592053131293674232738074i'))=
0 -7.34725866057318428016340868623722892521743368377951878300459469616931...e-33i

Interesting is the following: (using double precision)
.000000000000000e+00 - 7.347258660573146e-33i

I assume this is a bug is it?

Thank you for your answer.
A short further comment:
From Mathematica
N[BesselY[0, 40*I], {1000, 1000}] =
-5.34306132305811469773184318452720415....e-19+1.4894774793419899924224591570721184444946....e16 i
From Matlab
bessely(0,mp('40*i')) = -2.546479089470325 + 1.48947747934199e+16i
bessely(0,mp('40*i')) = 7.114780385429948677221705e-09 + 14894774793419899.92422459i

Normal Double Precision:
bessely(0,40*i) = -5.343061323058116e-19 + 1.489477479341990e+16i

Is there any information on this problem? Or is it fixed already, did not check it lately.
Hello Michael,

There were no updates to this issue recently.
(I have been busy with updating dense matrices engine).

Will resume bessely code update shortly. Sorry for the delay.

I have just released Linux version of toolbox with fixed bessely: http://goo.gl/btsBt
(you use Linux version, right?)

Now it provides accurate result up to the last digit of precision specified (and last digit is correctly rounded).
The main drawback - function is slow for large arguments. Will improve this in next updates.

Please test it in your conditions.
Thank you for the very fast update, but I am using Windows.
I will test it right away, when the release for windows is available.
Compilation for Windows just has been finished: http://goo.gl/jcoexM

Would appreciate your feedback.
Thank you. It gives the correct results now. Performance is indeed 1000x slower with 34 Digits precision compared to normal double, but this is not a big deal for me at the moment.
That is because toolbox switches to multiprecision mode (if necessary) - to compute the bessel accurately.

Now toolbox is using adaptive algorithm for series summation - every operation is analysed for cancellation/rounding errors. 
If it is of unacceptable level - precision is boosted up to compensate the issue. At the end we re-arrange the terms to guarantee the minimum accuracy loss overall [1]. 

It will became much faster in future, when we will add asymptotics for large arguments and (probably) simplify the algorithm.

Hope it is useful for your work.

Btw, our system solver has been updated today with faster decompositions under the hood. You might see some changes in computing time.

[1]. J. Demmel. Y. Hida. Fast and accurate floating point summation with application to computational geometry.