1
$\begingroup$

I'm working on a project for hard-sphere scattering, and to solve the for the phase function $\overline{\delta}_l(k,r)$ I use the variable phase equation

$$ \frac{d\overline{\delta}_l(r,k)}{dr} = -\frac{U(r)}{k} \left[\hat{j}_l(kr)\cos\overline{\delta}_l(r,k) - \hat{n}_l(kr)\sin\overline{\delta}_l(r,k)\right]^2 $$

To find the phase function you integrate the above equation. In Mathematica, based on the specifications of the project, I wrote the following

Vhs[V0_, a_, r_] := Piecewise[{{V0, Abs[r] <= a}, {0, Abs[r] > a}}]; 
(*potential*)
U[V0_, a_, r_] := (2*mu/hbar^2)*Vhs[V0, a, r]
mathSoln = ParametricNDSolveValue[{Odelta'[r] == -U[5000, 3, r]/
k*(r*k*SphericalBesselJ[l, r*k]*Cos[Odelta[r]] - 
r*k*SphericalBesselY[l, r*k]*Sin[Odelta[r]])^2, Odelta[0.1] == 0}, 
Odelta, {r, 0.1, 6}, {k, l}]; (*integrate variable phase eqn*)
k = {0.1, 1, 5, 10}; (*specify k vals*)
Table[Plot[Evaluate[Table[mathSoln[k[[i]], l][r], {l, 0, 5, 1}]], 
{r, 0.1, 6}, PlotRange -> Full, Frame -> True, 
PlotLabel -> "k = " <> ToString[k[[i]]], 
PlotStyle -> {{Dashed, Red}, {Thin, Orange}, {Dashing[Small], 
Yellow}, {Thick, Green}, {Dotted, Blue}, {Thin, Purple}}, 
PlotLegends -> {"l = 0", "l = 1", "l = 2", "l = 3", "l = 4", "l = 5", 
FrameLabel -> {"r", "phase eqn"}}], {i, 1, 4, 1}] (*plot for specific l and 
k*)

The above gives me the following plots 4 plots for 4 different values of k. l is varied from 0 to 5. The x-axis is r and the y axis is values of the phase equation.

Now, for the next part I need to find the phaseshift $\delta_l(k)$ from the phase function $\overline{\delta}_(r,k)$ by using

$$\delta_l(k) = \lim_{r \rightarrow \infty} \overline{\delta}_l(r,k)$$

What I wanted to do was write something like this

Clear[k];
Limit[Evaluate[Table[mathSoln[k, l][r], {l, 0, 5, 1}]], {r -> \[Infinity]}]

In order to take the limit of the parametric function generated from ParametricNDSolveValue. However, this does not evaluate the limit at all. Is there a way to take the limit of parametric functions produced by ParametricNDSolveValue? I believe that the issue probably lies in how I am writing this, or that I probably need to use something else other than Limit[] to take the limit in this case.

$\endgroup$
  • 2
    $\begingroup$ Your ParametricNDSolveValue only integrates out to r == 6, so the answer is no, because extrapolation would be used (if the limit were computed numerically), which is a completely unreliable way to evaluate the limit. My answers here] and here use changes of variables to allow integration to infinity. I still wonder if numerical error might be great at infinity. Perhaps you can set $\frac{d\overline{\delta}_l(r,k)}{dr}$ equal to zero and solve for potential limiting values that way. $\endgroup$ – Michael E2 Nov 8 '18 at 3:26
  • $\begingroup$ Looks like Odelta'[r] gets set to zero when r=a, based on both the figures and the definition of Vhs. Could you just integrate to r=3 and use Odelta[3] as your answer? $\endgroup$ – Chris K Nov 8 '18 at 4:41
3
$\begingroup$

I corrected the code so that the result was visible in all versions. To find the limit, we can construct another parametric function and make sure that the values of the function for r = 6 and r = 12 coincide. Therefore, you can use as a limit mathSoln12[k,l][12]

Vhs[V0_, a_, r_] := Piecewise[{{V0, Abs[r] <= a}, {0, Abs[r] > a}}];
(*potential*)
U[V0_, a_, r_] := (2*mu/hbar^2)*Vhs[V0, a, r]
mu = 1; hbar = 1; mathSoln = 
 ParametricNDSolveValue[{Odelta'[
     r] == -U[5000, 3, r]/
      k*(r*k*SphericalBesselJ[l, r*k]*Cos[Odelta[r]] - 
        r*k*SphericalBesselY[l, r*k]*Sin[Odelta[r]])^2, 
   Odelta[0.1] == 0}, 
  Odelta, {r, 0.1, 6}, {k, 
   l}];(*integrate variable phase eqn*); Table[
 Plot[Evaluate[Table[mathSoln[k, l][r], {l, 0, 5, 1}]], {r, 0.1, 6}, 
  PlotRange -> Full, Frame -> True, PlotLabel -> Row[{"k = ", k}], 
  PlotStyle -> {{Dashed, Red}, {Thin, Orange}, {Dashing[Small], 
     Yellow}, {Thick, Green}, {Dotted, Blue}, {Thin, Purple}}, 
  PlotLegends -> {"l = 0", "l = 1", "l = 2", "l = 3", "l = 4", 
    "l = 5", FrameLabel -> {"r", "phase eqn"}}], {k, {0.1, 1, 5, 
   10}}] (*plot for specific l and k*)

fig1 Determine the function in double scale

mathSoln12 = 
  ParametricNDSolveValue[{Odelta'[
      r] == -U[5000, 3, r]/
       k*(r*k*SphericalBesselJ[l, r*k]*Cos[Odelta[r]] - 
         r*k*SphericalBesselY[l, r*k]*Sin[Odelta[r]])^2, 
    Odelta[0.1] == 0}, Odelta, {r, 0.1, 12}, {k, l}];

Calculate the limit

dkl12 = 
 Table[{k, l, mathSoln12[k, l][12]}, {l, 0, 5, 1}, {k, {0.1, 1, 5, 10}}]


{{{0.1, 0, -0.299}, {1, 0, -2.99}, {5, 0, -14.95}, {10, 
   0, -29.8998}}, {{0.1, 1, -0.00846089}, {1, 1, -1.74196}, {5, 
   1, -13.446}, {10, 1, -28.3625}}, {{0.1, 2, -0.0000519763}, {1, 
   2, -0.834287}, {5, 2, -12.009}, {10, 2, -26.8586}}, {{0.1, 
   3, -1.48589*10^-7}, {1, 3, -0.286878}, {5, 3, -10.6398}, {10, 
   3, -25.3882}}, {{0.1, 4, -2.72718*10^-10}, {1, 4, -0.0597636}, {5, 
   4, -9.33914}, {10, 4, -23.9515}}, {{0.1, 5, -3.35067*10^-13}, {1, 
   5, -0.00708}, {5, 5, -8.10836}, {10, 5, -22.5487}}} 

Calculate the difference of limits at `r=12` and `r=6`


Table[{k, l, mathSoln12[k, l][12] - mathSoln12[k, l][6]}, {l,
   0, 5, 1}, {k, {0.1, 1, 5, 10}}]

Out[]= {{{0.1, 0, 0.}, {1, 0, 0.}, {5, 0, 0.}, {10, 0, 0.}}, {{0.1, 
   1, 0.}, {1, 1, 0.}, {5, 1, 0.}, {10, 1, 0.}}, {{0.1, 2, 0.}, {1, 2,
    0.}, {5, 2, 0.}, {10, 2, 0.}}, {{0.1, 3, 0.}, {1, 3, 0.}, {5, 3, 
   0.}, {10, 3, 0.}}, {{0.1, 4, 0.}, {1, 4, 0.}, {5, 4, 0.}, {10, 4, 
   0.}}, {{0.1, 5, 0.}, {1, 5, 0.}, {5, 5, 0.}, {10, 5, 0.}}}

Therefore, the function mathSoln12[k, l][12] reaches a constant value already at r = 6. With a margin, you can take as a limit mathSoln12[k, l][12].

$\endgroup$
  • $\begingroup$ Thanks a lot! With this, I was able to move towards the direction I needed to (plotting the limit for a range of k and decreasing the step size) and compare it to the analytic results. $\endgroup$ – Jomy Blue Nov 8 '18 at 17:25

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.