3
$\begingroup$

I'm new here. Actually, Mathematica is so cool, I can just write about 5 lines to do the same thing comparing to over 100 lines in Python. But the speed is very slow. I wonder if there is any method that can improve the simulation speed? Like using different algorithm, do not use mathematica or do not use the function it already has, etc.

Here is my code. I was about to set the width of each step to 0.001. But I cannot get the result for a long while. So I gave up and change it to 0.1. It also takes about an hour.

λ = 1000;
z = 100;
Array[f, {100, 201, 201}];
f = Table[
   1/2*Abs[1/(I*λ)*
       NIntegrate[
        If[x0^2 + y0^2 < 25, 
         Exp[I*2*Pi/λ*Sqrt[(x - x0)^2 + (y - y0)^2 + z^2]]/
          Sqrt[(x - x0)^2 + (y - y0)^2 + z^2], 0], {x0, -10, 
         10}, {y0, -10, 10}]]^2, {x, -10, 10, 0.1}, {y, -10, 10, 0.1}];
data = Table[
   f[[x - 100.1, y - 100.1]], {x, -10, 10, 0.1}, {y, -10, 10, 0.1}];
ListPlot3D[data, Mesh -> None, InterpolationOrder -> 3, 
 ColorFunction -> "SouthwestColors"]

Thanks for your help!

$\endgroup$
  • 4
    $\begingroup$ With step size 0.001, you try to generate a matrix of size $20000 \times 20000$. That means that 400 million integrals have to be computed. On my machine, a single integral takes 6 ms, so the full computation would need 666 hours, hence almost a month. But your problem has a rotational symmetry: The value of the integral depends only on the length of the vector {x,y}. So in principle, it suffices to compute the integrals only for points of the form {x,0} and to rotate the result about the origin. That should save you a lot of time. $\endgroup$ – Henrik Schumacher Jul 12 '18 at 5:46
8
$\begingroup$

This is about seven times as fast as your code on my computer:

g[x0_, y0_, x_, y_, z_, lambda_] := With[{sqrt = Sqrt[(x - x0)^2 + (y - y0)^2 + z^2]}, Exp[2 I Pi sqrt/lambda]/sqrt]
h[x_, y_, z_, lambda_] := 0.5 Abs[1/(I \[Lambda]) NIntegrate[g[x0, y0, x, y, z, lambda], {x0, y0} \[Element] Disk[{0, 0}, 5]]]^2

\[Lambda] = 1000;
z = 100;
Table[
   h[x, y, z, \[Lambda]],
   {x, -10, 10, 0.5},
   {y, -10, 10, 0.5}
   ]; // AbsoluteTiming

It looks different because I rewrote to more clearly see what I was optimizing, but it should be the same. More or less all the gain in speed is from getting rid of the If and instead specifying that NIntegrate should only integrate over the region where the function is nonzero. This is that part:

{x0, y0} \[Element] Disk[{0, 0}, 5]

You were integrating over a square and then left it up to NIntegrate to handle the discontinuous change where the function goes to zero. This way is faster.

I also did a micro-optimization in that I used With to compute a quantity that was previously computed twice only once, it is a minor thing but it also helps with the clarity of the code.

$\endgroup$
  • 1
    $\begingroup$ Also, when your machine has multiple cores you could use ParallelTable instead of Table as a simple way to cut calculation time. $\endgroup$ – Ruud3.1415 Jul 12 '18 at 16:32
  • $\begingroup$ Thanks! But there is a small issue occurred: I cannot get the right number in x and y axis if I just simply use ListPlot3D to plot it(the same issue on my own code). $\endgroup$ – FeverDream Jul 13 '18 at 0:43
  • $\begingroup$ oh my bad, I know the answer. It is just because I forgot to put x and y into the table $\endgroup$ – FeverDream Jul 13 '18 at 1:05

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.