# Efficient way to evaluate numerically product sines functions

I'm trying to perform integration of the following matrix.

contL = 1; contC = 1;
K = ConstantArray[0, {20*20, 20*20}];
For[m = 1, m <= 20, m++,
For[n = 1, n <= 20, n++,
For[p = 1, p <= 20, p++,
For[q = 1, q <= 20, q++,
K[[contL, contC]] = m^2 p^2 Sin[m x] Sin[n y] Sin[p x] Sin[q y];
contC++
];
];
contC = 1;
contL++;
];
]
NIntegrate[K, {x, -1/2, 1/2}, {y, -1/2, 1/2}]


Someone knows how to evaluate that integration faster? It is taking heavy time to evaluate. I already tried to use Method -> "Trapezoidal", but it yet doesn't improved the computation time.

Thank you in advance.

Regards.

• Your array is huge. Have you tried to evaluate the integral symbolically? – AccidentalFourierTransform Feb 4 '20 at 2:56
• @AccidentalFourierTransform, yes. It was slower than NIntegrate. Further I need evaluate 9 arrays like that. Would you know another approach? – Professor P. Cosmo Klunk Feb 4 '20 at 3:01
• Why should I avoid the For loop in Mathematica? – Michael E2 Feb 4 '20 at 3:52
• Re @AccidentalFourierTransform's suggestion: Consider Integrate[ m^2 p^2 Sin[m x] Sin[n y] Sin[p x] Sin[q y], {x, -1/2, 1/2}, {y, -1/2, 1/2}] and Integrate[ m^2 n^2 Sin[m x] Sin[n y] Sin[m x] Sin[n y], {x, -1/2, 1/2}, {y, -1/2, 1/2}] – Michael E2 Feb 4 '20 at 3:59

## 1 Answer

First of all, I observed that your matrix can be represented as KroneckerProduct of smaller matrixes which are themselves KroneckerProducts (or TensorProducts) of vectors:

Ux = Table[m^2 Sin[m x], {m, 1, k}];
Vx = Table[p^2 Sin[p x], {p, 1, k}];
Uy = Table[Sin[n y], {n, 1, k}];
Vy = Table[Sin[q y], {q, 1, k}];
Ax = KroneckerProduct[Ux, Vx];
Ay = KroneckerProduct[Uy, Vy];
T = KroneckerProduct[Ax, Ay];

K == T


True

The integrations in x and y are independent of each other, so we could compute Ax and Ay first, requiring $$2 k^2$$ NIntegrates instead of $$k^4$$ in your original code. However, we can actually do the integrations symbolically first and numericize the resulting expressions. If you want that really fast, you can create a CompiledFunction for the numericization; this will also guarantee that Ax and Ay are packed arrays so that KroneckerProduct[Ax, Ay] evaluates quickly.

Creating the compiled function:

Block[{a, b, m, n},
With[{
mmcode = N@Integrate[Sin[m x] Sin[m x], {x, a, b}],
mncode = N@Integrate[Sin[m x] Sin[n x], {x, a, b}]
},
cf = Compile[{{m, _Real}, {n, _Real}, {a, _Real}, {b, _Real}},
If[m == n,
If[m == 0, 0., mmcode],
If[m == 0 || n == 0, 0., mncode]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];
]
];


Computing Ay, Ax, and T:

First@AbsoluteTiming[
Ay = Outer[cf[#1, #2, -1/2, 1/2] &, Range[k], Range[k], 1];
Ax = Ay KroneckerProduct[#, #] &[Range[1., k]^2];
T = KroneckerProduct[Ax, Ay];
]


0.001463

I did not wait for the result for integrals over K in case of k = 20, but for k = 5, it took about 3.8639 seconds on my machine, while my method took only 0.000056 seconds, a speed-up of almost $$69\,000$$. Also, T should be precise up to machine precision as there is no numerical integration employed at all (so all errors stem only from the very few floating point operations +, -, *, from the approximations of Sin and Cos and from division by integers (have a look at mncode).

• Kronecker products sure are neat; a lot of things simplify greatly as soon as you recognize the appropriate structure applies. – J. M.'s torpor Feb 4 '20 at 8:20
• Perfect, @Henrik! Thank you so much! – Professor P. Cosmo Klunk Feb 4 '20 at 22:40
• You're welcome! =) – Henrik Schumacher Feb 4 '20 at 23:00