1
$\begingroup$

I am trying to plot a normalized polar plot for the following function with different values of $a$

$$\left\lvert \sum_{n=1}^\infty i^n (2n+1) \frac {P_n^1(cos(\theta))}{\sqrt{\frac{\pi k a}{2}}[-H_{n+\frac{3}{2}}^2 (ka) + \frac{n+1}{ka}H_{n+\frac{1}{2}}^2(ka)]} \right\rvert^2 $$

where $P_n^1(cos(\theta))$ is the associated Legendre polynomial and $H_{n+\frac{3}{2}}^2 (ka)$ and $H_{n+\frac{1}{2}}^2(ka)$ are Hankel function of 2nd kind. Here $k=2\pi$ and the value of $a$ vary as $a = [.25, .05, 2, 10, 20]$

I can get plots upto $a=2$ with $n=100$ but for $a=10, 20$ I am having difficulty plotting. I am running into issues where machine precision is lost and the norm for Hankel function becomes too big. This is my attempt below

k = 2 \[Pi] ;
Pr[a_, m_, \[Theta]_] := 
Abs[Sum[I^n (2 n + 1) LegendreP[n, 1, Cos[\[Theta]]]/(
 Sqrt[(\[Pi] k a )/
   2] (-HankelH2[n + 3/2, k a ] + (n + 1)/(k a)
      HankelH2[n + 1/2, k a ])), {n, 1, m}]]^2;

ap05m = 2; normp05 = FindMaximum[Pr[.05, ap05m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
ap25m = 7 ; normp25 = FindMaximum[Pr[.25, ap25m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
a2m = 100; norm2 = FindMaximum[Pr[2.0, a2m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
a10m = 150; norm10 = FindMaximum[Pr[10.0, a10m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
a20m = 150; norm20 = FindMaximum[Pr[20.0, a20m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];

PolarPlot[{Pr[.05, ap05m,\[Theta]]/normp05[[1]], 
           Pr[.25, ap25m, \[Theta]]/normp25[[1]],
           Pr[2, a2m, \[Theta]]/norm2[[1]], 
           Pr[10, a10m, \[Theta]]/norm10[[1]], 
           Pr[20, a20m, \[Theta]]/norm20[[1]]}, {\[Theta], 0, 2 \[Pi]},
PolarAxes -> True, PlotRange -> Automatic, 
PolarGridLines -> Automatic, PolarTicks -> {"Degrees", Automatic}, 
PolarAxesOrigin -> {0, 1}, PlotLegends -> "Expressions"]

The norm for $a=10$ and $a=20$ are

a10m = 150; norm10 = FindMaximum[Pr[10.0, a10m, \[Theta]], {\[Theta], 0, 2 \[Pi]}]

${7.24673*10^{27}, {\theta -> 0.121351}}$

a20m = 150; norm20 = FindMaximum[Pr[20.0, a20m, \[Theta]], {\[Theta], 0, 2 \[Pi]}]

${1.11063*10^{79}, {\theta -> 0.0809369}}$

These normalization factors are too big and as a result I don't see the graphs for $a=10$ and $a=20$

This is what the graph is supposed to look like given by the professor

enter image description here

This is what my current graph looks like

enter image description here

I have tried $Chop[]$ function, but it didn't work. I would appreciate any help in plotting the normalized graph for $a = 10, 20$ on the same graph with $a = .05, .25, 2$.

$\endgroup$
  • 1
    $\begingroup$ The problem with your code is that you recompute the prefactor for each value of angle again and again. The idea would be to precompute the prefactors once and store them in the table. Or, if a minimal code modification is required you can use memoization option of MA reference.wolfram.com/language/tutorial/…. $\endgroup$ – yarchik Oct 11 at 8:21
  • $\begingroup$ @yarchik that's a great idea! I didn't know that! However, i still don't get the graphs for a = 10 and a = 20. $\endgroup$ – Rumman Oct 11 at 21:50
4
$\begingroup$

I am running into issues where machine precision is lost and the norm for Hankel function becomes too big.

Exactly. So you already know the answer: Don't use machine precision but ensure your computations are numerically sound by giving Mathematica enough digits to work with. I believe this answer about controlling precision is a good start.

In general, you have to understand that as soon as you say, e.g. 0.25, you told Mathematica that this is in machine precision. As long as you stick to infinitely accurate numbers like 1/4, you can tell Mathematica how precise it should evaluate numbers in certain algorithms later.

That being said, in the following I use ListPolarPlot and create the list of values myself. The adaptive plotting algorithms of Mathematica usually include many points to ensure your curve is smooth and calculating with higher precision takes a lot of time. We can speed this up by using ParallelTable but it will still require some minutes. I use the definitions of k and Pr you have already given.

as = {5/100, 1/4, 2, 10, 20};
ams = {2, 7, 100, 150, 150};
Block[{$MaxExtraPrecision = 100},
 norms = First@
     FindMaximum[Pr[##, \[Theta]], {\[Theta], 0, 2 \[Pi]}, 
      WorkingPrecision -> 30] & @@@ Transpose[{as, ams}]
 ]

First, note that I converted all your numbers to rationals and only in the final call to FindMaximum, I ask Mathematica to use a working precision of 30. Also, note that during the computation, Mathematica needs more precision and that's why we set $MaxExtraPrecision. If you don't set it, Mathematica will tell you that this setting is too low in a warning message.

Now, we can create the list of radii. For each curve, one list of radii, and we do this computation in parallel as it takes some minutes. With ParallelEvaluate, we distribute the setting of $MaxExtraPrecision to the parallel kernels.

ParallelEvaluate[$MaxExtraPrecision = 100];
data = ParallelTable[
     Pr[#1, #2, \[Theta]]/#3, {\[Theta], 0, 2 Pi, Pi/80}] & @@@ 
   Transpose[{as, ams, norms}];

As you can see, in above code blocks, I used @@@ and ## (or #1). Familiarize yourself how this works by looking at this example

f[##, {#2}] & @@@ {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}

and checking the documentation. Don't just copy the code.

After that, you can plot your data and recreate what your teacher had:

ListPolarPlot[data, Joined -> True, 
  PolarAxes -> True,
  PlotRange -> Automatic, 
  PolarGridLines -> Automatic, 
  PolarTicks -> {"Degrees", Automatic},
  PolarAxesOrigin -> {0, 1},
  PlotLegends -> (StringTemplate["a = ``\[Lambda]"] /@ N[as]), 
  PlotStyle -> {Blue, Red, Green, Magenta, Black}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks a lot for your response. I am going reiterate your explanation in my own words so that you can correct/check my understanding of your explanation. -Mathematica by default uses Machine Precision if the values are using decimal $\endgroup$ – Rumman Oct 15 at 3:47
  • $\begingroup$ - This was giving me precision issue. So you changed the list of radius values to "infinitely accurate precision" by using "/" and you controlled the precision later on in the calculation - I understand your notations. I played around with "@@@" and "##", "#1" Thanks for the example. That was helpful. - You used ParallelTable to evaluate the 5 sets of data $\endgroup$ – Rumman Oct 15 at 4:02
  • $\begingroup$ @Rumman Yes, exactly. There are different ways of controlling precision/accuracy (see 0.25``20 or 0.25`20 notation) but you can read about them in several posts here or in the documentation. Just look up Precision, Accuracy, N, SetAccuracy, ... and the standard options WorkingPrecision, AccuracyGoal, and the links therein. $\endgroup$ – halirutan Oct 15 at 6:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.