Efficient way to evaluate numerically product sines functions

Question

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 04 '20 at 02:56
@AccidentalFourierTransform, yes. It was slower than NIntegrate. Further I need evaluate 9 arrays like that. Would you know another approach? — Riobaldo Tatarana, Feb 04 '20 at 03:01
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 04 '20 at 03:59

Henrik Schumacher · Accepted Answer · 2020-04-12T13:40:40.833

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 missing motivation, Feb 04 '20 at 08:20

Efficient way to evaluate numerically product sines functions

1 Answers1