Message Boards Message Boards

GROUPS:

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

Posted 1 month ago
455 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
Answer
9 Replies

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

Posted 1 month 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
Answer

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 1 month 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
Answer
Posted 1 month 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 1 month 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
Answer

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 1 month 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
Answer
Posted 27 days 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
Answer
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