# How can I improve the speed of computing of numerical integration and eigenvalue?

First, let me explain what I am calculating.

I have the following 2D-equation:

$$w(x,y) = \sum_{m}^{M}\sum_{n}^{N}A_{mn}\sin\left(\frac{2m\pi x}{a}\right)\sin\left(\frac{2n\pi y}{b}\right)$$

where $$A_{mn}$$ are coefficients to be determined, $$a$$ and $$b$$ are constants.

After some differentiations and manipulations, I end-up with a linear system that I need to integrate and them get the eigenvalues.

The first matrix I'm calling it $$\mathbf{K}$$ and the another one is the $$\mathbf{M}$$ matrix. The eigenvalues equation is $$(\mathbf{K} -\lambda \mathbf{M})A_{mn}=\mathbf{0}$$.

$$\mathbf{M}$$ and $$\mathbf{K}$$ are $$(MN,MN)$$ in dimension and they are symmetrics.

Now, let's go to my code.

$$\mathbf{M}$$ can be defined as

mNx = Table[(8 π^2)/
a^2 m p Cos[(m π x)/((1/2) a)] Sin[(
n π y)/((1/2) b)] Cos[(p π x)/((1/2) a)] Sin[(
q π y)/((1/2) b)], {m, 1, nL}, {n, 1, nC}, {p, 1, nL}, {q, 1,
nC}];
mNx = Flatten@mNx;
mNx = ArrayReshape[mNx, {nL*nC, nL*nC}];


And $$\mathbf{K}$$ can be defined as

m11 = Table[(32 π^4)/
a^4 m^2 p^2 Sin[(m π x)/((1/2) a)] Sin[(
n π y)/((1/2) b)] Sin[(p π x)/((1/2) a)] Sin[(
q π y)/((1/2) b)], {m, 1, nL}, {n, 1, nC}, {p, 1, nL}, {q, 1,
nC}];
m11 = Flatten@m11;
m11 = ArrayReshape[m11, {nL*nC, nL*nC}];


I need to integrate each element of these matrices multiplied by $$N_x(x,y)$$ and $$D(x,y)$$ in $$\mathbf{M}$$ and $$\mathbf{K}$$, respectively, as can be seen below. These functions come from a FEM solution (they are interpolation function). To integrate then, I am using Gauss Quadrature, as follow.

nof = 15;
gqx = GaussianQuadratureWeights[nof, 0, a/2];
gqy = GaussianQuadratureWeights[nof, 0, b/2];
xi = gqx[[All, 1]];
yi = gqy[[All, 1]];
wix = gqx[[All, 2]];
wiy = gqy[[All, 2]];

mM = ConstantArray[0, {nL*nC, nL*nC}];
mK = ConstantArray[0, {nL*nC, nL*nC}];

For[i = 1, i <= nL*nC, i++,
For[j = i, j <= nL*nC, j++,

func[x_, y_] = mNx[[i, j]];
Ten1 = Table[
Nx[xi[[x]], yi[[y]]] func[xi[[x]], yi[[y]]], {x, 1, nof}, {y, 1,
nof}];
Ten2 = Table[wix[[x]] wiy[[y]], {y, 1, nof}, {x, 1, nof}];

Res = Dot[Flatten@Ten1, Flatten@Ten2];
If[Sqrt[Res^2] < 1, Res = 0, 0];
If[i == j, mM[[i, j]] = Res, mM[[i, j]] = Res; mM[[j, i]] = Res];

func[x_, y_] = m11[[i, j]];
Ten1 = Table[
D11[xi[[x]], yi[[y]]] func[xi[[x]], yi[[y]]], {x, 1, nof}, {y, 1,
nof}];
Ten2 = Table[wix[[x]] wiy[[y]], {y, 1, nof}, {x, 1, nof}];
Res = Dot[Flatten@Ten1, Flatten@Ten2];
If[Sqrt[Res^2] < 1, Res = 0, 0];
If[i == j, mK[[i, j]] = Res, mK[[i, j]] = Res; mK[[j, i]] = Res];

]
]


After calculating $$\mathbf{M}$$ and $$\mathbf{K}$$, I calculate the eigenvalue. I just need the smallest one.

mM = SparseArray[mM];
mK = SparseArray[mK];

eig = Eigenvalues[{mK, mM}, -1];


For $$M=N=10$$ my laptop it is spending

185.013 s

I am considering $$M=N=10$$ here from previous studies, but I need to perform an eigenvalue convergence.

To facilitate, it follows the whole code.

PS 1) As $$N_x(x,y)$$ and $$D(x,y)$$ are interpolating functions, in this code I am considering hypothetical functions.

PS 2) I am trying to improve my code because I need to perform 500+ times with different parameters that comes from the FEM code.

Thank you all in advance.

<< NumericalDifferentialEquationAnalysis;
AbsoluteTiming[
Nx[x_, y_] = x*y^2;
D11[x_, y_] = x^3*y^2;
nL = 10; nC = 10; a = 1; b = 1;
mNx = Table[(8 π^2)/
a^2 m p Cos[(m π x)/((1/2) a)] Sin[(n π y)/((1/
2) b)] Cos[(p π x)/((1/2) a)] Sin[(q π y)/((1/
2) b)], {m, 1, nL}, {n, 1, nC}, {p, 1, nL}, {q, 1, nC}];
mNx = Flatten@mNx;
mNx = ArrayReshape[mNx, {nL*nC, nL*nC}];
m11 = Table[(32 π^4)/
a^4 m^2 p^2 Sin[(m π x)/((1/2) a)] Sin[(n π y)/((1/
2) b)] Sin[(p π x)/((1/2) a)] Sin[(q π y)/((1/
2) b)], {m, 1, nL}, {n, 1, nC}, {p, 1, nL}, {q, 1, nC}];
m11 = Flatten@m11;
m11 = ArrayReshape[m11, {nL*nC, nL*nC}];
nof = 15;
gqx = GaussianQuadratureWeights[nof, 0, a/2];
gqy = GaussianQuadratureWeights[nof, 0, b/2];
xi = gqx[[All, 1]];
yi = gqy[[All, 1]];
wix = gqx[[All, 2]];
wiy = gqy[[All, 2]];
mM = ConstantArray[0, {nL*nC, nL*nC}];
mK = ConstantArray[0, {nL*nC, nL*nC}];
For[i = 1, i <= nL*nC, i++,
For[j = i, j <= nL*nC, j++, func[x_, y_] = mNx[[i, j]];
Ten1 =
Table[Nx[xi[[x]], yi[[y]]] func[xi[[x]], yi[[y]]], {x, 1,
nof}, {y, 1, nof}];
Ten2 = Table[wix[[x]] wiy[[y]], {y, 1, nof}, {x, 1, nof}];
Res = Dot[Flatten@Ten1, Flatten@Ten2];
If[Sqrt[Res^2] < 1, Res = 0, 0];
If[i == j, mM[[i, j]] = Res, mM[[i, j]] = Res; mM[[j, i]] = Res];
func[x_, y_] = m11[[i, j]];
Ten1 =
Table[D11[xi[[x]], yi[[y]]] func[xi[[x]], yi[[y]]], {x, 1,
nof}, {y, 1, nof}];
Ten2 = Table[wix[[x]] wiy[[y]], {y, 1, nof}, {x, 1, nof}];
Res = Dot[Flatten@Ten1, Flatten@Ten2];
If[Sqrt[Res^2] < 1, Res = 0, 0];
If[i == j, mK[[i, j]] = Res, mK[[i, j]] = Res;
mK[[j, i]] = Res];]];
mM = SparseArray[mM];
mK = SparseArray[mK];
eig = Eigenvalues[{mK, mM}, -1];][]

Out= 154.176

• Have you tried using NIntegrate instead of your hand-coded Gauss quadrature? The For loops often spell trouble for speed in Mathematica. – Marius Ladegård Meyer Apr 29 '20 at 21:49
• @MariusLadegårdMeyer, thank you for responding. Unfortunately, using NIntegrate is huge slow. Using just NIntegrate I waited 30 min+ and Abort due to long time. Using Method -> "Trapeizodal", a well known fast method, I waited 30 min+ and nothing... Would you have another advice? – Professor P. Cosmo Klunk Apr 30 '20 at 1:24
• You have nested For loops, each doing $M N$ operations, for $M^2 N^2$ in total. That is where to focus on optimising. Anything that does not depend on i and j shouldn't be inside those inner loops, and you need to cache results and use Compile where possible. Basically this code needs completely rewriting in Mathematica style, because this is not. I suspect Outer may also be very useful. – SPPearce Apr 30 '20 at 5:24
• For Eigensystem/values use Method -> "Arnoldi"` – user21 Apr 30 '20 at 7:31
• @KraZug thank you for the advices. I'll try them. – Professor P. Cosmo Klunk Apr 30 '20 at 13:06