Taking the limit of a parametric function produced by ParametricNDSolveValue

Question

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.

Answer 1

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*)

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].