5

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.

user21
  • 39,710
  • 8
  • 110
  • 167
Xiaoyu Liu
  • 113
  • 8
  • 1
    Can you post a working, minimal example of what you want? – user21 Feb 08 '18 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. – Xiaoyu Liu Feb 08 '18 at 16:34
  • Do you want that x in the NIntegrate of your specific example? – Chris K Feb 08 '18 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. – Xiaoyu Liu Feb 08 '18 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. – Xiaoyu Liu Feb 08 '18 at 17:10

2 Answers2

7

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["NDSolve`FEM`"]
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`}
user21
  • 39,710
  • 8
  • 110
  • 167
  • Thank you, but I also need grid information to perform an integration. How to extract that? – Xiaoyu Liu Feb 08 '18 at 19:39
  • @XiaoyuLiu, see update (and Documentation for ElementMesh) – user21 Feb 08 '18 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? – Xiaoyu Liu Feb 09 '18 at 09: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]? – Xiaoyu Liu Feb 09 '18 at 09:48
  • @XiaoyuLiu, you can use Head[u[x]]["ValuesOnGird"]. – user21 Feb 09 '18 at 09:49
  • Sorry, when and where should I use that? Out of NDEigensystem mathematica doesn't know the definition of u[x] any more. – Xiaoyu Liu Feb 09 '18 at 09:52
  • @XiaoyuLiu, see update. – user21 Feb 09 '18 at 09:54
  • This is good, but where is "ElementMesh"? I think I need at least step size or the length of the segments. And by the way I tested on discrete dot product, it is at least two orders faster. However, the interpolating process I used to extract the discrete numbers are slow (Using Table). If I could get discretized eigenfunctions directly, that would be the best solution. – Xiaoyu Liu Feb 09 '18 at 10:06
5

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}
*)
Michael E2
  • 235,386
  • 17
  • 334
  • 747