Using the result of NIntegrate as a potential in NDEigensysem

I am trying to find the eigenvalues and eigenstates of a Hamiltonian with a potential which cannot be analytically defined. In particular, it is the integral of the product between two functions:

V[r_,d_] := NIntegrate[Vint[r-s] X[s,d],{s,0,Infinity}]

where Vint[r_]:=(r^2 - c^2)/(r^2 + c^2) 1/(r^2 + b^2)^2 with b=0.0023 and c=0.4878, and X[s_,d_]:=1/(Sqrt[2 \[Pi]] d) Exp[-(1/2) s^2/d^2]

The integral seems to give no problem and I can successfully plot it as a function of r. Nevertheless, when I try to put it inside NDEigensystem with a fixed value of d:

Sol = NDEigensystem[{-D[u[r], {r, 2}] + V[r, 0.1] u[r], u[R] == 0,
u == 0}, u[r], {r, 0, R}, 2, Method -> {"PDEDiscretization" -> {"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}];

I get two possible outputs depending on the value of d:

NIntegrate: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option.

or

NIntegrate: The integrand -((3.98942 E^(-<<18>> s^2) (-0.237949+(r-s)^2))/((5.29*10^-6+(r-s)^2)^2 (0.237949 +(r-s)^2))) has evaluated to non-numerical values for all sampling points in the region with boundaries {{1.,0.}}.

Moreover, when I try to set a value for MinRecursion in the definition of V, I always get the second error.

Am I doing something wrong or it is just not possible to use such a function inside NDEigensystem?

• Please provide the complete minimal problem, so that readers can run it. Also, it might be better to use different symbols for V[r_,c_] and V[r-s] to avoid possible confusion. Also, try V[r_?NumericQ, c_?NumericQ]. Nov 10 '21 at 14:16
• What values of d are you using typically, and what value of MinRecursion? Nov 10 '21 at 15:17
• I tried with d between 0.01 and 1. Sorry, I forgot to write that. Nov 10 '21 at 15:19
• Values for R and MaxRecursion? Nov 10 '21 at 15:24
• R is about 1000, and I did not specify any MaxRecursion. I also tried with ?NumericQ as you suggested but I get the second error again. Nov 10 '21 at 16:21

This seemingly simple problem is numerically challenging due to the diversity of scales involved, or order R, d, b, and

Sqrt[-(2 b^3 (b + c)^2)/((b^2 - 2 b c - c^2) π)]

To gain an understanding of problem, begin by plotting the potential, V.

Plot[V[r, .1], {r, 0, R}, PlotRange -> All, ImageSize -> Large,
AxesLabel -> {r, V}, LabelStyle -> {15, Bold, Black}] Thus, there are two families of solutions, those trapped at the bottom of the deep well near r = 0, and those filling the entire range, {0 , R}.

As stated in the question, V is the integral of Vint[r-s] X[s, d]. Plotting Vint shows that it is delta-function-like for small b. Not surprizingly, Nintegrate encounters difficulties, unless MaxRecursion is very large. This explains the first error message in the question. The second error message may be due to a symbolic call by NDEigensystem to V during initialization. In any case, redefining V as

V[r_?NumericQ, d_?NumericQ] := NIntegrate[Vint[r - s] X[s, d], {s, 0, Infinity}]

is good practice. The integral over the delta-function-like Vint is.

Integrate[(r^2 - c1^2)/(r^2 + c1^2)/(r^2 + b1^2)^2, {r, -Infinity, Infinity},
Assumptions -> b1 > 0 && c1 > 0]

(* ((b1^2 - 2 b1 c1 - c1^2) π)/(2 b1^3 (b1 + c1)^2) *)

Now, first consider eigenfunctions trapped in the bottom of the well. By trial and error, I found that NDEigensystem does not find those modes unless the minimum value of the potential is shifted to near zero (and the eigenvalues shift accordingly). Solve

R = 0.02;
sol = NDEigensystem[{-D[u[r], {r, 2}] + (V[r, 0.1] -
1/(Sqrt[2 π] .1) ((b^2 - 2 b c - c^2) π)/(2 b^3 (b + c)^2)) u[r], u[R] == 0,
u == 0}, u[r], {r, 0, R}, 2,   Method -> {"PDEDiscretization" ->
{"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> R 10^-2}}}];
%[] + 1/(Sqrt[2 π] .1) ((b^2 - 2 b c - c^2) π)/(2 b^3 (b + c)^2)

(* {-5.10788*10^8, -5.10096*10^8} *)

Plot[Evaluate@sol[], {r, 0, R}, PlotRange -> All, ImageSize -> Large,
AxesLabel -> {r, u}, LabelStyle -> {15, Bold, Black}] Accurately computing solutions filling the range {0, R} is exceedingly slow using the code in the question. However, using the delta-function approximation to Vint makes such computations practical.

Vdel[r_?NumericQ, d_?NumericQ] :=
Piecewise[{{X[r, d] ((b^2 - 2 b c - c^2) π)/(2 b^3 (b + c)^2), r < 3}}, 0]

(Underflow messages otherwise occur for large values of r.)

R = 1;
sol1 = NDEigensystem[{-D[u[r], {r, 2}] + Vdel[r, 0.1] u[r], u[R] == 0, u == 0},
u[r], {r, 0, R}, 2, Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> R 10^-6}}}];
%[]

(* {37.7529, 150.994} *)

Plot[Evaluate@sol1[], {r, 0, R}, PlotRange -> All, ImageSize -> Large,
AxesLabel -> {r, u}, LabelStyle -> {15, Bold, Black}, PlotPoints -> 1000] The exceeding small MaxCellMeasure is required to resolve rapid oscillations near the well. For larger R, the wavelengths of the eigenfunctions away from the well become correspondingly longer, and the eigenfunctions approximated by (n Pi/R)^2, with n a positive integer.

• May be of advantage to evaluate the potential analytically @bbgodfrey . Vint[r_] = (r^2 - c^2)/(r^2 + c^2) 1/(r^2 + b^2)^2 /. {b -> 23/10000, c -> 4878/10000} // Together // ExpandAll // Apart; X[s_, d_] = 1/(Sqrt[2 \[Pi]] d) Exp[-(1/2) s^2/d^2]; vxe = Vint[r - s] X[s, d] // Expand; V[r_, d_] = Integrate[#, {s, 0, Infinity}, Assumptions -> {1/100 < d < 1, 0 < r < 1000}] & /@ vxe // ComplexExpand[Re@#, TargetFunctions -> {Re, Im}] & // Simplify[#, Assumptions -> {1/100 < d < 1, 0 < r < 1000}] & Nov 11 '21 at 6:20
• @Akku14 Very impressive. This eliminates the first error message in the question. By the way, I am surprised that each of the three terms in vxe must be integrated separately. Nov 11 '21 at 13:46