How to improve the performance of this code

Question

I have implemented a code in Mathematica 9 to simulate a scattering problem, and I got really disappointed about its performance when compared to Matlab. Since I am a newbie in Mathematica, would someone give me clues here?

ClearAll["Global`*"];
integrand[k_, P0_, P1_, rho_, l_] = Module[{x, h},
   x = P1*l + P0*(1. - l) - rho;
   h = HankelH2[0, Norm[x, 2]*k];
   h
   ];
innerInt[k_, P0_, P1_, rho_] := Module[{d, r, v, x},
   d = Norm[P0 - P1, 2];(* distancia entre os pontos *)
   v = integrand[k, P0, P1, rho, x];
   r = NIntegrate[v, {x, 0, 1}]*d;
   r
   ];
solveSystemCylinder[segsStart_, segsEnd_, Ei_, k_, eta_] := 
  Module[{colocPts, r, x, t, eiVals},
   colocPts = (segsStart + segsEnd)/2;
   t[x_] := MapThread[innerInt[k, #1, #2, x] &, {segsStart, segsEnd}];
   r = Map[t, colocPts];
   eiVals = Map[Ei[#, k] &, colocPts];
   N[LinearSolve[r, eiVals]]
   ];
Ei[p_, k_] := Exp[I*k*First[p]];

k := 2*Pi/Lambda;
Lambda := 1;
raio := Lambda;
eta := 1;
ND := 30;(* num de divisões *)
pts := N[raio*{Cos[#], Sin[#]} & /@ (Range[0, ND]*(2*Pi/ND))];
pts2 := ({0, 2*Lambda} + #) & /@ pts;
allPts := Join[pts, pts2];

startPts := Take[pts, Length[pts] - 1];
endPts := Take[pts, -(Length[pts] - 1)];

startPts2 := Take[pts2, Length[pts2] - 1];
endPts2 := Take[pts2, -(Length[pts2] - 1)];

(*resolve*)
Js = solveSystemCylinder[Join[startPts, startPts2], 
   Join[endPts, endPts2], Ei, k, eta];

One thing that might be slowing you down is the erratic uses of SetDelayed and Set. At the very least, your constants should not be defined with := and your functions should be defined with. You can use Timing to tell how long various function calls take -- figuring out which parts are slow will go a ways towards understanding what can be optimized. — bill s, Dec 01 '13 at 22:58
As Alex has defined Ei as a function, does he really need to pass it to solveSystemCylinder[]? It might make more sense stylistically and clarity wise to make Ei[p_,k_] a local function of solveSystemCylinder[] similar to t[x_]. Doing so doesn't appear to affect performance, but sometimes these kind of things can clarify one's approach to a problem. — Jagra, Dec 02 '13 at 02:27

Cassini · Answer 1 · 2013-12-02T00:09:21.593

1

I've made a number of changes to your code which may speed things up, but I honestly can't say with certainty. Rather than enumerating the changes, I'll just list the code here:

ClearAll["Global`*"];
integrand[k_?NumericQ, P0_?NumericQ, P1_?NumericQ, rho_?NumericQ, 
   l_?NumericQ] :=
  Module[{x, h},
   x = P1*l + P0*(1. - l) - rho;
   h = HankelH2[0, Norm[x]*k]
   ];
innerInt[k_?NumericQ, P0_?NumericQ, P1_?NumericQ, rho_?NumericQ] :=
  Module[{d, r},
   d = Norm[P0 - P1];(*distancia entre os pontos*)
   r = NIntegrate[integrand[k, P0, P1, rho, x], {x, 0, 1}, 
      Method -> {Automatic, "SymbolicProcessing" -> 0}]*d
   ];
solveSystemCylinder[segsStart_, segsEnd_, Ei_, k_, eta_] := 
  Module[{colocPts, r, x, t, eiVals},
   colocPts = (segsStart + segsEnd)/2;
   t[x_] := 
    MapThread[innerInt[k, #1, #2, x] &, {segsStart, segsEnd}];
   r = Map[t, colocPts];
   eiVals = Map[Ei[#, k] &, colocPts];
   N[LinearSolve[r, eiVals]]
   ];
Ei[p_, k_] := Exp[I*k*First[p]];
k = 2*Pi/Lambda;
Lambda = 1;
raio = Lambda;
eta = 1;
ND = 30;(*num de divisões*)
pts := N[raio*{Cos[#], Sin[#]} & /@ (Range[0, ND]*(2*Pi/ND))];
pts2 := ({0, 2*Lambda} + #) & /@ pts;
allPts := Join[pts, pts2];
startPts := Take[pts, Length[pts] - 1];
endPts := Take[pts, -(Length[pts] - 1)];
startPts2 := Take[pts2, Length[pts2] - 1];
endPts2 := Take[pts2, -(Length[pts2] - 1)];
(*resolve*)
Js = solveSystemCylinder[Join[startPts, startPts2], 
   Join[endPts, endPts2], Ei, k, eta];

Some of the changes were cosmetic, mostly because they made it easier for me to read.

Hope this helps...

edited Dec 02 '13 at 00:09

answered Dec 01 '13 at 23:41

Cassini

5,556
5
28
38

Yes, Ei is certainly a function!! I will try your suggestions – Alex Dec 01 '13 at 23:46
I also forced innerInt to only accept numeric quantities and I set it running, but I have no idea what a "reasonable" amount of execution time is. It's going now... – Cassini Dec 02 '13 at 00:04
r = NIntegrate[v, {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}]*d is not working. Mathematica complains of singularities – Alex Dec 02 '13 at 00:05
forcing integrand to accept only numerical values does not work... Mathematica complains it is receiving non numerical values. – Alex Dec 02 '13 at 00:07
1

your vector quantities (P0..) fail the NumericQ test, just use NumericQ on x – george2079 Dec 02 '13 at 12:54
for completeness if you want to do the numeric test on the vectors use : PO_List?(And@@NumericQ/@#&). I dont think it will help performance here though – george2079 Dec 02 '13 at 15:21
1

worth reading re: the "SymbolicEvaluation" option http://mathematica.stackexchange.com/questions/9971/the-difference-between-symbolicprocessing-0-and-restricting-the-function-de – george2079 Dec 02 '13 at 15:51

How to improve the performance of this code

1 Answers1