# Integrate over one variable of a 2D interpolating function returned from NDEigensystem

I'm trying to implement the answer for "Integrate only one variable of a 2D interpolating function" (https://mathematica.stackexchange.com/a/161962/73672) but for interpolating functions returned from NDEigenSystem it isn't working.

First I solved for my required eigenfunctions and stored them as "funs":

ClearAll["Global*"];
ClearAll[vals, funs, schröd];
A = 0.025;
Subscript[V, 0] = 1;
d = 2;
schröd = -A*d^2 π*D[ψ[n, φ], {φ, 2}] +
A/(4 π) (φ*φ*ψ[n, φ] +
2 I*D[ψ[n, φ], {n, 1}] -
D[ψ[n, φ], {n, 2}]) -
Subscript[V,
0] ((Cos[2 π*d*n]) + Cos[φ] - 20) ψ[
n, φ];
Subscript[n, min] = -1/2; Subscript[n, max] = 1/2;
Subscript[φ, min] = -π;
Subscript[φ, max] = π;
Ω =
Rectangle[{Subscript[n, min],
Subscript[φ, min]}, {Subscript[n, max],
Subscript[φ, max]}];

{vals, funs} =
NDEigensystem[{schröd,
PeriodicBoundaryCondition[ψ[n, φ],
Subscript[φ, min] <= φ <=
Subscript[φ, max] && n == Subscript[n, max],
FindGeometricTransform[{{Subscript[n, min],
Subscript[φ, min]}, {Subscript[n, min],
Subscript[φ, max]}}, {{Subscript[n, max],
Subscript[φ, min]}, {Subscript[n, max],
Subscript[φ, max]}}][[2]]],
PeriodicBoundaryCondition[Exp[I 2 π n]*ψ[n, φ],
Subscript[n, min] <= n <= Subscript[n, max] && φ ==
Subscript[φ, max],
FindGeometricTransform[{{Subscript[n, min],
Subscript[φ, min]}, {Subscript[n, max],
Subscript[φ, min]}}, {{Subscript[n, min],
Subscript[φ, max]}, {Subscript[n, max],
Subscript[φ,
max]}}][[2]]]}, ψ, {n, φ} ∈
Rectangle[{Subscript[n, min],
Subscript[φ, min]}, {Subscript[n, max],
Subscript[φ, max]}], 2];
Plot3D[{Evaluate[Abs[{funs[[1]][n, φ]  }^2]]}, {n,
Subscript[n, min], Subscript[n, max]}, {φ,
Subscript[φ, min], Subscript[φ, max]},
PlotRange -> All,
PlotLabel ->
"Eig  fun    1 ", AxesLabel -> Automatic]


The answer to the other question was as follows:

da =
Flatten[
Table[
{t, tau, N @ Sin[2 (t + 3 tau)] Exp[-2 t - tau]},
{t, 0, 2, 2/100}, {tau, 0, 5, 5/100}],
1];
f = Interpolation @ da;
{{x1, x2}, {y1, y2}} = f["Domain"];
intx = Integrate[f[x, y], x] /. x -> x2;

nintx[y_?NumericQ] := Module[{x}, NIntegrate[f[x, y], {x, x1, x2}]];

Plot[nintx[y], {y, y1, y2}, PlotRange -> All]


However, trying to implement this myself gives an error after the third line, I believe because it reads n2 as 0.5:

f = funs[[1]];
{{n1, n2}, {φ1, φ2}} = f["Domain"];
intn = Integrate[f[n, φ], n] /. n -> n2;


General::ivar: 0.5 is not a valid variable.
Integrate::ilim: Invalid integration variable or limit(s) in 0.5.

• funs doesn't evaluate with your code. Couldn't find definition of function \[Psi] Aug 11 '20 at 15:55
• Where does funs come from? Have you checked that n doesn't have the value .5? Aug 11 '20 at 16:28
• I've made some edits to the code and it is returning funs now. @UlrichNeumann are you referring to the fact that I am returning $\psi$ rather than $\psi(n,\varphi)$? Aug 11 '20 at 16:50
• @m_goldberg funs is just what I'm storing the eigenfunctions returned from NDEigenSystem in. n is just a variable, but it looks like n2 gets set to 0.5 because that's the domain of the eigenfunction Aug 11 '20 at 16:53

Integrate won't integrate an InterpolatingFunction over an ElementMesh, which is what NDEigenSystem returns:

f@"ElementMesh"

(*  NDSolveFEMElementMesh[{{-0.5, 0.5}, {-3.14159,
3.14159}}, {NDSolveFEMQuadElement["<" 408 ">"]}]  *)


We can integrate with Simpson's rule, since the mesh consists of quadratic QuadElements which have a node at the midpoints of the edges of the rectangles and in this case are evenly spaced:

nn = f["Grid"][[All, 1]] // DeleteDuplicates // Sort;
pp = f["Grid"][[All, 2]] // DeleteDuplicates // Sort;

(* check even spacing *)
{dn} = nn // Differences // DeleteDuplicates

(*  {0.0625}  *)

simp = With[{y = Outer[f, nn, pp]},
(Total@y[[{1, -1}]] + 2 Total@y[[3 ;; -3 ;; 2]] + 4 Total@y[[2 ;; -2 ;; 2]])*
dn/3
];
intn = Interpolation@Transpose@{pp, simp}


Compare with NIntegrate:

nint = Table[
NIntegrate[f[n, phi],
Evaluate@Flatten@{n, First@f@"Domain"}], {phi,
f["Grid"][[All, 2]] // DeleteDuplicates // Sort}];
(simp - nint)/nint // Abs // Max

(*  2.00292*10^-15  *)


Try

intn[phi_?NumericQ] := NIntegrate[f[n, phi], {n, n1, n2}]

Table[{phi, intn[phi]}, {phi,Subdivide[\[CurlyPhi]1, \[CurlyPhi]2, 5]}]
(*{{-3.14159, -0.00205648 - 0.00607021 I}, {-1.88496, -0.0253336 -0.0747785 I}, {-0.628319, -0.142423 -0.420398 I}, {0.628319, -0.142423 -0.420397 I}, {1.88496, -0.0253317 -0.0747728 I}, {3.14159, -0.00201903 - 0.00595968 I}}*)
`

Result is complex!