3
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ Your array is huge. Have you tried to evaluate the integral symbolically? $\endgroup$ – AccidentalFourierTransform Feb 4 '20 at 2:56
  • $\begingroup$ @AccidentalFourierTransform, yes. It was slower than NIntegrate. Further I need evaluate 9 arrays like that. Would you know another approach? $\endgroup$ – Professor P. Cosmo Klunk Feb 4 '20 at 3:01
  • 2
    $\begingroup$ Why should I avoid the For loop in Mathematica? $\endgroup$ – Michael E2 Feb 4 '20 at 3:52
  • 1
    $\begingroup$ 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}] $\endgroup$ – Michael E2 Feb 4 '20 at 3:59
9
$\begingroup$

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

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

Your Answer

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

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