NIntegrate: how to speed up code?

Question

a = 4;(* nodes  Х  (-a,a)*)
b = 4;(* nodes У  (-b,b)*)
n = (2 a + 1) ( 2 b + 1);(*all nodes *)
x0 = 2;
y0 = 2;
k = 1/2 // N;
Ax = k/(x0)^2;
Ay = k/(y0)^2;
EE = 2 10^5;
μ = 0.3;
h = 0.8;
Dc = EE*h^3/12/(1 - μ^2) ;
ϕ = {Flatten[Table[Exp[-Ax (x - xi)^2] Exp[-Ay (y - yi)^2], {xi, -a, 
                                a}, {yi, -b, b} ]]};
ϕT = Transpose[ϕ];
solK = Dc/2 Laplacian[(ϕT . ϕ),{x,y}];
K1 = NIntegrate[solK, {x, -a, a}, {y, -b, b}, 
     Method -> {"MultidimensionalRule", "Generators" -> 5, 
                "SymbolicProcessing" -> 0}, PrecisionGoal -> 3, AccuracyGoal -> 3];

Hi! I want find deformation (strain) energy for plane by pressure .In one of the steps I need to use NIntegrate function so problem with speed of my code when I take (a=6, b=6) it works too slowly! What can I do to inprove code?

Answer 1

The code in the Question runs slowly, because it evaluates n^2 integrals. However, all the integrands are of the form

Dc Exp[-Ax (x - a1)^2] Exp[-Ay (y - b1)^2] Exp[-Ax (x - a2)^2] Exp[-Ay (y - b2)^2];

This generic term can be integrated symbolically in several seconds

Integrate[Laplacian[%, {x, y}], {x, -a, a}, {y, -b, b}]

(* Dc (E^((-2*a^2 - a1^2 - a2^2 - 4*b^2 - b1^2 - b2^2)/8)*Sqrt[Pi]*
  (-((2*a + a1 + a2)*E^((-4*a*(a1 + a2) + 8*b^2 + (b1 + b2)^2)/16)*Erf[(2*b - b1 - b2)/4]) - 
   (2*a - a1 - a2)*E^((4*a*(a1 + a2) + 8*b^2 + (b1 + b2)^2)/16)*Erf[(2*b - b1 - b2)/4] + 
   E^((b*(b - b1 - b2))/4)*(E^((4*a^2 + (a1 + a2)^2)/16)*(-((b1 + b2)*(-1 + E^((b*(b1 + b2))/2))) + 2*b*(1 + E^((b*(b1 + b2))/2)))*
      Erf[(-2*a - a1 - a2)/4] - (2*a + a1 + a2)*E^((-4*a*(a1 + a2) + (2*b + b1 + b2)^2)/16)*Erf[(2*b + b1 + b2)/4]) - 
   E^((b*(b - b1 - b2))/4)*(E^((4*a^2 + (a1 + a2)^2)/16)*(-((b1 + b2)*(-1 + E^((b*(b1 + b2))/2))) + 2*b*(1 + E^((b*(b1 + b2))/2)))*
      Erf[(2*a - a1 - a2)/4] + (2*a - a1 - a2)*E^((4*a*(a1 + a2) + (2*b + b1 + b2)^2)/16)*Erf[(2*b + b1 + b2)/4])))/4 *)

which can be evaluated by substitution to obtain the desired numerical values. Note that k has been set to 1/2, rather than 0.5 as in the Question.

Sidelight

Defining

ϕx = Table[Exp[-Ax (x - xi)^2], {xi, -a, a}];
ϕy = Table[Exp[-Ay (y - yi)^2], {yi, -b, b}];

ϕT . ϕ is given by

Outer[Times, Flatten[Outer[Times, ϕx, ϕy]], Flatten[Outer[Times, ϕx, ϕy]]];

Further, if Flatten is omitted, which affect only the order of array elements, this last expression becomes

Outer[Times, ϕx, ϕy, ϕx, ϕy];

or, equivalently (up to a reordering of array elements),

Outer[Times, ϕx, ϕx, ϕy, ϕy];

Answer 2

A quick and dirty use of Parallelize gave me a 2.7x speed improvement on my machine (4-core CPU w/ hyperthreading). But since NIntegrate itself cant be parallelized, I used Map to do the trick

K1=Parallelize[NIntegrate[#,{x,-a,a},{y,-b,b},Method->{"MultidimensionalRule","Generators"->5,"SymbolicProcessing"->0},PrecisionGoal->3,AccuracyGoal->3]&/@solK]

hope this helps :)