# Using NIntegrate with interpolating function

I am relatively new to Mathematica and have been trying to use the NDEigensystem command to work with some quantum systems. This is actually inline with a previous question I had asked about a week back. I am able to use the NDEigensystem command to get the eigen values and eigen functions, here is the basic model I'm testing this with:

m2 = 0.5;
ℏ = 1;
w = 0.5;
\[ScriptCapitalO]2 = -ℏ^2/(2 m2) Laplacian[u[x, y], {x, y}] +
1/2 m2 w^2 (x^2 + y^2) u[x, y];

{vals, funs} =
NDEigensystem[{\[ScriptCapitalO]2,
DirichletCondition[u[x, y] == 0, True]},
u[x, y], {x, -10, 10}, {y, -10, 10}, 28,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.5}}}}];


Using the answer to my previous question I am able to get the eigenfunctions as functions as shown below; however, I am still unable to use them with NIntegrate for some reason:

funs2 = Function[{x, y}, #] & /@ funs;
\[Psi]1 = funs2[[1]];
\[Psi]2 = funs2[[2]];
Ans = NIntegrate[\[Psi]1 x \[Psi]2,{x,-10,10},{y,-10,10}]


The last command doesn't return anything and just gives me back the command line again. I am not certain where I am going wrong. Would be very grateful for any help.

Have a look at the documentation. This is from the ref page of NDEigensystem

{vals, funs} =
NDEigensystem[-Laplacian[u[x], {x}], u[x], {x, 0, \[Pi]}, 4]

NIntegrate[#^2, {x, 0, \[Pi]}] & /@ funs

(* {1., 1., 0.999995, 1.} *)


Note that the argument u[x] to NDEigensystem tells NDEigensystem that the resulting interpolating functions will also have the independent variable x as in InterpolatingFunction[data][x]. If u were given as an argument the result would be InterpolatingFunction[data].

• Thanks a lot, sorry for the late response, had fallen a little ill. This definitely works Aug 4, 2020 at 3:38
• @rahulmenon, no worries. I hope you are feeling better. Aug 4, 2020 at 4:59
Clear["Global*"]

m2 = 1/2;
ℏ = 1;
w = 1/2;
\[ScriptCapitalO]2 = -ℏ^2/(2 m2) Laplacian[u[x, y], {x, y}] +
1/2 m2 w^2 (x^2 + y^2) u[x, y];

{vals, funs} =
NDEigensystem[{\[ScriptCapitalO]2, DirichletCondition[u[x, y] == 0, True]},
u[x, y], {x, -10, 10}, {y, -10, 10}, 28,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" -> \
{"MaxCellMeasure" -> 0.5}}}}];

funs2 = Function[{x, y}, #] & /@ funs;

ψ1[x_?NumericQ, y_?NumericQ] := funs2[[1]][x, y];
ψ2[x_?NumericQ, y_?NumericQ] := funs2[[2]][x, y];


Verifying that the functions evaluate

#[1, 1] & /@ {ψ1, ψ2}

(* {-0.21977, 0.104243} *)


The integrand is

Plot3D[
ψ1[x, y]*ψ2[x, y], {x, -10, 10}, {y, -10, 10},
PlotRange -> All, AxesLabel -> Automatic,
PlotPoints -> 50, MaxRecursion -> 3]


The min and max are

#[{ψ1[x, y]*ψ2[x, y], -5 < x < 5, -5 < y < 5},
{x, y}] & /@ {NMinimize, NMaximize}

(* {{-0.0482587, {x -> -0.397854, y -> 1.35926}}, {0.0482587, {x -> 0.397854,
y -> -1.35926}}} *)


From the symmetry, the integral is expected to be near zero

Ans = NIntegrate[ψ1[x, y]*ψ2[x, y],
{x, -10, 10}, {y, -10, 10}] // Quiet

(* 6.8028*10^-13 *)


Alright, I seem to have a solution that works for now, but I do planning on using a loop and this might get a little hard to incorporate into the same, so if anyone has any better ideas please do let me know. I explicitly defined an Integrand as a function of x and y and was able to integrate the functions after that

Integrand[x_,y_] = \[Psi]1[x,y] x \[Psi]2[x,y]
NIntegrate[Integrand[x, y], {x, -10, 10}, {y, -10, 10}]
`

Using this method Mathematica returned a value for this numerical integration.