1
$\begingroup$

I need to DensityPlot an NIntegrate where near the maximum piont. So I need to use FindMaximumto find the my plot range, but it has been running for more than 48 hours and has not given any results if I don't define starting points and intervals, if I define the starting point and the interval, then the time will decrease as the interval decreases. So I think my integral function needs to speed up, I used Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0} method, it does improve the speed, but when I run to find the maximum, it obviously gives a wrong answer and some errors, so I use "LocalAdaptive", But it made the program run for a few hours too.

g[xo,zo] is the function that I want the maximum value and plot the density diagram, Its physical meaning is the wave field intensity of light. It's defined by the following functions.

f1[xo_?NumericQ, zo_?NumericQ] := 
  NIntegrate[(E^(I k zo)  E^(((I k) xo^2)/(2 zo))
        YO[x] E^(((I k) x^2)/(2 zo)) E^(-(((I k) x xo)/zo)))/
    Sqrt[I λ zo], {x, -30000, 30000}, 
   Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
   MinRecursion -> 3, WorkingPrecision -> 10];

g[xo_?NumericQ, zo_?NumericQ] := Abs[f1[xo, zo]^2];

FindMaximum[{g[xo, zo], -300. <= xo <= 300. && 4000000. <= zo <= 5000000}, {{xo, 0}, {zo, 4720000}}]

The above program in the case that I know the maximum interval with Findmaximum can calculate the result in the time I can receive. It can output

{1634.28, {xo -> 2.19213*10^-10, zo -> 4.71831*10^6}}  

It spend 3563.0903136 seconds. And DensityPlot takes Longer Longer time.

DensityPlot[g[xo, zo], {zo, 4698310, 4738310}, {xo, -200, 200}, 
 PlotLegends -> Automatic, ColorFunction -> "Rainbow", Mesh -> None, 
 MaxRecursion -> 5, PlotPoints -> 30, PlotRange -> All]

But if I don't know the FindMaximum interval, I need to:

FindMaximum[{g[xo, zo], -300. <= xo <= 300. && 0. <= zo <= 5000000}, {{xo, 0}, {zo, 4720000}}]

And it took me dozens of hours and it kept running so I stopped the calculation.

So First:
I want to know if there is any NIntegrate method that allows me to find the maximum value accurately and within the time not too long, even if it's a little bit less accurate.

Second:
If there are some ways to make DensityPlot faster? I tried some other ways to integrate and Plot quickly. It took me some times, but I failed because I was not very good at mathematics.Like Henrik Schumacher's answer

The designed focus is zo= 4.72 mm = 4.72*10^6 nm, so the focus should be around (xo=0,zo=4720000) which makes me know the interval of FindMaximum .

Ps: what I did is the content of a document and I want to repeat it.the document link

The code definition of YO[x] is as follows, it contains a solution of the partial differential equation, and it is proved that the solution is completely correct.

Clear["Global`*"]
t = AbsoluteTime[];
θ = 1.1/1000;
λ = 1.2398/19.5 ;
f = 4.72 10^6;
Subscript[δ, 1] = 1.274/10^6;
Subscript[δ, 2] = 4.304/10^6;
Subscript[bt, 1] = 5.254/10^9;
Subscript[bt, 2] = 2.435/10^7;
χ1 = -2 Subscript[δ, 1] + 2 I Subscript[bt, 1];
χ2 = -2 Subscript[δ, 2] + 2 I Subscript[bt, 2];
Δχ = χ1 - χ2;
k = 2 (π/λ);
Table[β[b] = -((2 (b x) Sin[θ])/f) - ((b x)/f)^2, {b, -5,5}];
Table[χ[a] = (Δχ (1-(-1)^Abs[a]))/(2 I a π), {a, -10, -1}];
Table[χ[a] = (Δχ (1-(-1)^Abs[a]))/(2 I a π), {a, 1, 10}];
Table[χ[a] = (χ1 + χ2)/2, {a, 0, 0}];
eqns = 
  Join[
    Table[
      (Sin[θ] + (h x)/f) D[Subscript[Y, h][x, z], x] + 
       Cos[θ] D[Subscript[Y, h][x, z], z] == 
      ((I Pi) 
        (β[h] Subscript[Y, h][x, z] + 
          Sum[χ[h - l] Subscript[Y, l][x, z], {l, -5, 5}]))/λ,{h, -5, 5}],
    Table[Subscript[Y, h][x, 0] == If[h == 0, 1, 0], {h, -5, 5}]];
s = NDSolve[eqns, 
   Table[Subscript[Y, h], {h, -5, 5}], {x, 0, 30000}, {z, 0, 30000}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MaxPoints" -> 100, "MinPoints" -> 100, 
       "DifferenceOrder" -> 4}, "TemporalVariable" -> z}];

Clear[YO]
YO[x_?NumericQ /; (0 > x >= -30000)] := YO[-x]
YO[x_?NumericQ] = 
  First[Sum[E^((I h π)/(λ f) x^2) Subscript[Y, h][x, 13500], {h, -5, -1}] /. s];

This is the picture should:
enter image description here

This is what the picture I have make, and I just calculate h=-1 of YO[x] for calculating {h,-5,-1} is too long to calculate. enter image description here

$\endgroup$
2
  • $\begingroup$ To @shrocat Haven't you seen the error message of NDSolve for defining the YO: NDSolve::bcart: Warning: An insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution. >> . Therefore it makes no sense to calculate further. Give additional boundary conditions. $\endgroup$
    – Akku14
    Commented Jan 15, 2021 at 6:17
  • $\begingroup$ @Akku14 It sure have the error message, but I think it is fine because using YO[x] I can calculate the same result as the literature, The problem is that's too slow, If I do difference integration useing other languages like Matlab to Plot, it reduced accuracy but increased speed. So I wonder if I can adjust the parameter of the Method of NIntegrate to reduce the accuracy within allowable limits for speeding up calculation. $\endgroup$
    – shrocat
    Commented Jan 15, 2021 at 10:27

0

Your Answer

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