# Getting the overlap of NDEigenfunctions of different problems

I am solving a number of Schrödinger eigenvalue problems for an array of different potentials, and I would like to calculate, in a quick, efficient and natural way, the overlap between the different eigenfunctions.

As a representative example, consider the ground-state eigenfunctions produced by a dipole and a quadrupole potential on the unit disk, as obtained with NDEigensystem, i.e.

{{e1}, {f}} = NDEigensystem[
{
-Laplacian[u[x, y], {x, y}] + 15 x u[x, y],
DirichletCondition[u[x, y] == 0, True]
},
u, {x, y} ∈ Disk[], 1
]


and

{{e2}, {g}} = NDEigensystem[
{
-Laplacian[u[x, y], {x, y}] + 15 (x^2 - y^2) u[x, y],
DirichletCondition[u[x, y] == 0, True]
},
u, {x, y} ∈ Disk[], 1
]


I would like to calculate the integral $\int_D f(\mathbf r)^* g(\mathbf r) \mathrm d \mathbf r$. The naive way to do this would be to just put them into a NIntegrate construct and and let Mathematica do its thing,

NIntegrate[Conjugate[f[x, y]] g[x, y], {x, y} ∈ Disk[]]


but this produces the error

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small

and generally it feels like there should be a better way to deal with this object that makes better use of the internal InterpolatingStructure mesh held by $f$ and $g$. I will be doing a lot of these as an internal part of an algorithm, so I would like to streamline this as much as possible. What are the best ways to achieve this?

## 1 Answer

One thing you can do is to use the same mesh for all computations:

Needs["NDSolveFEM"]
mesh = ToElementMesh[Disk[]];
{{e1}, {f}} =
NDEigensystem[{-Laplacian[u[x, y], {x, y}] + 15 x u[x, y],
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} \[Element] mesh, 1];
{{e2}, {g}} =
NDEigensystem[{-Laplacian[u[x, y], {x, y}] + 15 (x^2 - y^2) u[x, y],
DirichletCondition[u[x, y] == 0, True]},
u, {x, y} \[Element] mesh, 1];
NIntegrate[Conjugate[f[x, y]] g[x, y], {x, y} \[Element] mesh]
0.8704637241510582


You should note though the giving a mesh directly to NIntegrate will disable adaptive refinement in NIntegrate. If the input are interpolating functions that should not be too much of a problem; in fact, when the mesh in the interpolation function and the one generated and refined in NIntegrate` do not line up at the boundary that can result in problems. But your original mesh should be able to 'capture' the places in your domain where you expect action.

• But would adaptive refinement in NIntegrate actually help if you have interpolating functions? Wouldn't that just try to make up more information than is actually there? – Emilio Pisanty Feb 27 '18 at 14:29
• @EmilioPisanty, I added a note to address your comment. See if that helps. – user21 Feb 27 '18 at 14:37