# How to evaluate that integral accurately and faster?

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?

• Did you mean $n\pi$ rather than $p\pi$?
– JimB
Feb 12, 2020 at 23:19
• Hello, @JimB! Oh, I saw, I edited. Feb 12, 2020 at 23:30
• 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. Feb 13, 2020 at 0:02
• By the way, capital $K$ is a reserved word, so using it as a variable name could cause problems. Feb 13, 2020 at 0:06
• 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. Feb 13, 2020 at 1:11

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.

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).

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;


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

Needs["NDSolveFEM"];
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}]
};

{
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},
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}]