1. General
  2. Błędy
0
Rozpoczęte

Calculation of Hankel Function

Michael_ 10 lat temu • Ostatnio zmodyfikowane przez Pavel Holoborodko 10 lat temu • 9
Hello.

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

besselj(0,mp('71.62771527592053131293674232738074i'))+1i*bessely(0,mp('71.62771527592053131293674232738074i'))=
-6.103515625e-05 - 0.0003108494982263580776736011003369421i

which is quite different. Now with higher precision

mp.Digits(200) %Higher Precision
besselj(0,mp('71.62771527592053131293674232738074i'))+1i*bessely(0,mp('71.62771527592053131293674232738074i'))=
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)
besselj(0,71.62771527592053131293674232738074i)+1i*bessely(0,71.62771527592053131293674232738074i)=
.000000000000000e+00 - 7.347258660573146e-33i

I assume this is a bug is it?

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

Normal Double Precision:
bessely(0,40*i) = -5.343061323058116e-19 + 1.489477479341990e+16i
Rozpoczęte
Pavel Holoborodko 10 lat temu
Michael_ 10 lat temu
Hello.

Is there any information on this problem? Or is it fixed already, did not check it lately.
Pavel Holoborodko 10 lat temu
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.
Pavel Holoborodko 10 lat temu
Michael,

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.
Michael_ 10 lat temu
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.
Pavel Holoborodko 10 lat temu
Compilation for Windows just has been finished: http://goo.gl/jcoexM

Would appreciate your feedback.
Michael_ 10 lat temu
Hello.
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.
Pavel Holoborodko 10 lat temu
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.

Customer support service by UserEcho