2
$\begingroup$

I'm running Mathematica 13.0, student edition, on Windows. In physics, the total scattering cross-section of plane waves off a hard sphere of radius $1$ is

$$\sigma(k) = \frac{4 \pi}{k^2} \sum_{l=0}^\infty (2l+1) \left|\frac{j_l(k)}{h_l(k)} \right|^2$$

where $j_l$ is a spherical Bessel function and $h_l$ is a spherical Hankel function, both of the first kind.

I am hoping to evaluate this in Mathematica and to make a plot of $\sigma(k)$ against $k$ from $.1$ to $20$. However, I'm running into serious difficulties - Mathematica seems to hang, even if I truncate the upper limit of the sum from $\infty$ down to just $30$.

For example, for $k=1$,

SBJ = SphericalBesselJ;
SH1 = SphericalHankelH1;
(4 Pi)/1^2 NSum[(2 l + 1) Abs[SBJ[l, 1]/SH1[l, 1]]^2, {l, 0, 30}]

hangs for a quite long time before spitting the error "Summand (or its derivative) [complicated expression for summand] is not numerical at point l=15".

Truncating the upper limit to just $10$ seems to run instantaneously, but I expect that to give a poor approximation at larger values of $k$; $k=20$, say.

I was under the impression that spherical Bessel and Hankel functions are actually rather simple - for example, $h_0(k) = -i \frac{e^{ik}}{k}$, $h_1(k) = (-\frac{i}{k^2}-\frac{1}{k})e^{ik}$, and so on. That is, the Hankel functions are just complex exponentials multiplying complex polynomials in $\frac{1}{k}$. I could see trouble brewing for small $k$, but above I was taking $k=1$.


What's causing the massive slowdown in calculating the sum, even when truncating the upper bound of the sum to a small value like $l_{max}=30$? How might I speed this up?

$\endgroup$

2 Answers 2

4
$\begingroup$

I changed Abs[...]^2 to ... Conjugate[...] and NSum to Sum[...]//N

SBJ = SphericalBesselJ;
SH1 = SphericalHankelH1;
(4 Pi)/1^2   Sum[(2 l + 1)  # Conjugate[#] &[SBJ[l, 1]/SH1[l, 1]] , {l, 0, 30}] // N
(*10.6262 + 4.36409*10^-17 I*)
$\endgroup$
2
  • $\begingroup$ Thank you for the answer, this exactly solved my hanging issue, and I've now got a nice plot that fits with my rough intuition (decays from $4\pi$ to $2\pi$). Is there an easy way to see why this was so helpful? $\endgroup$
    – user196574
    Jul 7 at 7:03
  • 2
    $\begingroup$ @user196574 Thanks! Very often Abs[] shows numerically problems, that's why I try to avoid if possible. Don't know why NSum fails. $\endgroup$ Jul 7 at 7:14
4
$\begingroup$

Abs is no problem, if you avoid interpolation of NSum and calculate NSum at the exact integer l values with the help of NSumTerms .

SBJ = SphericalBesselJ;
SH1 = SphericalHankelH1;
(4 Pi)/1^2 NSum[(2 l + 1) Abs[SBJ[l, 1]/SH1[l, 1]]^2, {l, 0, 30}, 
  NSumTerms -> 31]

(*   10.6262   *)
$\endgroup$
6
  • $\begingroup$ Interesting. Why was there interpolation happening for a finite sum in the first place? That seems to be bad in every case I can think of. $\endgroup$
    – QuantumDot
    Jul 7 at 15:42
  • $\begingroup$ NSum only provides an approximation. "By default NSum uses 15 terms at the beginning before approximating the tail." If the tail is not sufficiently well-behaved, additional terms are required using the option NSumTerms ("number of terms to use before extrapolation"). The behavior of the tail is tested by observing the effect of varying the value of NSumTerms. $\endgroup$
    – Bob Hanlon
    Jul 7 at 16:00
  • $\begingroup$ @BobHanlon , NSum provides an aproximation for infinite sums or, if NSumTerms is lower than terms for finite sums. With l from 0 to 30 und NSumTerms =31, there is no tail. Regard rp = Reap[ NSum[(2 l + 1) Abs[SBJ[l, 1]/SH1[l, 1]]^2, {l, 0, 30}, NSumTerms -> 31, EvaluationMonitor :> Sow[l]]][[2, 1]] which gives all integer l. $\endgroup$
    – Akku14
    Jul 7 at 18:18
  • $\begingroup$ Not quite shure anymore, but in this cases does: {0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., \ 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., \ 29., 30.} $\endgroup$
    – Akku14
    Jul 7 at 18:26
  • 1
    $\begingroup$ @Akku14 - Look at ListPlot[Block[{k = 0}, NSum[ 1/(n + 1), {n, 0, 30}, NSumTerms -> #, EvaluationMonitor :> k++]; {#, k}] & /@ Range[15, 30], Frame -> True, FrameLabel -> {NSumTerms, "Number of terms"}] $\endgroup$
    – Bob Hanlon
    Jul 7 at 21:20

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.