2
$\begingroup$

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];][[1]]

Out[2]= 154.176
$\endgroup$
7
  • 2
    $\begingroup$ Have you tried using NIntegrate instead of your hand-coded Gauss quadrature? The For loops often spell trouble for speed in Mathematica. $\endgroup$ – Marius Ladegård Meyer Apr 29 '20 at 21:49
  • $\begingroup$ @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? $\endgroup$ – Professor P. Cosmo Klunk Apr 30 '20 at 1:24
  • 2
    $\begingroup$ 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. $\endgroup$ – SPPearce Apr 30 '20 at 5:24
  • $\begingroup$ For Eigensystem/values use Method -> "Arnoldi" $\endgroup$ – user21 Apr 30 '20 at 7:31
  • $\begingroup$ @KraZug thank you for the advices. I'll try them. $\endgroup$ – Professor P. Cosmo Klunk Apr 30 '20 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.