0
$\begingroup$

I've the following integral:

$$\int_{-a/2}^{a/2}Cf(x)\sin{\frac{m\pi(x+0.5a)}{a}}\sin{\frac{p\pi(x+0.5a)}{a}}dx$$

where $$f(x) =53.62\cos^4\left(\frac{7\pi\sqrt{x^2}}{9}\right)+3.04\sin^4\left(\frac{7\pi\sqrt{x^2}}{9}\right)+2.54\sin^2\left(\frac{14\pi\sqrt{x^2}}{9}\right)$$

$$m=1,2,\dots,M$$ $$p=1,2,\dots,P$$. $$C=\text{constant}$$

I'm using the following code

f[x] = 53.62 Cos[(7 \[Pi] Sqrt[x^2])/9]^4 + 
   3.04 Sin[(7 \[Pi] Sqrt[x^2])/9]^4 + 
   2.54 Sin[(14 \[Pi] Sqrt[x^2])/9]^2;
nL = 20; nC = 20;
K = ConstantArray[0, {nL*nC, nL*nC}];
a=1;b=1;
Do[
  Do[
   Do[
    Do[

     k = (m - 1) nC + n;
     l = (p - 1) nC + q;

     If[n == q, 
      K[[k, l]] = (\[Pi]^2 b m^2 p^2)/(2 a^4)
         NIntegrate[
         f[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(
           p \[Pi] (x + 0.5 a))/a], {x, -a/2, a/2}, 
         Method -> "Trapezoidal"], K[[k, l]] = 0]

     , {q, 1, nC}],
    {p, 1, nL}],
   {n, 1, nC}],
  {m, 1, nL}];

Does anyone know another faster and accurate approach?

Thank you all in advance.

Improve this question
$\endgroup$
13
  • $\begingroup$ Did you mean $n\pi$ rather than $p\pi$? $\endgroup$
    – JimB
    Feb 12, 2020 at 23:19
  • $\begingroup$ Hello, @JimB! Oh, I saw, I edited. $\endgroup$
    – Joaquim José da Silva Xavier
    Feb 12, 2020 at 23:30
  • $\begingroup$ The integral can be calculated symbolically. First, change $\sqrt{x^2}$ to $x$ three places. Then change $0.5a$ to $a/2$ two places. Then do the integral as 3 integrals with different integrands. Do each integration only once, outside the loops, of course. $\endgroup$
    – LouisB
    Feb 13, 2020 at 0:02
  • 2
    $\begingroup$ By the way, capital $K$ is a reserved word, so using it as a variable name could cause problems. $\endgroup$
    – LouisB
    Feb 13, 2020 at 0:06
  • 1
    $\begingroup$ Notice that your integrals depend on $m$ and $p$, but you recalculate them every time you change $q$. Also, in the integrals are symmetric in $m$ and $p$, so do not recalculate. It should be faster if you create a symmetric array of the values of the integrals, then use those matrix elements in q-p-n-m loop. $\endgroup$
    – LouisB
    Feb 13, 2020 at 1:11

4 Answers 4

Reset to default
2
$\begingroup$

Here is a symmetric array example that will reduce the total integration time. For simplicity this example is based on the $f(x)$ in your original post. 

Clear[f, amat, m, p]

f[x] = 53.62 Cos[(7 π Sqrt[x^2])/9]^4 + 
   3.04 Sin[(7 π Sqrt[x^2])/9]^4 + 
   2.54 Sin[(14 π Sqrt[x^2])/9]^2;
nL = 20;

amat = SymmetrizedArray[{m_, p_} :> NIntegrate[f[x] 
      Sin[(m π (x + 0.5 a))/a] Sin[(p π (x + 0.5 a))/a],
    {x, -a/2, a/2},  Method -> "Trapezoidal" ],
  {nL, nL}, Symmetric[{1, 2}]]

The code creates a symmetric array amat. The code should be executed before the nested $q-p-n-m-$loop. Then, inside the nested loop, simply use amat[ m, p ] instead of NIntegrate[ ... ]. amat must be recalculated only when f[x] changes.

We can see the elements of amat by evaluating amat // Normal // MatrixForm. We notice that half of the elements are almost zero.

Note that if your f[x] is always a linear combination of the same functions every time, and I think it must be, you should find a way to apply NIntegrate to those functions outside of all the loops and then calculate linear combinations of the integrals inside the nested loops where the coefficients are known.

Improve this answer
$\endgroup$
2
$\begingroup$

For some weird reason that I cannot understand at the moment, this one gets rid of the error messages and is this faster.

amat = SymmetrizedArray[{m_, p_} :> Plus[
     NIntegrate[
      f[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(p \[Pi] (x + 0.5 a))/
         a], {x, -a/2, 0}, 
      Method -> {"GaussKronrodRule", "Points" -> 7}],
     NIntegrate[
      f[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(p \[Pi] (x + 0.5 a))/
         a], {x, 0, a/2}, 
      Method -> {"GaussKronrodRule", "Points" -> 7}]
     ], {nL, nL}, Symmetric[{1, 2}]];

Should also be much more accurate (no time to check that).

Improve this answer
$\endgroup$
1
$\begingroup$

Your integral can be computed in analytic form. No need for numerics. Then you can tabulated it very easy.

Notice that the following definition

g[1]=53.62;
g[2]=3.04;
g[3]=2.54;
fp[x_]:=g[1] Cos[(7π)/9 x]^4 + g[2] Sin[(7π)/9 x]^4+ g[3] Sin[(7π)/9 2 x]^2;

is equivalent to your function.

I just saw an identical idea in the comments section. All credits to LouisB.

Improve this answer
$\endgroup$
0
$\begingroup$ 
Needs["NDSolve`FEM`"];
G = 6.894745 10^9;
E1 = 26.25 G;
E2 = 1.49 G;
G12 = 1.04 G;
nu12 = 0.28;
nu21 = (E2*nu12)/E1;
t = 0.0050 .0254;
a = 1; b = 1;
u0 = .01;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}};
Qon = Inverse[Son];

Do[
 Do[
  angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0, 
     angle1}, {-angle0, -angle1}, {angle0, 
     angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0, 
     angle1}, {-angle0, -angle1}, {angle0, 
     angle1}, {-angle0, -angle1}, {angle0, angle1}};

  num = Dimensions[angles][[1]];
  h = num*t;
  pos = Table[0, num + 1];
  pos[[1]] = -h/2;
  For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];
  mA = ConstantArray[0, {3, 3}];
  mB = ConstantArray[0, {3, 3}];
  mD = ConstantArray[0, {3, 3}];

  For[i = 1, i <= num, i++,
   T0 = angles[[i, 1]] ;
   T1 = angles[[i, 2]] ;
   theta[x_] := ((2/a) (T1 - T0) Sqrt[x^2] + T0) (Pi/180);
   m = Cos[theta[x]];
   n = Sin[theta[x]];
   Q11 = Qon[[1, 1]];
   Q12 = Qon[[1, 2]];
   Q22 = Qon[[2, 2]];
   Q66 = Qon[[3, 3]];
   Qxx = m^4*Q11 + n^4*Q22 + 2*m^2*n^2*Q12 + 4*m^2*n^2*Q66;
   Qyy = n^4*Q11 + m^4*Q22 + 2*m^2*n^2*Q12 + 4*m^2*n^2*Q66;
   Qxy = m^2*n^2*Q11 + m^2*n^2*Q22 + (m^4 + n^4)*Q12 + -4*m^2*n^2*Q66;
   Qss = m^2*n^2*Q11 + m^2*n^2*Q22 - 2*m^2*n^2*Q12 + (m^2 - n^2)^2*Q66;
   Qxs = m^3*n*Q11 - m*n^3*Q22 + (m*n^3 - m^3*n)*Q12 + 
     2*(m*n^3 - m^3*n)*Q66;
   Qys = m*n^3*Q11 - m^3*n*Q22 + (m^3*n - m*n^3)*Q12 + 
     2*(m^3*n - m*n^3)*Q66;
   Qoff = {{Qxx, Qxy, Qxs}, {Qxy, Qyy, Qys}, {Qxs, Qys, Qss}};
   mA = mA + Qoff*(pos[[i + 1]] - pos[[i]]);
   mB = mB + Qoff*(pos[[i + 1]]^2 - pos[[i]]^2);
   mD = mD + Qoff*(pos[[i + 1]]^3 - pos[[i]]^3);

   ];

  mB = mB/2;
  mD = mD/3;

  A11[x] = mA[[1, 1]]; A12[x] = mA[[1, 2]]; A16[x] = mA[[1, 3]];
  A22[x] = mA[[2, 2]]; A26[x] = mA[[2, 3]]; A66[x] = mA[[3, 3]];
  D11[x] = mD[[1, 1]]; D12[x] = mD[[1, 2]]; D16[x] = mD[[1, 3]];
  D22[x] = mD[[2, 2]]; D26[x] = mD[[2, 3]]; D66[x] = mD[[3, 3]];

  Nx[x_, y_] = 
   A11[x] D[u[x, y], {x, 1}] + A12[x] D[v[x, y], {y, 1}] + 
    A16[x] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

  Ny[x_, y_] = 
   A12[x] D[u[x, y], {x, 1}] + A22[x] D[v[x, y], {y, 1}] + 
    A26[x] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

  Nxy[x_, y_] = 
   A16[x] D[u[x, y], {x, 1}] + A26[x] D[v[x, y], {y, 1}] + 
    A66[x] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

  PDEs = {
    D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], 
    D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x, 1}]
    };

  gammaD =
   {
    DirichletCondition[u[x, y] == -u0/2, x == a/2], 
    DirichletCondition[u[x, y] == u0/2, x == -a/2]
    };

  omega = Rectangle[{-a/2, -b/2}, {a/2, b/2}];

  mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];

  {U, V} =
   NDSolveValue[
    {
     PDEs == {0, 0},
     gammaD,
     DirichletCondition[v[x, y] == 0, y == b/2],
     DirichletCondition[v[x, y] == 0, y == -b/2]
     }
    ,
    {u, v}, {x, y} \[Element] mesh
    ];

  Nx[x_, y_] = 
   A11 D[U[x, y], {x, 1}] + A12 D[V[x, y], {y, 1}] + 
    A16 (D[U[x, y], {y, 1}] + D[V[x, y], {x, 1}]);

  Ny[x_, y_] = 
   A12 D[U[x, y], {x, 1}] + A22 D[V[x, y], {y, 1}] + 
    A26 (D[U[x, y], {y, 1}] + D[V[x, y], {x, 1}]);

  Nxy[x_, y_] = 
   A16 D[U[x, y], {x, 1}] + A26 D[V[x, y], {y, 1}] + 
    A66 (D[U[x, y], {y, 1}] + D[V[x, y], {x, 1}]);

  nL = 20; nC = 20;
  K = ConstantArray[0, {nL*nC, nL*nC}];
  M = ConstantArray[0, {nL*nC, nL*nC}];
  Do[
   Do[
    Do[
     Do[

      k = (m - 1) nC + n;
      l = (p - 1) nC + q;

      If[n == q,
       K[[k, l]] = (\[Pi]^4 b m^2 p^2)/(2 a^4)
           NIntegrate[
           D11[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(
             p \[Pi] (x + 0.5 a))/a], {x, -a/2, a/2}, 
           Method -> "Trapezoidal"]
         + (\[Pi]^4 n^2 q^2)/(2 b^3)
           NIntegrate[
           D22[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(
             p \[Pi] (x + 0.5 a))/a], {x, -a/2, a/2}, 
           Method -> "Trapezoidal"], K[[k, l]] = 0]

      , {q, 1, nC}],
     {p, 1, nL}],
    {n, 1, nC}],
   {m, 1, nL}];
  , {angle0, 0, 90, 1}]
 , {angle1, 0, 90, 1}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.