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
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.
ParametricNDSolveValue
only integrates out tor == 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$Odelta'[r]
gets set to zero whenr=a
, based on both the figures and the definition ofVhs
. Could you just integrate tor=3
and useOdelta[3]
as your answer? $\endgroup$