0
$\begingroup$

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];
Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ 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. $\endgroup$
    – bill s
    Commented Dec 1, 2013 at 22:58
  • $\begingroup$ 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. $\endgroup$
    – Jagra
    Commented Dec 2, 2013 at 2:27

1 Answer 1

Reset to default
1
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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