Getting the overlap of NDEigenfunctions of different problems

Question

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?

Answer 1

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

Needs["NDSolve`FEM`"]
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.