5
$\begingroup$

I want to simulate the fresnel diffraction through circular aperture.But mathematica is taking forever to calculate.Here is my code.

  f[x_?NumericQ, y_?NumericQ, d_?NumericQ, p_?NumericQ, t_?NumericQ] := 
 NIntegrate[(\[Zeta]*(Exp[
       I*20000 ((d^2 + \[Zeta]^2)^(0.5) + ((p^2 + (x - \[Zeta]*
                  Cos[t])^2 + (y - \[Zeta]*
                  Sin[t])^2)^(0.5)))]))/(d^2 + \[Zeta]^2)^(0.5), {\
\[Zeta], 0, 0.008}]
w[x_?NumericQ, y_?NumericQ, d_?NumericQ, p_?NumericQ] := 
 NIntegrate[f[x, y, d, p, t], {t, 0, 2*Pi}]
g[x_, y_, d_, p_] := (Abs[w[x, y, d, p]])^2
Plot3D[g[x, y, 0.004, 20], {x, -100, 100}, {y, -100, 100}, 
 PlotRange -> Full]

What's the way out? *Specs:i5 7th gen intel processor,12 gb ram

$\endgroup$
9
$\begingroup$

Let's compile the integrand into a Listable CompiledFunction:

cintegrand = Block[{x, y, d, p, ζ, t},
   With[{code = 
      N[(ζ (Exp[I 20000 ((d^2 + ζ^2)^(1/2) + ((p^2 + (x - ζ Cos[t])^2 + (y - ζ Sin[t])^2)^(1/2)))]))/(d^2 + ζ^2)^(1/2)]
     },
    Compile[{{x, _Real}, {y, _Real}, {d, _Real}, {p, _Real}, {ζ, _Real}, {t, _Real}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True
     ]
    ]
   ];

Next, pick a Gauss quadrature rule for high order polynomials (the integrand is very smooth but quite oscillatory):

{pts, weights, errweights} = NIntegrate`GaussRuleData[9, $MachinePrecision];

Divide the integration domain into $m \times n$ rectangles and set up the quadrature points and weights:

m = 20;
n = 20;
ζdata = Partition[Subdivide[0., 0.008, m], 2, 1];
tdata = Partition[Subdivide[0., 2 Pi, n], 2, 1];
{ζ, t} = Transpose[Tuples[{Flatten[ζdata.{1. - pts, pts}], Flatten[tdata.{1. - pts, pts}]}]];
ω = Flatten[KroneckerProduct[
    Flatten[KroneckerProduct[Differences /@ ζdata, weights]],
    Flatten[KroneckerProduct[Differences /@ tdata, weights]]
    ]];

Define the mapping that parameterizes the graph of OP's function w:

W = {x, y} \[Function] {x, y, Abs[cintegrand[x, y, 0.004, 20., ζ, t].ω]^2};

Compute the values of W on a $101 \times 101$ grid:

R = 2.;
data = Outer[W, Subdivide[-R, R, 100], Subdivide[-R, R, 100]]; // AbsoluteTiming // First

28.2668

Plot the result

ListPlot3D[Flatten[data, 1], PlotRange -> All, AxesLabel -> {"x", "y", "w"}]

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Could you please recomend a source for learning such techniques for reducing time of computation,like the use of block,listablecompiledfunction etc..? $\endgroup$ – user157588 May 24 '18 at 13:26
  • $\begingroup$ This is a good start. The rest is basically experience, some basics in numerical Mathematics and computer science, reading a lot on StackExchange,... and experience of course. $\endgroup$ – Henrik Schumacher May 24 '18 at 13:32
  • $\begingroup$ Oh, and I would not call Block a way of reducing computational time. It is just a scoping construct that works a bit more predictably in conjunction with Compile than Module. In general, I would advise Module. $\endgroup$ – Henrik Schumacher May 24 '18 at 13:34
3
$\begingroup$

TIP: In Mathematica use exact values,because have infinty precision.

f[x_?NumericQ, y_?NumericQ, d_?NumericQ, p_?NumericQ] := 
NIntegrate[(ζ*(Exp[I*20000 ((d^2 + ζ^2)^(1/2) + ((p^2 + (x - ζ*Cos[t])^2 + 
(y - ζ*Sin[t])^2)^(1/2)))]))/(d^2 + ζ^2)^(1/2), {ζ, 0, 8/1000}, {t, 0, 2*Pi}, 
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}];

g[x_, y_, d_, p_] := Abs[f[x, y, d, p]]^2;

n = 1/6; (*The smaller the value , 3D plot is smoother,but CPU time will increase !!! *)

ListPlot3D[Partition[Flatten[Table[{x, y, g[x, y, 4/1000, 20]},
{x, -2, 2, n}, {y, -2, 2, n}]], 3]]
(*In Range: -2<x<2 and -2<y<2.If you increase range domain CPU time will increase !!! *)

enter image description here

On my a 100$old laptop calculation time is about 24 minutes.

|improve this answer|||||
$\endgroup$

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.