1
$\begingroup$ 
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?

| improve this question | |
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Apr 23 '15 at 17:29
  • $\begingroup$ Independent of speed, are you getting the answer your expect? You are computing a very large array of integrals. $\endgroup$ – bbgodfrey Apr 23 '15 at 17:41
  • $\begingroup$ yes,i get answer that is need $\endgroup$ – Arthuros Apr 23 '15 at 17:44
  • $\begingroup$ yes, with a very large array of integrals. I just want to make it faster.I know that Mathematica have function Compile can i do smth without it function (Because I'm not good at Parallel Computing) to make code faster or this is better way to use Compile ? $\endgroup$ – Arthuros Apr 23 '15 at 18:37
  • $\begingroup$ note that solK is symmetric, so you shouldn't integrate redundant terms. $\endgroup$ – george2079 Apr 23 '15 at 20:38
2
$\begingroup$

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];
| improve this answer | |
$\endgroup$
  • $\begingroup$ Dear bbgodfrey,thank you very much for useful ideas! $\endgroup$ – Arthuros Apr 23 '15 at 20:57
0
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thak you!Good idea! $\endgroup$ – Arthuros Apr 23 '15 at 21:24

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.