# How to improve the performance of this code

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. Dec 1, 2013 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. Dec 2, 2013 at 2:27

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

• Yes, Ei is certainly a function!! I will try your suggestions
– Alex
Dec 1, 2013 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... Dec 2, 2013 at 0:04
• r = NIntegrate[v, {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}]*d is not working. Mathematica complains of singularities
– Alex
Dec 2, 2013 at 0:05
• your vector quantities (P0..) fail the NumericQ test, just use NumericQ on x Dec 2, 2013 at 12:54
• worth reading re: the "SymbolicEvaluation" option mathematica.stackexchange.com/questions/9971/… Dec 2, 2013 at 15:51