# Is there a way to let NDEigensystem output discrete eigenfunctions rather than InterpolatingFunctions?

I am new to Mathematica and I encountered a speed problem.

I want to do tons of integrations of multiplications of two eigenfunctions (or matrix elements in physical terms). The speed now is quite unsatisfying.

As I only require numerical results, I would like to know how to speed up this process. The most efficient way I could think of is to get the discretized eigenfunctions. But I would also like to know how and at which stage discretization of the eigenfunctions could speed up the whole calculation most.

For example, the problem is

{eige, eigv}=NDEigensystem[SomeSystem];
result[[i,j]]=NIntegrate[eigv[[i]] eigv[[j]]];


Here i and j are different solutions corresponding to eige[[i]] and eige[[j]] respectively.

The stages I could think of are:

1. When solving the differential eigensystem using NDEigensystem, (get discretized eigv at NDEigensystem)
2. When solved, transform the eigenfunction of the form InterpolatingFunction into discrete ones, (discretize the eigv which is now in InterpolatingFunction before NIntegration and then integrate by dot product.)
3. When using NIntegration.

And, I would be very thankful if you could also teach me how to do it in these ways.

A specific example is below:

num = 10;
AbsoluteTiming[{eige, eigv} = NDEigensystem[-Laplacian[u[x], {x}], u[x], {x, 0, \[Pi]}, num];]
AbsoluteTiming[Sum[NIntegrate[eigv[[i]] x eigv[[j]], {x, 0, \[Pi]}], {i, 1,num}, {j, 1, num}]]


When num is as large as 100, the result is super slow.

Editted

After some efforts trying to understand the stages taken by NDEigensystem, I give up letting the function itself output discretized data. Although it seems this function contains a normalization process (or stage) after interpolation stage, I do not think I could realize arbitrary integration as fast as the internal function do. I found we can discretize the solution by hand. See below

Clear[x, u];
num = 40;
delta = 0.01;
{eige, eigv} = NDEigensystem[-Laplacian[u[x], {x}], u[x], {x, 0, \[Pi]}, num];
Sum[NIntegrate[eigv[[i]] x eigv[[j]], {x} \[Element] (eigv[[1]] /. x -> "ElementMesh")], {i, 1, num}, {j, 1, num}] // Timing
Sum[NIntegrate[eigv[[i]] x eigv[[j]], {x, 0, \[Pi]}], {i, 1, num}, {j, 1, num}] // Timing
AbsoluteTiming[
xRange = Range[0, \[Pi], delta];
eigvv = Table[eigv, {x, xRange}];
Print[Dimensions[eigvv], Dimensions[xRange]];
Sum[eigvv[[All, i]].(eigvv[[All, j]] xRange) delta, {i, 1, num}, {j, 1, num}]]
(*
{30.5625, 61.5422}
some warnings
{41.1875, 61.5243}
{315,40}{315}
{0.158992, 61.8818}
*)


The improvement is almost two orders.

• Can you post a working, minimal example of what you want? Commented Feb 8, 2018 at 16:18
• Hi, I have added a specific example. The resulted eigenfunction is rather smooth, so I think this code should have a better performance even though the integration is performed many times. Commented Feb 8, 2018 at 16:34
• Do you want that x in the NIntegrate of your specific example? Commented Feb 8, 2018 at 16:40
• Yeah, I hope there is a general solution for this. That x is an example, and you can even generalize to any function f[x] to insert between eigenfunctions and also the eigenfunction could be a vector itself. Commented Feb 8, 2018 at 16:44
• I tried the solution from another post but it is weird that it is not actually calculating. Sum[NIntegrate[eigv[[i]] x eigv[[j]], x \[Element] eigv[[2]]["ElementMesh"]], {i, 1, num}, {j, 1, num}]. Besides, it seems that from the FullForm, eigenfunctions are already discretized. Commented Feb 8, 2018 at 17:10

You can extract the values from the grid:

{vals, funs} =
NDEigensystem[-Laplacian[u[x], {x}], u, {x, 0, \[Pi]}, 10];
funs[[3]]["ValuesOnGrid"]


Extract the mesh with:

funs[[3]]["ElementMesh"]


If you want to stay with your original formulation u[x] than you can use:

Head[funs[[3]][x]]["ValuesOnGird"]


Or another alternative:

Needs["NDSolveFEM"]
mesh = ToElementMesh[Line[{{0}, {Pi}}]];
num = 10;
{vals, eigv} =
NDEigensystem[-Laplacian[u[x], {x}], u, Element[{x}, mesh], num];
funs[[3]]["ValuesOnGrid"]
AbsoluteTiming[
Sum[NIntegrate[eigv[[i]][x] x eigv[[j]][x], Element[{x}, mesh]], {i,
1, num}, {j, 1, num}]]
{1.091801, 19.03787410869329}

• Thank you, but I also need grid information to perform an integration. How to extract that? Commented Feb 8, 2018 at 19:39
• @XiaoyuLiu, see update (and Documentation for ElementMesh) Commented Feb 8, 2018 at 19:49
• By the way in my example, I used u[x] instead of u as the second argument of NDEigensystem and in that way I cannot extract "ValuesOnGrid" and "ElementMesh". The result is always in the form of interpolating function. Can you explain me why is that? Commented Feb 9, 2018 at 9:41
• I tested a little bit on this difference and found they are not the same. When using u[x] it is pretty much like a function. The most important difference is when you want to extract discrete points. Try the following:{vals, funs} = NDEigensystem[-Laplacian[u[x], {x}], u, {x, 0, \[Pi]}, 10]; {eige, eigv} = NDEigensystem[-Laplacian[u[x], {x}], u[x], {x, 0, \[Pi]}, num]; vec1 = Table[eigv[[3]], {x, 0, \[Pi], 0.01}] vec2 = Table[funs[[3]], {x, 0, \[Pi], 0.01}]. Now the question is what is the correct way to extract ValuesOnGrid from my original definition with u[x]? Commented Feb 9, 2018 at 9:48
• @XiaoyuLiu, you can use Head[u[x]]["ValuesOnGird"]. Commented Feb 9, 2018 at 9:49

In response to a comment:

Thank you, but I also need grid information to perform an integration....

There are built-in methods for integrating over an ElementMesh, which are automatically invoked if you use the ElementMesh as the integration domain:

Clear[x, u];
num = 40;
AbsoluteTiming[{eige, eigv} = NDEigensystem[-Laplacian[u[x], {x}], u[x], {x, 0, π}, num];]
AbsoluteTiming[
Sum[
NIntegrate[eigv[[i]] x eigv[[j]], {x} ∈ (eigv[[1]] /. x -> "ElementMesh")],
{i, 1, num}, {j, 1, num}]]
(*
{0.019178, Null}
{21.7475, 61.5422}
*)


Compare with approach in the OP:

AbsoluteTiming[{eige, eigv} = NDEigensystem[-Laplacian[u[x], {x}], u[x], {x, 0, π}, num];]
AbsoluteTiming[
Sum[NIntegrate[eigv[[i]] x eigv[[j]], {x, 0, π}],
{i, 1, num}, {j, 1, num}]]
(*
{0.018287, Null}

...NIntegrate warnings...

{35.0064, 61.5243}
*)