Message Boards Message Boards

GROUPS:

HankelH1 real part for large values of n differs from Python result

Posted 11 months ago
2520 Views
|
9 Replies
|
5 Total Likes
|

Dear all,

when evaluating "HankelH1[n,z]" in Mathematica the resulting real part significantly differs from the result I obtain when I use "scipy.special.hankel1(n,z)" in Python for large values of n. For example, "HankelH1[5,3]" and "scipy.special.hankel1(5,3)" are in sufficient agreement, however for "HankelH1[16,3]" and "scipy.special.hankel1(16,3)" the result significantly differs. Does anybody have an idea why that is the case?

Thank you very much.

Kind regards

POSTED BY: D D
9 Replies

You say: "the result significantly differs" .Post a value: "scipy.special.hankel1(16,3)" ?

POSTED BY: Mariusz Iwaniuk
Posted 11 months ago

Hello,

thank you for your reply. Sorry that I did not add in the first post. Here you see the values. Please focus on the real part. In the first line is the result from Mathematica, while in the second line in each paragraph the result obtained from Python. Note that in Python "j" corresponds to "I" in Mathematica and "e" in Python corresponds to "10^" in Mathematica.

HankelH1[5, 3] // N == 0.0430284 - 1.90595 I

scipy.special.hankel1(5,3) == (0.04302843487704767-1.9059459538286736j)


HankelH1[16, 3] // N == 2.74882*10^-11 - 7.36866*10^8 I

scipy.special.hankel1(16,3) == (-1.0049113641502094e-08-736865544.0285703j)
POSTED BY: D D

I think in such cases, it is better to make a plot and compare:

Plot[{Log@Abs@Re@HankelH1[x, 3], Log@Abs@Im@HankelH1[x, 3]}, {x, 0, 
  20}, WorkingPrecision -> 100]

enter image description here

it has likely to do with some precision issues. For example compare:

HankelH1[16.0, 3]
N[HankelH1[16, 3]]
1.98576*10^-6 - 7.36866*10^8 I
2.74882*10^-11 - 7.36866*10^8 I

For the first example, machine precision is not enough to get the correct answer.

POSTED BY: Sander Huisman
Posted 11 months ago

Hello,

thank you very much for your reply and the nice idea with the plot. So, you clearly demonstrated that it has to do something with precision issues. Nevertheless, I wonder why in both cases, no result coincidences with the Python result. The physical problem I am solving suggests that the Python result is the correct one. However, this might only be the case because the physical problem I am solving was designed by a Python user, whereas I am a beginner in Python. Maybe I can find out which value is "true". I hope we all agree that the real part significantly differs, which causes problems in my further calculations.

POSTED BY: D D
Posted 11 months ago

I tried the online Keisan calculator. The result agrees with WL

hankelH1(16, 3)
(* 2.748824970048593080895E-11 -736865544.0285708677881i  *)

N[HankelH1[16, 3], 25]
(* 2.7488249700485930808946*10^-11 - 7.368655440285708677881434*10^8 I *)
POSTED BY: Rohit Namjoshi
Posted 11 months ago

Hello,

thank you for your reply. Then it seems that "N[HankelH1[16, 3], 25]" is the way to go. Nevertheless, the result from Python

scipy.special.hankel1(16,3) == (-1.0049113641502094e-08-736865544.0285703j)

where the real part significantly differs, leads to a more meaningful physical solution in my specific problem. Nevertheless, as I stated before, this might be only the case due to the fact that the problem was designed by a Python user.

In general, I trust more in Mathematica because there I feel more confident. However, I wonder if there is really a bug in Python for such a basic function? It is hard to believe for me. Maybe there are different conventions? However, for small n everything seems to be fine which makes me really mistrustful.

POSTED BY: D D

If scipy just uses machine precision arithmetic it might just simply not be enough to calculate with the correct precision.

Since the real part is determined by the Bessel function of the second kind:

N[BesselJ[16, 3], 25]

You can further check that in scipy.

POSTED BY: Sander Huisman
Posted 11 months ago

Ok, thank you!

Then I have now a problem to think about. I need to build up an at least 16x16 square matrix, whose elements are composed of HankelH1. Hereinafter, I need to find the roots of the determinant. Since my approach in Mathematica was too slow regarding the calculation time, I was told to use Python in order to accomplish this numerical problem. Now, Python is seemingly not accurate enough -- at least in my approach. Exporting the data from Mathematica would not help in this case, because I need the determinant as a function of a variable "k" in order to optimize it via "scipy.optimize.root". However, I think at this level I should directly contact the one who designed the task and happily, as I was sure, in Mathematica is, as always, everything fine! :-)

POSTED BY: D D
Posted 11 months ago

Dear all,

I just want to add further information: You can alternatively in Python use "mpmath". For example: mpmath.hankel1(16, 3) == (2.74882497004859e-11 - 736865544.028571j)

Hope this helps others who encounter the same issue.

Kind regards.

POSTED BY: D D
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract