Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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];
share|improve this question
2  
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 1 '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 2 '13 at 2:27
add comment

1 Answer

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

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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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