3
$\begingroup$

I have a function acf2[x,y] that describes the auto-correlation function of a simulated real surface, with a peak at the origin (acf2[0,0] = 1), and need to find the maximum distance from the origin that the function value first drops to acf2[x,y] = 1/5. A simple example of the function is

$$acf2(x,y)=\frac{1}{4} \left(e^{-31250 i \pi x}+e^{31250 i \pi x}+e^{-62500 i \pi y}+e^{62500 i \pi y}\right)$$

for a surface of size $-32\times 10^{-6} \leq x,y \leq 32\times 10^{-6}$.

Using the Nmaximize function, in the form

NMaximize[{Sqrt[x^2 + y^2], Re[acf2[x, y]] >= 1/5 && lowerbound <= x <= upperbound && lowerbound <= y <= upperbound}, {x, y}]

where upperbound and lowerbound define the boundaries of the surface, results in mathematica obtaining a position for which acf2[x,y] >= 1/5 near the boundaries of the surface. This is the expected result, but not what I am trying to find. Instead, I need to obtain the furthest distance the function first drops to 1/5. In other words, the furthest distance on the first countour line of acf[x,y] = 1/5 that encloses the origin.

The way I have attempted to solve this is by saying "all the values of acf2[x,y] on the straight line between the maximum position and the origin must be >= 1/5". Attempting to implement this I have converted the function to cylindrical coordinates

acf2polar[r_, theta_] = acf2[x, y] /. {x -> r*Cos[theta], y -> r*Sin[theta]};

and added my condition as a constraint:

NMaximize[{r, Re[acf2polar[Interval[{0, r}], theta]] >= 1/5 }, {r, theta}]

using Interval[{0,r}] to try and specify I need to constraint to hold for the whole line between r and the origin.

However, running this returns an error that the constraints aren't valid and that 'constraints should be equalities, inequalities, or domain specifications involving the variables'.

Am I on the right lines here, or is there another method I should be using?

Thanks in advance.

$\endgroup$
  • $\begingroup$ It might be worth doing this in some sense "by eye". Consider this: ContourPlot[ Re[acf2[x, y]], {x, -32 10^-6, 32 10^-6}, {y, -32 10^-6, 32 10^-6} , ContourShading -> False, Contours -> {1/5}]. This generates a contour plot where the contours are the curves where Re[acf2[x, y]] == 1/5. It looks to me like you have to move in the x direction to move toward the farthest point at which the function drops to 1/5. $\endgroup$ – march Sep 22 '17 at 16:32
  • $\begingroup$ Hi march, thanks for the reply. I've already had a look at using ContourPlot, and unfortunately trying to extract [x,y] values from the plotted points doesn't give very accurate results, due to the limited number of points plotted on the graph. The application of the problem requires a fairly high degree of accuracy, and so I think a mathematical approach is required. $\endgroup$ – L.Todhunter Sep 25 '17 at 14:27
1
$\begingroup$

One possible solution is to find the maximal radius for each angle, and then optimize over the angle. I also repeat your definitions for completeness.

acf2[x_, y_] := 1/4 (E^(-31250 I π x) + E^(31250 I π x) + 
                     E^(-62500 I π y) + E^(62500 I π y));
acf2polar[r_, theta_] := acf2[x, y] /. {x -> r*Cos[theta], y -> r*Sin[theta]};
InThatDirection[theta_?NumericQ] := s /. FindRoot[Re[acf2polar[s, theta] - 1/5], {s, 0}, 
   Method -> {"Newton", "StepControl" -> {"TrustRegion",
              "StartingScaledStepSize" -> 0.001}}];

The function InThatDirection finds the first root of acf2polar[s,theta]-1/5 for a given value of theta. Note that I had to change the trust region, because otherwise FindRoot would skip over the first root (so keep this in mind when you apply this approach to other correlation functions).
Now we can maximize the radius over all angles:

NMaximize[{InThatDirection[theta], 0 < theta < 2 π}, theta, Method -> "RandomSearch"]

{0.0000225546, {theta -> 1.64376*10^-9}}

Or if you want to get the point:

p = #[[1]] {Cos[theta], Sin[theta]} /. #[[2]] &@%

{0.0000225546, 3.70743*10^-14}

$\endgroup$
  • $\begingroup$ This solves the problem. Thank you! $\endgroup$ – L.Todhunter Oct 27 '17 at 10:23
  • $\begingroup$ Good to hear it worked for you! You can also accept the answer to express that. $\endgroup$ – M. Stern Oct 27 '17 at 14:53

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.