✕ SortBy[Tally[Round[Differences[zeros2Terms], 0.0001]], Last][[-1]]

Let us look at and next to each other. One observes an obvious difference between the two curves: in , ones sees triples of maximum-minimum-maximum, with all three extrema being negative. This situation does not occur in .

✕ GraphicsRow[{Plot[fGolden[x], {x, 0, 14 Pi}, PlotLabel -> Subscript[f, \[Phi]], ImageSize -> 320], Plot[f2Terms[x], {x, 0, 14 Pi}, PlotLabel -> Subscript[f, 2], ImageSize -> 320]}]

Here is a plot of after a zero for 1,000 randomly selected zeros. Graphically, one sees many curves crossing zero exactly at the first zero of . Mouse over the curves to see the underlying curve in red to better follow its graph.

✕ Show[Plot[f2Terms[t - #], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ RandomSample[zeros2Terms, 1000], PlotRange -> All, Frame -> True, GridLines -> {{0}, {3}}, ImageSize -> 400] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]

The distance from the origin to the first zero of is indeed the most commonly occurring distance between zeros.

✕ {#, N[#]} &@(2 x /. Solve[f2Terms[x] == 0 \[And] 1 < x < 2, x] // Simplify)

Let us note that adding more terms does not lead to more zeros. Consider the function .

✕ f2TermsB[x_] := Cos[GoldenRatio^2 x]

✕ zeros2TermsB = findZeros[f2TermsB, 10^5];

It has many more zeros up to a given (large) than the function .

✕ Max[zeros2TermsB]/Max[zerosGolden]

The reason for this is that the addition of the slower oscillating functions ( and ) effectively “converts” some zeros into minimum-maximum-minimum triples all below (above) the real axis. The following graphic visualizes this effect.

✕ Plot[{fGolden[x], Cos[GoldenRatio^2 x]}, {x, 10000, 10036}, Filling -> Axis]

Seeing Envelopes in Families of Curves

Now let’s come back to our function composed of three cos terms.

If we are at a given zero of , what will happen afterward? If we are at a zero that has distance to its next zero , how different can the function behave in the interval ? For the function , there is a very strong correlation between the function value and the zero distances, even if we move to the function values in the intervals , .

✕ Table[Histogram3D[{#[[1, 1]], fGolden[#[[-1, 2]] + #[[-1, 1]]/2]} & /@ Partition[Transpose[{differencesGolden, Most[zerosGolden]}], k, 1], 100, PlotLabel -> Row[{"shift:\[ThinSpace]", k - 0.5}], AxesLabel -> {HoldForm[d], Subscript[f, \[Phi]], None}], {k, 4}] // Partition[#, 2] &

Plotting the function starting at zeros shows envelopes. Using a mouseover effect allows one to highlight individual curves. The graphic shows a few clearly recognizable envelopes. Near to them, many of the curves cluster. The four purple points indicate the intersections of the envelopes with the axis.

✕ Show[Plot[fGolden[# + t], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ RandomSample[zerosGolden, 1000], Epilog -> {Directive[Purple, PointSize[0.012]], Point[{#, 0} & /@ {1.81362, 1.04398, 1.49906, 3.01144}]}, PlotRange -> All, Frame -> True, GridLines -> {{0}, {-1, 3}}, ImageSize -> 400] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]

The envelope that belongs to zeros that are followed with nearly reaching explains the position of the largest maximum in the zero-distance distribution.

✕ zerosGolden[[1]]

Using transcendental roots, one can obtain a closed-form representation of the root that can be numericalized to any precision.

✕ Solve[fGolden[x] == 0 \[And] 1/2 < x < 3/2, x]

✕ root3 = N[ Root[{Cos[#1] + Cos[GoldenRatio #1] + Cos[GoldenRatio^2 #1] &, 0.9068093955855631129`20}], 50]

The next graphic shows the histogram of the zero distances together with the just-calculated root.

✕ Histogram[differencesGolden, 1000, ChartStyle -> Opacity[0.4], PlotRange -> {{1.7, 1.9}, All}, GridLines -> {{2 zerosGolden[[1]]}, {}}, GridLinesStyle -> Purple]

We select zeros with a spacing to the next zero that are in a small neighborhood of the just-calculated zero spacing.

✕ getNearbyZeros[zeros_, z0_, \[Delta]_] := Last /@ Select[Transpose[{Differences[zeros], Most[zeros]}], z0 - \[Delta] < #[[1]] < z0 + \[Delta] &]

✕ zerosAroundPeak = getNearbyZeros[zerosGolden, 2 root3, 0.001]; Length[zerosAroundPeak]

Plotting these curves shifted by the corresponding zeros such that they all have the point in common shows that all these curves are indeed locally (nearly) identical.

✕ Show[Plot[fGolden[# + t], {t, -3, 6}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ RandomSample[zerosAroundPeak, 100], PlotRange -> All, Frame -> True, GridLines -> {{0}, {-1, 1, 3}}, ImageSize -> 400] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]

Curves that start at the zero function value and attain a function value of ≈3 will soon move apart. The next graphic shows the mean value of spread of the distances after a selected distance. This spread is the reason that the zero distances that appear after 2 root3 are not peaks in the distribution of zero distances.

✕ Module[{zeroPosis, data}, zeroPosis = Select[Last /@ Select[Transpose[{differencesGolden, Range[Length[zerosGolden] - 1]}], 2 root3 - 0.001 < #[[1]] < 2 root3 + 0.001 &], # < 99000 &]; data = Table[{j, Mean[#] \[PlusMinus] StandardDeviation[#]} &[ differencesGolden[[zeroPosis + j]]], {j, 100}]; Graphics[{RGBColor[0.36, 0.51, 0.71], Line[{{#1, Subtract @@ #2}, {#1, Plus @@ #2}} & @@@ data]}, AspectRatio -> 1/2, Frame -> True]]

Instead of looking for the distance of successive roots of , one could look at the roots with an arbitrary right-hand side . Based on the above graphics, the most interesting distribution might occur for . Similar to the two-summand case, the distance between the zeros of the fastest component, namely , dominates the distribution. (Note the logarithmic vertical scale.)

✕ Histogram[ Differences[findZeros[fGolden, 10^5, -1]], 1000, {"Log", "Count"}, PlotRange -> All]

The above plot of seemed to show that the smallest values attained are around –1.5. This is indeed the absolute minimum possible; this follows from the fact that the summands lie on the Cayley surface.

✕ Minimize[{x + y + z, x^2 + y^2 + z^2 == 1 + 2 x y z \[And] -1 <= x <= 1 \[And] -1 <= y <= 1 \[And] -1 <= z <= 1}, {x, y, z}]

I use findZeros to find some near-minima positions. Note the relatively large distances between these minima. Because of numerical errors, one sometimes gets two nearby values, in which case duplicates are deleted.

✕ zerosGoldenMin = findZeros[fGolden, 100, -3/2] // DeleteDuplicates[#, Abs[#1 - #2] < 0.1 &] &;

Close to these absolute minima positions, the function takes on two universal shapes. This is clearly visible by plotting in the neighborhoods of all 100 minima.

✕ Show[Plot[fGolden[# + t], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ zerosGoldenMin, GridLines -> {{}, {-3/2, 3}}]

Interlude II: Using Rational Approximations of ϕ

The structure in the zero-distance distribution is already visible in relatively low-degree rational approximations of .

✕ fGoldenApprox[x_] = fGolden[x] /. GoldenRatio -> Convergents[GoldenRatio, 12][[-1]]

✕ zerosGoldenApprox = findZeros[fGoldenApprox, 100000];

✕ Histogram[Differences[zerosGoldenApprox], 1000, PlotRange -> All]

For small rational values of , in , one can factor the expression and calculate the roots of each factor symbolically to more quickly generate the list of zeros.

✕ getFirstZeros[f_[x_] + f_[\[Alpha]_ x_] + f_[\[Beta]_ x_] , n_] := Module[{sols, zs}, sols = Select[x /. Solve[f[x] + f[\[Alpha] x] + f[\[Beta] x] == 0, x] // N // ExpandAll // Chop, FreeQ[#, _Complex, \[Infinity]] &]; zs = getZeros[#, Ceiling[n/Length[sols]]] & /@ sols; Sort@Take[Flatten[zs], n]]

✕ getZeros[ConditionalExpression[a_. + C[1] b_, C[1] \[Element] Integers], n_] := a + b Range[n]

✕ Table[Histogram[ Differences[ getFirstZeros[ Cos[x] + Cos[\[Alpha] x] + Cos[\[Alpha]^2 x] /. \[Alpha] -> Convergents[GoldenRatio, d][[-1]], 100000]], 1000, PlotRange -> All, PlotLabel -> Row[{"convergent" , " \[Rule] ", d}]], {d, 2, 5}] // Partition[#, 2] &

For comparison, here are the corresponding plots for .

✕ Table[Histogram[ Differences[ getFirstZeros[ Sin[x] + Sin[\[Alpha] x] + Sin[\[Alpha]^2 x] /. \[Alpha] -> Convergents[GoldenRatio, d][[-1]], 100000]], 1000, PlotRange -> All, PlotLabel -> Row[{"convergent" , " \[Rule] ", d}]], {d, 2, 5}] // Partition[#, 2] &

The Extrema of f ϕ (x)

I now repeat some of the above visualizations for the extrema instead of the zeros.

✕ findExtremas[f_, n_] := findZeros[f', n]

✕ extremasGolden = Prepend[findExtremas[fGolden, 100000], 0.];

Here is a plot of together with the extrema, marked with the just-calculated values.

✕ Plot[fGolden[x], {x, 0, 10 Pi}, Epilog -> {Darker[Red], Point[{#, fGolden[#]} & /@ Take[ extremasGolden, 30]]}]

The extrema in the plane shows that extrema cluster around .

✕ Graphics[{PointSize[0.001], RGBColor[0.36, 0.51, 0.71], Opacity[0.1], Point[N[{#, fGolden[#]} & /@ extremasGolden]]}, AspectRatio -> 1/2, ImageSize -> 400, Frame -> True]

The distribution of the function values at the extrema is already visible in the last graphic. The peaks are at –3/2, –1 and 3. The following histogram shows how pronounced the density increases at these three values.

✕ Histogram[fGolden /@ extremasGolden, 1000, GridLines -> {{-3/2, -1, 3}, {}}, PlotRange -> All]

In a 3D histogram, one sees a strong correlation between the position of the extrema mod and the function value at the extrema.

✕ Histogram3D[{Mod[#, 2 Pi], fGolden[#]} & /@ extremasGolden, 100, AxesLabel -> {\[CapitalDelta]x, Subscript[f, \[Phi]], None}]

If one separates the minima and the maxima, one obtains the following distributions.

✕ Function[lg, Histogram3D[{Mod[#, 2 Pi], fGolden[#]} & /@ Select[extremasGolden, lg[fGolden''[#], 0] &], 100, AxesLabel -> {\[CapitalDelta]x, Subscript[f, \[Phi]], None}]] /@ {Less, Greater}

The function values of extrema ordinates are strongly correlated to the function values of successive extrema.

✕ Table[Histogram3D[{#[[1]], #[[-1]]} & /@ Partition[fGolden /@ extremasGolden, j, 1], 100, PlotLabel -> Row[{"shift:", j - 1}]], {j, 2, 4}]

A periodogram of the function values at the extrema shows a lot of structure. The various periodicities arising from the three cosine terms and their interference terms become visible. (The equivalent curve for the slope at the zeros is much noisier.) For all , , the periodogram of the function values at the extrema is a relatively smooth curve with only a few structures.

✕ Periodogram[fGolden /@ extremasGolden, PlotStyle -> RGBColor[0.88, 0.61, 0.14], Filling -> Axis]

Here again are two graphics that show how the three terms contribute to the function value at the maximum. The value of the term at the extrema is quite limited.

✕ summandValuesE = {Cos[#], Cos[GoldenRatio #], Cos[GoldenRatio^2 #]} & /@ extremasGolden;

✕ Graphics3D[{PointSize[0.002], RGBColor[0.36, 0.51, 0.71], Point[Union@Round[summandValuesE, 0.01]]}, PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, Axes -> True, AxesLabel -> {Subscript[f, \[Phi]], Cos[GoldenRatio x], Cos[GoldenRatio^2 x]}]

✕ Histogram[#, 200] & /@ Transpose[summandValuesE]

The two disconnected supports arise from the contributions at the minima and maxima. The contribution at the maxima is quite localized.

✕ Histogram[Cos[GoldenRatio^2 #], 200, PlotRange -> {{-1, 1}, All}] & /@ {Select[extremasGolden, fGolden''[#] > 0 &], Select[extremasGolden, fGolden''[#] < 0 &]}

I zoom into the second graphic of the last result.

✕ Histogram[ Cos[GoldenRatio^2 #] & /@ Select[extremasGolden, fGolden''[#] < 0 &], 200, PlotRange -> {{0.94, 1}, All}]

Here is the distribution of the function values at the minima and maxima.

✕ {Histogram[fGolden /@ Select[extremasGolden, fGolden''[#] > 0 &], 200], Histogram[fGolden /@ Select[extremasGolden, fGolden''[#] < 0 &], 200]}

Reduced to the interval , one sees a small, smooth peak around .

✕ Histogram[Mod[extremasGolden, 2 Pi], 200]

And this is the overall distribution of the distances between successive extrema. Compared with the zeros, it does not show much unexpected structure.

✕ Histogram[Differences[extremasGolden], 1000, PlotRange -> All]

Higher differences do show some interesting structure.

✕ Histogram[Differences[extremasGolden, 2], 1000, PlotRange -> All]

In an analogy to the zeros, I also show the pair correlation function. (In this context, this is also called a peak-to-peak plot.)

✕ Histogram3D[{#2 - #1, #3 - #2} & @@@ Partition[extremasGolden, 3, 1], 100, PlotRange -> All]

And here are the shapes of near the extrema for 1,000 randomly selected extrema.

✕ Show[Plot[fGolden[# + t], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ \ RandomSample[extremasGolden, 1000], PlotRange -> All, Frame -> True, GridLines -> {{0}, {3}}, ImageSize -> 400] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]

For use in the next section and in the next example, I form a list of the zeros and extrema in the order in which they are on the real axis.

✕ zerosAndExtremaGolden = Sort[Join[{#, "zero"} & /@ zerosGolden, {#, "extrema"} & /@ extremasGolden]];

Using this list, one can make a correlation plot of the slope at a zero with the height of the following maximum. One sees quite a strong correlation.

✕ Histogram3D[{fGolden'[#1[[1]]], fGolden[#2[[1]]]} & @@@ Cases[Partition[zerosAndExtremaGolden, 2, 1], {{_, "zero"}, {_, "extrema"}}], 120]

There is even a strong correlation between the slope at a zero and the second next extrema value.

✕ Histogram3D[{fGolden'[#[[1, 1]]], fGolden[Cases[#, {_, "extrema"}][[2, 1]]]} & /@ Take[Cases[ Partition[zerosAndExtremaGolden, 5, 1], {{_, "zero"}, __}], 60000], 120]

Interlude III: The Complex Zeros of f ϕ (x)

All zeros of are on the real axis, while does have complex zeros—exactly at the maximum in the negative maximum-minimum-maximum triples. A plot of shows zeros as peaks. These pairs of zeros with nonvanishing imaginary parts result in relatively large spacing between successive roots in . (The appearance of these complex-valued zeros does not depend on , being irrational, but also occurs for rational , .)

✕ Plot3D[Evaluate[-Log[Abs[fGolden[x + I y]]]], {x, 0, 12}, {y, -1, 1}, PlotPoints -> {60, 60}, BoxRatios -> {3, 1, 1}, MeshFunctions -> {#3 &}, ViewPoint -> {1.9, -2.4, 1.46}, MeshStyle -> RGBColor[0.36, 0.51, 0.71], AxesLabel -> {x, y, Subscript[f, \[Phi]][x + I y]}]

A contour plot showing the curves of vanishing real and imaginary parts has the zeros as crossings of curves of different color.

✕ ContourPlot[ Evaluate[# == 0 & /@ ReIm[fGolden[x + I y]]], {x, 0, 24}, {y, -1, 1}, PlotPoints -> {60, 60}, AspectRatio -> 1/3, Epilog -> {Purple, PointSize[0.012], Point[ ReIm[ z /. Solve[ fGolden[z] == 0 \[And] -2 < Im[z] < 2 \[And] 0 < Re[z] < 24]]]}]

Similar to the above case along the real axis, here is the path of convergence in Newton iterations.

✕ Module[{ppx = 111, ppy = 21, f = Evaluate[(# - fGolden[#]/fGolden'[#])] &, data, vals, color}, data = Table[ NestList[f, N[x0 + I y0], 20], {y0, -1, 1, 2/ppy}, {x0, 0, 30, 30/ppx}]; vals = First /@ Select[Reverse[ SortBy[Tally[Flatten[Round[Map[Last, data, {2}], 0.001]]], Last]], #[[2]] > 50 &]; (color[#] = Blend[{{0, RGBColor[0.88, 0.61, 0.14]}, {1/2, RGBColor[0.36, 0.51, 0.71]}, {1, RGBColor[0.561, 0.69, 0.19]}}, RandomReal[]]) & /@ vals; Graphics3D[{Thickness[0.001], Opacity[0.5], Map[If[MemberQ[vals, Round[Last[#], 0.001]], {color[Round[Last[#], 0.001]], BSplineCurve[MapIndexed[Append[#, #2[[1]]] &, ReIm[#]]]}, {}] &, data, {2}]}, BoxRatios -> {3, 1, 2}, ViewPoint -> {0.95, -3.17, 0.70}]]

Now that we have the positions of the zeros and extrema, we can also have a look at the complex zeros in more detail. In the above plots of over the complex plane, we saw that complex zeros occur when we have the situation maximum-minimum-maximum with all three function values at these negative extrema. Using the zeros and the extrema, we can easily find the positions of these extrema triples.

✕ tripleExtrema = SequenceCases[ zerosAndExtremaGolden, {{_, "extrema"}, {_, "extrema"}, {_, "extrema"}}][[All, 2, 1]];

About one-seventh of all these roots have a nonvanishing imaginary part.

✕ Length[tripleExtrema]/(Length[tripleExtrema] + Length[zerosGolden]) // N

The function value of the maximum of all of these consecutive triple extrema are indeed all negative and never smaller than –1.

✕ Histogram[fGolden /@ tripleExtrema, 100]

A log-log histogram of the absolute values of these middle maxima suggests that over a large range, a power law relation between their frequencies and their function values holds.

✕ Histogram[ Sort[-Select[fGolden /@ tripleExtrema, Negative]], {"Log", 100}, {"Log", "Count"}]

Interestingly, the distance between the middle maxima is narrowly concentrated at three different distances.

✕ tripleExtremaDifferences = Differences[tripleExtrema];

✕ Histogram[tripleExtremaDifferences, 100]

The next three plots zoom into the last three localized structures.

✕ Function[x, Histogram[Select[tripleExtremaDifferences, x - 1 < # < x + 1 &], 100]] /@ {5, 7, 12}

Here are the three typical shapes that belong to these three classes of distances. I show 100 randomly selected and shifted pieces of .

✕ extremaGroups = With[{L = Sort[Transpose[{tripleExtremaDifferences, Most[tripleExtrema]}]]}, Function[x, {x, Select[L, x - 1 < #[[1]] < x + 1 & ]}] /@ {5, 7, 12}];

✕ Function[{\[CapitalDelta], l}, Show[Plot[Evaluate[fGolden[# + t]], {t, -2, \[CapitalDelta] + 2}, PlotRange -> All, PlotLabel -> Row[{HoldForm[\[CapitalDelta]x == \[CapitalDelta]]}], PlotStyle -> Directive[Thickness[0.001], RGBColor[0.36, 0.51, 0.71], Opacity[0.6]]] & /@ RandomSample[l[[All, 2]], 100]]] @@@ extremaGroups

I use the position of all middle maxima as starting values for a numerical root-finding procedure to get the nearby complex roots.

✕ cRoots = If[Im[#] < 0, Conjugate[#]] & /@ z /. FindRoot[Evaluate[fGolden[z] == 0], {z, # + 1/2 I}] & /@ tripleExtrema;

Large imaginary parts of the complex zeros occur at real parts that are nearby multiples of π.

✕ Histogram3D[{Mod[Re[#], 2 Pi], Im[#]} & /@ cRoots, 100, AxesLabel -> {Im[x], Mod[Re[x], 2 Pi]}]

The function value at a middle maximum is strongly correlated with the magnitude of the imaginary part of the nearby complex root.

✕ Histogram3D[Transpose[{fGolden /@ tripleExtrema, Im[cRoots]}], 100, AxesLabel -> {Subscript[f, \[Phi]], Im[z]}]

And here are the distributions of imaginary parts and the differences of the imaginary parts of consecutive (along the real axis) zeros.

✕ {Histogram[Im[cRoots], 100], Histogram[Differences[Im[cRoots]], 100]}

If one splits the complex roots into the three groups from above, one obtains the following distributions.

✕ Column[GraphicsRow[{Histogram[Im[#], 100], Histogram[Differences[Im[#]], 100]} &[ If[Im[#] < 0, Conjugate[#]] & /@ z /. FindRoot[Evaluate[fGolden[z] == 0], {z, # + 1/2 I}] & /@ #2[[ All, 2]]], ImageSize -> 360, PlotLabel -> Row[{HoldForm[\[CapitalDelta]x == #1]}]] & @@@ extremaGroups]

As a function of , the complex zeros of periodically join the real axis. The following graphic shows the surface in the -space in yellow/brown and the complex zeros as blue tubes. We first plot the surface and then plot the curves as mesh lines on this surface.

✕ ContourPlot3D[ Re[Cos[x + I y] + Cos[\[Alpha] (x + I y)] + Cos[\[Alpha]^2 (x + I y)]], {x, 0, 12}, {y, 0, 1}, {\[Alpha], 1, 2}, Contours -> {0}, BoxRatios -> {3, 1, 1}, ContourStyle -> {Directive[RGBColor[0.88, 0.61, 0.14], Opacity[0.3], Specularity[RGBColor[0.36, 0.51, 0.71], 10]]}, BoundaryStyle -> None, ImageSize -> 600, PlotPoints -> {160, 80, 80}, MaxRecursion -> 0, MeshStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.006]], Mesh -> {{0}}, MeshFunctions -> {Function[{x, y, \[Alpha]}, Im[Cos[x + I y] + Cos[\[Alpha] (x + I y)] + Cos[\[Alpha]^2 (x + I y)]]]}, AxesLabel -> {x, y, \[Alpha]}]

The next natural step (and the last one for this section) would be to look at the defined implicitly through . As for any value of , we will have (possibly infinitely) many ; we want to construct the Riemann surface for . A convenient way to calculate Riemann surfaces is through solving the differential equation for the defining function. Through differentiation, we immediately get the differential equation.

✕ f\[Alpha][z_, \[Alpha]_] := Cos[z] + Cos[\[Alpha] z] + Cos[\[Alpha]^2 z]

✕ Solve[D[f\[Alpha][z[\[Alpha]], \[Alpha]], \[Alpha]] == 0, Derivative[1][z][\[Alpha]]] // Simplify

✕ rhs[z_, \[Alpha]_] := -(( z (Sin[\[Alpha] z] + 2 \[Alpha] Sin[\[Alpha]^2 z]) )/( Sin[z] + \[Alpha] (Sin[\[Alpha] z] + \[Alpha] Sin[\[Alpha]^2 z])))

As starting values, I use along a segment of the real axis.

✕ ICs\[Alpha]x = Cases[Flatten[ Table[{#, x} & /@ (\[Alpha] /. Quiet[Solve[ f\[Alpha][N[1, 40] x, \[Alpha]] == 0 \[And] 0 < \[Alpha] < 2, \[Alpha]]]), {x, 0, 6, 6/100}], 1], {_Real, _}];

These are the points that we will use for the numerical differential equation solving.

✕ cp = ContourPlot[ f\[Alpha][x, \[Alpha]] == 0, {x, 0, 6}, {\[Alpha], 0, 1.9}, Epilog -> {Purple, Point[Reverse /@ ICs\[Alpha]x]}, FrameLabel -> {x, \[Alpha]}]

Starting at a point , we move on a semicircle around the origin of the complex plane. To make sure we stay on the defining Riemann surface, we use the projection method for the numerical solution. And we change variables from , where .

✕ Monitor[rsf\[Alpha]Data = Table[With[{r = ICs\[Alpha]x[[k, 1]]}, nds\[Alpha] = NDSolveValue[{Derivative[1][z][\[CurlyPhi]] == I r Exp[I \[CurlyPhi]] rhs[z[\[CurlyPhi]], r Exp[I \[CurlyPhi]]], z[0] == ICs\[Alpha]x[[k, 2]]}, z, {\[CurlyPhi], 0 , Pi }, WorkingPrecision -> 30, PrecisionGoal -> 6, AccuracyGoal -> 6, Method -> {"Projection", Method -> "StiffnessSwitching", "Invariants" -> f\[Alpha][z[\[CurlyPhi]], r Exp[I \[CurlyPhi]]]}, MaxStepSize -> 0.01, MaxSteps -> 10^5]; {{r, nds\[Alpha][0], nds\[Alpha][Pi]}, ParametricPlot3D[ Evaluate[{r Cos[\[CurlyPhi]], r Sin[\[CurlyPhi]], #[nds\[Alpha][\[CurlyPhi]]]}], Evaluate[Flatten[{\[CurlyPhi], nds\[Alpha][[1]]}]], BoxRatios -> {1, 1, 1}, Axes -> True, PlotRange -> All, PlotStyle -> Directive[Thickness[0.002], Darker[Blue]], ColorFunctionScaling -> False, ColorFunction -> (ColorData["DarkRainbow"][ Abs[1 - #4/ Pi]] &)] & /@ {Re, Im}}], {k, Length[ICs\[Alpha]x]}] // Quiet;, Text@ Row[{"path ", k, " of ", Length[ICs\[Alpha]x]}]]

Here is a plot of the real part of . One clearly sees how the sets of initially nearby zeros split at branch points in the complex plane.

✕ rsfRe = Show[rsf\[Alpha]Data[[All, 2, 1]], PlotRange -> {All, {0, All}, 10 {-1, 1}}, AxesLabel -> {Re[\[Alpha]], Im[\[Alpha]], Re[z]}, ViewPoint -> {-0.512, -3.293, 0.581}]

One can calculate the branch points numerically as the simultaneous solutions of and .

✕ branchPointEqs[z_, \[Alpha]_] = {f\[Alpha][z, \[Alpha]], D[f\[Alpha][z, \[Alpha]], z]}

✕ branchPoints = Union[Round[Select[#, Total[Abs[branchPointEqs @@ #]] < 10^-10 &], 10.^-6]] &@ ( Table[{z, \[Alpha]} /. FindRoot[Evaluate[branchPointEqs[z, \[Alpha]] == {0, 0}] , {z, RandomReal[{-8, 8}] + I RandomReal[{-8, 8}] }, {\[Alpha], RandomReal[{-5, 5}] + I RandomReal[{-5, 5}] }, PrecisionGoal -> 10] // Quiet, {20000}]);

Here are the branch points near the origin of the complex plane.

✕ Graphics[{RGBColor[0.36, 0.51, 0.71], Point[ReIm /@ Last /@ branchPoints]}, Frame -> True, PlotRange -> 3.5, FrameLabel -> {Re[\[Alpha]], Im[\[Alpha]]}, PlotRangeClipping -> True]

Representing the positions of the branch points in the above plot as vertical cylinders shows that the splitting indeed occurs at the branch points (we do not include the branch points with to have a better view of this complicated surface).

✕ Show[{rsfRe, Graphics3D[{CapForm[None], GrayLevel[0.3], Specularity[Purple, 20], Cylinder[{Append[ReIm[#], -10], Append[ReIm[#], 10]}, 0.02] & /@ Select[Last /@ branchPoints, Im[#] > 10^-2 &]}]}, PlotRange -> {{-2, 2}, {0, 2}, 8 {-1, 1}}, ViewPoint -> {-1.46, 2.87, 1.05}]

Finding More Envelopes

The possible “shapes” of fGolden near points where the function has an extremum and a value ≈±1 is quite limited. The three possible curves arise from the different ways to form the sum –1 from the values ±1 of the three summands.

✕ maxMinus1Value = Select[extremasGolden, Abs[fGolden[#] - (-1)] < 10^-4 &]; Length[maxMinus1Value]

✕ Take[maxMinus1Value, 12]

✕ Show[Plot[fGolden[# + t], {t, 0 - 2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ maxMinus1Value, PlotRange -> All, Frame -> True, GridLines -> {{0}, {3}}, ImageSize -> 400]

We compare these three curves with the possible curves that could be obtained from considering all possible phases between the three summands—meaning we consider all + + such that and .

✕ gGolden[x_, {\[CurlyPhi]1_, \[CurlyPhi]2_, \[CurlyPhi]3_}] := Cos[x + \[CurlyPhi]1] + Cos[GoldenRatio x + \[CurlyPhi]2] + Cos[GoldenRatio^2 x + \[CurlyPhi]3];

✕ phaseList = Table[ Block[{\[CurlyPhi]1 = RandomReal[{-Pi, Pi}]}, Check[Flatten[{\[CurlyPhi]1, {\[CurlyPhi]2, \[CurlyPhi]3} /. FindRoot[{Cos[\[CurlyPhi]1] + Cos[\[CurlyPhi]2] + Cos[\[CurlyPhi]3] == -1, -Sin[\[CurlyPhi]\ 1] - GoldenRatio Sin[\[CurlyPhi]2] - GoldenRatio^2 Sin[\[CurlyPhi]3] == 0}, {{\[CurlyPhi]2, RandomReal[{-Pi, Pi}]}, {\[CurlyPhi]3, RandomReal[{-Pi, Pi}]}}]}], {}]], {600}]; // Quiet

The next graphic shows the two curve families together.

✕ Show[ {Plot[Evaluate[ gGolden[x, #] &@ DeleteCases[Take[phaseList, All], {}]], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Opacity[0.3], Thickness[0.001]]], Plot[fGolden[t - #], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.88, 0.61, 0.14], Thickness[0.003], Opacity[0.3]]] & /@ maxMinus1Value}]

These curves now allow us to answer the original question about the location of the singularities in the distribution of the distances of successive zeros. Similar to the largest peak identified above arising from the bunching of curves with function values , the other three maxima arise from curves with local minima or maxima . We calculate some of these zeros numerically.

✕ minEnvelopeZeros = Union[Round[ Sort[\[Delta] /. Quiet[FindRoot[ fGolden[# + \[Delta]] == 0, {\[Delta], 1/2}]] & /@ Take[maxMinus1Value, 100]], 0.025]]

Indeed, the gridlines match the singularities precisely.

✕ Histogram[Differences[zerosGolden], 1000, PlotRange -> All, GridLines -> {Flatten[{2 zerosGolden[[1]], 2 minEnvelopeZeros}], {}}]

The unifying feature of all four singularities is their location at the zeros of envelope curves. Here are the curves of fGolden around 100 zeros for each of the four singularities.

✕ With[{L = {#[[1, 1]], Last /@ #} & /@ Reverse[SortBy[Split[Sort[Transpose[ {Round[differencesGolden, 0.002], Most[zerosGolden]}]], #1[[ 1]] === #2[[1]] &], Length]]}, Partition[Function[p, Show[ Plot[fGolden[# + t], {t, -2, 0 + 5}, PlotRange -> All, PlotLabel -> Row[{"x", "\[TildeTilde]", p}], GridLines -> {None, {-1, 3}}, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ Take[Select[L, Abs[#[[1]] - p] < 0.1 &, 1][[1, 2]], UpTo[100]]]] /@ {1, 1.5, 1.8, 3}, 2]] // Grid

These graphics show strikingly the common feature of these four groups of zero distances: either maxima cluster around or minima cluster around .

The curves in the plane that fulfill the two conditions are shown in the next contour plot.

✕ ContourPlot[ Evaluate[Function[sign, Derivative[1, {0, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0 /. (Solve[ gGolden[0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == sign 1, \[CurlyPhi]1] /. ConditionalExpression[x_, _] :> x /. _C :> 0)] /@ {-1, 1}], {\[CurlyPhi]2, 0, 2 Pi}, {\[CurlyPhi]3, 0, 2 Pi}, PlotPoints -> 60, ContourStyle -> RGBColor[0.36, 0.51, 0.71]]

Now we can also calculate a numerical approximation to the exact value of the position of the first singularity. We consider the envelope of all curves with the properties , . The envelope property is represented through the vanishing partial derivative with respect to .

✕ (Table[Round[ Mod[{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3} /. FindRoot[ Evaluate[{gGolden[ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == -1, Derivative[1, {0, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[0, {1, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0}], {\[CurlyPhi]1, RandomReal[{-Pi, Pi}]}, {\[CurlyPhi]2, RandomReal[{-Pi, Pi}]}, {\[CurlyPhi]3, RandomReal[{-Pi, Pi}]}, PrecisionGoal -> 12] , 2 Pi] /. x_?(Abs[#] < 10^-6 || Abs[2 Pi - #] < 10.^-6 &) :> 0, 10.^-6], {100}] // Quiet // Tally) /. x_Real?(Abs[# - Pi] < 0.01 &) :> Pi

The envelope curve of this family is a small modification of our original function—the sign of the first term is inverted.

✕ cosEqs = {gGolden[x, {Pi, Pi, 0}], gGolden[x, {0, Pi, Pi}], gGolden[x, {Pi, 0, Pi}]}

The positions of the first roots of the last equations are , and . Multiplying these numbers by two to account for their symmetric distance from the origin yields the zero distances seen above at , and .

✕ FindRoot[# == 0, {x, 1/2}] & /@ cosEqs

And, using the higher-precision root, we zoom into the neighborhood of the first peak to confirm our conjecture.

✕ With[{z0 = 2 x /. FindRoot[ Evaluate[Cos[x] - Cos[GoldenRatio x] - Cos[GoldenRatio^2 x] == 0], {x, 1/2}, PrecisionGoal -> 12]}, Histogram[ Select[Differences[ zerosGoldenRefined], -0.0003 < z0 - # < 0.0004 &] - z0, 100, "PDF", GridLines -> {{0}, {}}]]

The envelope arises from the fact that a small variation in can be compensated by appropriate changes in and .

✕ ((-\[Alpha] + \[Beta]*Sqrt[\[Alpha]^2 + \[Beta]^2 - 1])*Sin[x*\[Alpha]])/(\[Alpha]^2 + \[Beta]^2) + ((-\[Beta] - \[Alpha]*Sqrt[\[Alpha]^2 + \[Beta]^2 - 1])*Sin[x*\[Beta]])/ (\[Alpha]^2 + \[Beta]^2) // FullSimplify

✕ ((-\[Alpha] + \[Beta] Sqrt[\[Alpha]^2 + \[Beta]^2 - 1]) Sin[x \[Alpha]] - (\[Beta] + \[Alpha] Sqrt[\[Alpha]^2 + \[Beta]^2 - 1]) Sin[x \[Beta]])/(\[Alpha]^2 + \[Beta]^2)

✕ Manipulate[ Plot[Evaluate[{-Cos[x] + Cos[\[Alpha] x] + Cos[\[Beta] x], -Cos[x] + Cos[x \[Alpha]] + Cos[x \[Beta]] + (((-\[Alpha] + \[Beta] Sqrt[\[Alpha]^2 + \[Beta]^2 \ - 1]) Sin[x \[Alpha]] - (\[Beta] + \[Alpha] Sqrt[\[Alpha]^2 + \[Beta]^2 - 1]) Sin[ x \[Beta]])/(\[Alpha]^2 + \[Beta]^2)) \ \[CurlyPhi]1MinusPi} /. { \[Alpha] -> GoldenRatio, \[Beta] -> GoldenRatio^2} /. \[CurlyPhi]1 -> Pi + 0.1], {x, 0, 1}], {{\[CurlyPhi]1MinusPi, 0.3, HoldForm[\[CurlyPhi]1 - Pi]}, -0.5, 0.5, Appearance -> "Labeled"}, TrackedSymbols :> True]

These zeros of the equation do indeed match the position of the above singularities in the distribution of the zero distances.

We end this section by noting that envelope condition, meaning the vanishing of derivatives with respect to the phases , is a very convenient method for finding special families of curves. For instance, there are families of solutions such that many curves meet at a zero.

✕ findSolutions[eqs_] := Module[{fm}, While[fm = FindRoot[Evaluate[eqs], {\[CurlyPhi]1, RandomReal[{0, 2 Pi}]}, {\[CurlyPhi]2, RandomReal[{0, 2 Pi}]}, {\[CurlyPhi]3, RandomReal[{0, 2 Pi}]}, {d, RandomReal[{1/2, 5}]}] // Quiet; Not[ Total[Abs[(Subtract @@@ eqs) /. fm]] < 10^-6 \[And] 10^-3 < Abs[d /. fm] < 5] ] ; fm /. (\[CurlyPhi] : (\[CurlyPhi]1 | \[CurlyPhi]2 | \[CurlyPhi]3) \ -> \[Xi]_) :> (\[CurlyPhi] -> Mod[\[Xi], 2 Pi])]

Modulo reflection symmetry, there is a unique solution to this problem.

✕ SeedRandom[1]; fs1 = findSolutions[{gGolden[ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[0, {1, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, gGolden[d, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[0, {1, 0, 0}][gGolden][ d, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0}]

Plotting a family of curves with different values of shows that at the two zeros, the family of curves bunches up.

✕ Plot[Evaluate[ Table[gGolden[ x, {\[CurlyPhi]1 + \[Delta]\[CurlyPhi]1, \[CurlyPhi]2, \ \[CurlyPhi]3}], {\[Delta]\[CurlyPhi]1, -0.2, 0.2, 0.4/6}] /. fs1], {x, -3, 6}, GridLines -> ({{d}, {}} /. fs1)]

And here are some solutions that show curves that bunch at a zero and at an extremum.

✕ Function[sr, SeedRandom[sr]; fs2 = findSolutions[{gGolden[ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[0, {1, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[1, {0, 0, 0}][gGolden][ d, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[0, {1, 0, 0}][gGolden][ d, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0}]; Plot[Evaluate[ Table[gGolden[ x, {\[CurlyPhi]1 + \[Delta]\[CurlyPhi]1, \[CurlyPhi]2, \ \[CurlyPhi]3}], {\[Delta]\[CurlyPhi]1, -0.2, 0.2, 0.4/6}] /. fs2], {x, -3, 6}, GridLines -> ({{d}, {}} /. fs2)]] /@ {1, 6, 38}

We also look at a second family of envelope solutions: what are the most general types of curves that fulfill the envelope condition at an extremum? This means we have to simultaneously fulfill the following two equations.

✕ {Derivative[1, {0, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0, Derivative[0, {1, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0}

This in turn means we must have .

✕ cpExtEnv = ContourPlot[-GoldenRatio Sin[\[CurlyPhi]2] - GoldenRatio^2 Sin[\[CurlyPhi]3] == 0, {\[CurlyPhi]2, 0, 2 Pi}, {\[CurlyPhi]3, 0, 2 Pi}, PlotPoints -> 40]

Here is a contour plot of the curves in the plane where the extremum envelope condition is fulfilled.

✕ linesExtEnv = Cases[Normal[cpExtEnv], _Line, \[Infinity]];

✕ Show[ Plot[ Cos[x + 0] + Cos[GoldenRatio x + #1] + Cos[GoldenRatio^2 x + #2], {x, -5, 5}, PlotStyle -> Directive[Thickness[0.001], RGBColor[0.36, 0.51, 0.71], Opacity[0.5]], PlotRange -> All, GridLines -> {{}, {-3, -1, 1, 3}}] & @@@ \ RandomSample[Level[linesExtEnv, {-2}], UpTo[300]]]

We add a mouseover to the contour plot that shows a family of curves near to the envelope conditions.

✕ make\[CurlyPhi]2\[CurlyPhi]3Plot[{\[CurlyPhi]2_, \[CurlyPhi]3_}] := Column[{Plot[ Evaluate[ Table[gGolden[ x, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}], {\[CurlyPhi]1, \ -0.15, 0.15, 0.3/6}]], {x, -4, 4}, PlotLabel -> Subscript[\[CurlyPhi], 1] \[TildeTilde] 0], Plot[ Evaluate[ Table[gGolden[ x, {Pi + \[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}], {\ \[CurlyPhi]1, -0.15, 0.15, 0.3/6}]], {x, -4, 4}, PlotLabel -> Subscript[\[CurlyPhi], 1] \[TildeTilde] Pi]}, Dividers -> Center]

✕ Normal[cpExtEnv] /. Line[l_] :> (Tooltip[Line[ #], Dynamic[make\[CurlyPhi]2\[CurlyPhi]3Plot[Mean[#]]]] & /@ Partition[l, 2, 1] )

And the next graphic shows the distance between the nearest (to the origin) positive or negative roots.

✕ addHeight[{\[CurlyPhi]2_, \[CurlyPhi]3_}, pm_] := Module[{roots, \[Rho] }, roots = Sort[x /. Solve[gGolden[x, {0, \[CurlyPhi]2, \[CurlyPhi]3}] == 0 \[And] -6 < x < 6, x]] // Quiet; \[Rho] = If[pm === 1, Select[roots, # > 0 &, 1], -Select[Sort[-roots], # > 0 &, 1] ][[ 1]]; {{\[CurlyPhi]2, \[CurlyPhi]3, 0}, {\[CurlyPhi]2, \[CurlyPhi]3, \[Rho]}}]

✕ Graphics3D[{Opacity[0.2], Gray, Polygon[{{0, 0, 0}, {2 Pi, 0, 0}, {2 Pi, 2 Pi, 0}, {0, 2 Pi, 0}}], EdgeForm[], Opacity[1], Transpose[{{RGBColor[0.88, 0.61, 0.14], RGBColor[0.36, 0.51, 0.71]}, Table[ Polygon[Join[#1, Reverse[#2]] & @@@ Partition[addHeight[#, s] & /@ #[[1]], 2, 1]] & /@ linesExtEnv, {s, {1, -1}}]}]}, PlotRange -> All, Axes -> True, Lighting -> "Neutral", AxesLabel -> {Subscript[\[CurlyPhi], 2], Subscript[\[CurlyPhi], 3], \!\(\*SubscriptBox[\(\[Rho]\), \("\"\)]\)} ]

The Peak Positions of the Zero Distances

Now I can implement the following function, maximaPositions , to find the singularities of the distributions of the distances of successive zeros for functions of the form + for arbitrary real .

✕ maximaPositions[ a1_. Cos[x_] + a2_. Cos[\[Alpha]_ x_] + a3_. Cos[\[Beta]_ x_], x_] := 2*(Function[{s1, s2, s3}, Min[N[ x /. Solve[ s1 a1 Cos[x] + s2 a2 Cos[\[Alpha] x] + s3 a3 Cos[\[Beta] x] == 0 \[And] 0 < x < 10, x]]]] @@@ {{1, 1, 1}, {-1, 1, 1}, {1, -1, 1}, {1, 1, -1}})

Here are the positions of the peaks for .

✕ peaksGolden = maximaPositions[fGolden[x], x]

Zooming into the histogram shows again that the predicted peak positions match the observed ones well.

✕ Partition[Histogram[Differences[zerosGolden], 2000, PlotRange -> {{# - 0.02, # + 0.02}, All}, GridLines -> {peaksGolden, {}}, PlotRangeClipping -> True] & /@ peaksGolden, 2] // Grid

Here is a quick numerical/graphical check for a “random” function of the prescribed form with not identically one.

✕ testFunction[x_] := Cos[x] + 2 Cos[Sqrt[2] x] + Sqrt[3] Cos[Pi x]

✕ zerosTestFunction = findZeros[testFunction, 10^5];

✕ singPosFunction = maximaPositions[testFunction[x], x]

✕ Histogram[Differences[zerosTestFunction], 100, GridLines -> {singPosFunction, {}}]

We should check for an extended range of parameters if the conjectured formula for the position of the singularities really holds. So for a few hundred values of in , we calculate the histograms of the zero distances. This calculation will take a few hours. (In the interactive version, the cell is set to unevaluatable.)

✕ With[{p\[Alpha] = 200}, Monitor[ \[Alpha]Data = SparseArray[ Flatten[Table[({j, Round[p\[Alpha] #1] + 1} -> #2) & @@@ #[[ j]][[2]], {j, Length[#]}], 1]] & @ Table[ f\[Alpha][x_] := Cos[x] + Cos[\[Alpha] x] + Cos[\[Alpha]^2 x]; zList = findZeros[f\[Alpha], 20000, 0]; {\[Alpha], Tally[Round[Differences[zList], N[1/p\[Alpha]]]]};, {\[Alpha], 1/2, 2, 1/(120 Pi)}], Row[{"\[Alpha]=", N[\[Alpha]]}]]]

✕ rPlot\[Alpha] = ReliefPlot[Log[1 + Take[\[Alpha]Data, {2, -1}, {1, 1200}]], DataRange -> {{0, 1200 0.005}, {1/2, 2}}, FrameTicks -> True, FrameLabel -> {"\[CapitalDelta] zeros", "\[Alpha]"}, AspectRatio -> 1/3]

I overlay with the predicted positions of the singularities; they match perfectly.

✕ Monitor[\[Alpha]Sings = Table[{#, \[Alpha]} & /@ maximaPositions[N[Cos[x] + Cos[\[Alpha] x] + Cos[\[Alpha]^2 x]], x], {\[Alpha], 0.5, 2, 1/5/201}];, \[Alpha]]

Now, what do the positions of the singularities in the case of two parameters , in look like? Plotting the locations of a few thousand singularity positions in space clearly shows four intersecting surfaces.

✕ singGraphics3D = Graphics3D[{RGBColor[0.36, 0.51, 0.71], Sphere[Flatten[ Table[{\[Alpha], \[Beta], #} & /@ Quiet[maximaPositions[ N[Cos[x] + Cos[\[Alpha] x] + Cos[\[Beta] x]], x]], {\[Beta], 1/2 + Pi/1000, 5/2 + Pi/1000, 2/31}, {\[Alpha], 1/2 + Pi/1001, 5/2 + Pi/1001, 2/31}], 2], 0.015]}, PlotRange -> {{0.5, 2.5}, {0.5, 2.5}, {0, 6}}, BoxRatios -> {1, 1, 1}, Axes -> True, ViewPoint -> {2.83, 1.76, -0.61}, Method -> {"SpherePoints" -> 6}]

Instead of calculating many concrete tuples of positions, we can consider the positions of the singularities as functions of , in . Differentiating the last expression allows us to get a partial differential equation for each of the four equations that, given boundary conditions, we can solve numerically and then plot. Or instead of solving a partial differential equation, we can fix one of the two parameters and and obtain a family of parametrized ordinary differential equations. We carry this procedure out for the two sheets that don’t cross each other.

✕ smallestRoot[{pm1_, pm2_, pm3_}, \[Alpha]_, \[Beta]_] := Min@Select[DeleteCases[Table[ Quiet@ Check[x /. FindRoot[ pm1 Cos[x] + pm2 Cos[\[Alpha] x] + pm3 Cos[\[Beta] x] == 0, {x, RandomReal[{0, Pi}]}, PrecisionGoal -> 12], {}], {40}], {}], # > 0 &]

✕ makeGrid[{pm1_, pm2_, pm3_}, M_] := Show[Cases[Flatten[ {Table[ nds\[Alpha] = NDSolveValue[{D[ Cos[x[\[Alpha]]] + Cos[\[Alpha] x[\[Alpha]]] + Cos[\[Beta] x[\[Alpha]]], \[Alpha]] == 0, x[E] == smallestRoot[{pm1, pm2, pm3}, E, \[Beta]]}, x, {\[Alpha], 3/4, 3}, PrecisionGoal -> 12]; ParametricPlot3D[{\[Alpha], \[Beta], 2 nds\[Alpha][\[Alpha]]}, Evaluate[Flatten[{\[Alpha], nds\[Alpha][[1]]}]], PlotRange -> All, PlotStyle -> Gray], {\[Beta], 3/4 + Pi/1000, 3 + Pi/1000, 2/M}], Table[ nds\[Beta] = NDSolveValue[{D[ Cos[x[\[Beta]]] + Cos[\[Alpha] x[\[Beta]]] + Cos[\[Beta] x[\[Beta]]], \[Beta]] == 0, x[E] == smallestRoot[{pm1, pm2, pm3}, \[Alpha], E]}, x, {\[Beta], 3/4, 3}, PrecisionGoal -> 12]; ParametricPlot3D[{\[Alpha], \[Beta], 2 nds\[Beta][\[Beta]]}, Evaluate[Flatten[{\[Beta], nds\[Beta][[1]]}]], PlotRange -> All, PlotStyle -> Directive[RGBColor[0.88, 0.61, 0.14], Thickness[0.002]]], {\[Alpha], 3/4 + Pi/1000, 3 + Pi/1000, 2/M}] }], _Graphics3D], PlotRange -> All, AxesLabel -> {"\[Alpha]", "\[Beta]", None}] // Quiet

✕ Show[{singGraphics3D, makeGrid[{1, 1, 1}, 20], makeGrid[{-1, 1, 1}, 20]}]

The other two sheets could also be calculated, but things would become a bit more complicated because the initial conditions calculated with smallestRoot are no longer continuous functions of and . The following graphic visualizes this situation when new zeros suddenly appear that didn’t exist for the blue curve due to a small change in .

✕ Module[{\[Alpha] = 2.02, \[Beta] = 2.704}, Plot[{Cos[x] - Cos[\[Alpha] x] + Cos[\[Beta] x], Cos[x] - Cos[(\[Alpha] - 0.1) x] + Cos[\[Beta] x]}, {x, 0, 3}]]

Interestingly, one can generalize even further the above formula for the peaks in the distance between successive zeros and allow arbitrary phases in the three cos functions. Numerical experiments indicate that for many such cosine sums, one can just ignore the phases and find the smallest zeros of exactly the same equations as above.

✕ maximaPositionsGeneralized[ a1_. Cos[x_ + \[CurlyPhi]1_.] + a2_. Cos[\[Alpha]_ x_ + \[CurlyPhi]2_.] + a3_. Cos[\[Beta]_ x_ + \[CurlyPhi]3_.], x_] := maximaPositions[ a1 Cos[x] + a2 Cos[\[Alpha] x] + a3 Cos[\[Beta] x], x]

Here is another random example.

✕ fRandom[x_] := 1/2 Cos[x + 1] + Sqrt[2] Cos[GoldenRatio x + Pi/5] + Sqrt[3] Cos[Pi x + 2] zerosRandom = findZeros[fRandom, 10^5];

✕ singPosFunction = maximaPositionsGeneralized[fRandom[x], x]

✕ Histogram[Differences[zerosRandom], 100, GridLines -> {singPosFunction, {}}]

I modify the phases in the above function fRandom , and obtain the same position for the singularities.

✕ fRandom2[x_] := 1/2 Cos[x + Cos[3]] + Sqrt[2] Cos[GoldenRatio x + Log[2]] + Sqrt[3] Cos[Pi x + E] zerosRandom2 = findZeros[fRandom2, 10^5];

✕ Histogram[Differences[zerosRandom2], 100, GridLines -> {singPosFunction, {}}]

In degenerate cases, such as all phases being , meaning , the positions of the peaks in the distribution of the zero distances will sometimes be different from the just-given conjecture. This means that e.g. the position of the spacings of the extrema of the sin equivalent of are not described by the above equations.

Even More Peaks—Well, Sometimes

In the last section, we established the position of the peaks of the original MathOverflow post as twice the smallest zeros of . And for many random instances of and , these four numbers indeed characterize the peaks visible in the distribution of successive zero distances of . We saw this clearly in the plots in the last section that spanned various parameter ranges for and .

Earlier, I remarked that the original MathOverflow plot that used the function has some features special to this function, like the gaps in the distribution of the zero distances. At the same time, the function is generic in the sense that the zero-spacing distribution has exactly four peaks, as we expect for generic and : the four red curves in the above graphic for rPlotα . We consider a slight modification of fGolden where instead of , we use .

✕ \[Phi]Brass = (1 + Sqrt[3])/2; fBrass[x_] := Cos[x] + Cos[\[Phi]Brass x] + Cos[\[Phi]Brass^2 x]

We again calculate the first 100k zeros and their distances.

✕ zerosBrass = findZeros[fBrass, 10^5];

✕ differencesBrass = Differences[zerosBrass];

Interestingly, we now get a distribution with six maxima. Four of the maxima are again well described by the maxima position conjectured above.

✕ peaksBrass = maximaPositions[fBrass[x], x];

✕ Histogram[differencesBrass, 1000, GridLines -> {peaksBrass, {}}]

The square root term makes all these examples special. The function is ultimately only dependent on two, rather than three, algebraically independent terms.

✕ Cos[x] + Cos[(1 + Sqrt[n])/2 x] + Cos[(1 + Sqrt[n])/2 ^2 x] // ExpandAll // TrigExpand

✕ Collect[Subtract @@ Eliminate[{fSqrtN == %, Cos[x/4]^2 + Sin[x/4]^2 == 1, Cos[Sqrt[n] x/4]^2 + Sin[Sqrt[n] x/4]^2 == 1}, {Sin[x/4], Sin[Sqrt[n] x/4]}] , fP, Factor]

This is reflected in the fact that the triples , similar to the above equivalent for the golden ratio, are located on relatively simple 2D surfaces.

✕ Show[{ContourPlot3D[(-1 + x + 2 y^2)^2 + 4 (-1 + x - 2 x y^2) z^2 + 4 z^4 == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, PlotPoints -> 40, Mesh -> None, ContourStyle -> Opacity[0.5]], Graphics3D[{PointSize[0.005], Point[ Table[{Cos[x] , Cos[\[Phi]Brass x] , Cos[\[Phi]Brass^2 x]}, {x, 1., 10000., 1.}]]}]}]

A tempting and easy generalization of the positions of the maxima might be zero distances that follow the previously calculated four distances. This generalization is easy to implement.

✕ rootDistances[ a1_. Cos[x_] + a2_. Cos[\[Alpha]_ x_] + a3_. Cos[\[Beta]_ x_], x_, xMax_] := Function[{s1, s2, s3}, DeleteDuplicates[Differences[ Sort[ N[x /. Solve[s1 a1 Cos[x] + s2 a2 Cos[\[Alpha] x] + s3 a3 Cos[\[Beta] x] == 0 \[And] -xMax < x < xMax, x]]]]]] @@@ {{1, 1, 1}, {-1, 1, 1}, {1, -1, 1}, {1, 1, -1}}

Interestingly, that gives nearly the correct positions, but not quite.

✕ rdBrass = rootDistances[fBrass[x], x, 12]

✕ Function[c, Histogram[Select[differencesBrass, c - 0.005 < # < c + 0.005 &], 100, GridLines -> {Flatten[rdBrass], {}}, PlotRangeClipping -> True]] /@ {1.0977, 2.366}

Where are these additional maxima? We will not calculate their exact positions here. But a quick look at the neighborhood of zeros that have the peak distances shows clearly that there are additional families of envelope curves involved. Interestingly again, families with occur, and this time reflected and translated versions of the curve arise.

✕ With[{L = {#[[1, 1]], Last /@ #} & /@ Reverse[SortBy[Split[Sort[Transpose[ {Round[differencesBrass, 0.001], Most[zerosBrass]}]], #1[[ 1]] === #2[[1]] &], Length]]}, Function[p, Show[Plot[fBrass[# + t], {t, -10, 0 + 10}, PlotRange -> All, PlotLabel -> Row[{"x", "\[TildeTilde]", p}], GridLines -> {None, {-1, 1, 3}}, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.5]]] & /@ Take[Select[L, Abs[#[[1]] - p] < 0.1 &, 1][[ 1, 2]], UpTo[50]]]] /@ {1.098, 2.367}] /. l_Line :> Mouseover[l, {Red, Thick, l}]

We will not try to interpret the family of curves that form these singularities here and now; in the next section we will develop some general methods to identify the positions of the singularities on the zero distances.

For larger square roots, even stranger distributions with even more maxima can occur in the distribution of the zero spacings. Here is the case:

✕ \[Phi]Platinum = (1 + Sqrt[13])/2; fPlatinum[x_] := Cos[x] + Cos[\[Phi]Platinum x] + Cos[\[Phi]Platinum^2 x] zerosPlatinum = findZeros[fPlatinum, 10^5];

The zeros mod have multiple gaps.

✕ Histogram[Mod[zerosPlatinum, 2 Pi], 1000]

✕ PolarPlot[1 + fPlatinum[t]^2/6, {t, 0, 200 2 Pi}, PlotStyle -> Directive[Opacity[0.4], RGBColor[0.36, 0.51, 0.71], Thickness[0.001]], PlotPoints -> 1000, Prolog -> {RGBColor[0.88, 0.61, 0.14], Disk[]}]

✕ differencesPlatinum = Differences[zerosPlatinum]; peaksPlatinum = maximaPositions[fPlatinum[x], x]; Histogram[differencesPlatinum, 1000, GridLines -> {peaksPlatinum, {}}]

We end this section with the square root of 11, which has an even more complicated zero-spacing distribution.

✕ \[Phi]Rhodium = (1 + Sqrt[11]); fRhodium[x_] := Cos[x] + Cos[\[Phi]Rhodium x] + Cos[\[Phi]Rhodium^2 x] zerosRhodium = findZeros[fRhodium, 10^5]; peaksRhodium = maximaPositions[fRhodium[x], x] Histogram[Differences[zerosRhodium], 1000, GridLines -> {peaksRhodium, {}}, PlotRange -> All]

Zero Distances between Nonsuccessive Zeros

So far, we have looked at the distribution of the distances between successive zeros. In this section, we will generalize to the distribution between any pair of zeros. While this seems like a more general problem, the equations describing the possible distances will easily allow us to determine the peak positions for successive zero distances.

If we consider not only the distance between consecutive zeros but also between more distant zeros, we get a more complicated distribution of zero distances.

✕ hgd = Histogram[ zeroDistancesAll = Select[(zerosGolden[[#2]] - zerosGolden[[#1]]) & @@@ Sort[Flatten[Table[{k, k + j}, {k, 20000}, {j, 1, 50}], 1]], # < 20 &], {0.001}]

All of the peaks are sums of the four peaks seen in the distances between consecutive zeros. But the reverse is not true—not all sums of consecutive zero distances show up as peaks. We identify the sums that are observed.

✕ observedPeaks = Sort[First /@ Select[Tally[ Round[zeroDistancesAll, 0.001]], Last[#] > 50 &]];

✕ calculatedPeaks = Flatten[Table[Union[{Total[First /@ #], Sort[Last /@ #]} & /@ Tuples[ Transpose[{{1.8136187, 1.04397, 1.49906, 3.01144}, Range[4]}], {j}]], {j, 1, 10}], 1];

✕ nf = Nearest[(#[[1]] -> #) & /@ calculatedPeaks];

✕ peakSums = DeleteDuplicates[Sort[nf[#][[1]] & /@ observedPeaks]]

Here is the left half of the last histogram shown together with the peak sum positions.

✕ Show[{hgd, ListPlot[Callout[ {#1, 20 + 12 #}, Row[#2]] & @@@ peakSums, Filling -> Axis, PlotStyle -> PointSize[0.008]]}]

The variable ranges from 0 to infinity in and so is not suited to be used as a parametrization variable to show large ranges. The three expressions , and are all in the interval . Above, we plotted the summands in these cosine terms; now we plot the triples , where is the distance between the and the zero. This gives some interesting-looking curves.

✕ zeroPairIndices = Sort[Flatten[Table[{k, k + j}, {k, 25000}, {j, 1, 30}], 1]];

✕ Graphics3D[{PointSize[0.001], RGBColor[0.36, 0.51, 0.71], Point[ {Cos[GoldenRatio zerosGolden[[#1]]], Cos[GoldenRatio^2 zerosGolden[[#1]]], zerosGolden[[#2]] - zerosGolden[[#1]]} & @@@ zeroPairIndices]}, BoxRatios -> {1, 1, 2}, Axes -> True, PlotRange -> {{-1, 1}, {-1, 1}, {0, 8}}]

In the case of , we just get a curve; in the case of general and , we obtain an intricate-looking surface—for instance, for .

✕ zerosEPi = findZeros[Function[x, Cos[x] + Cos[E x] + Cos[Pi x]], 4 10^4];

✕ Graphics3D[{PointSize[0.001], RGBColor[0.36, 0.51, 0.71], Point[ {Cos[ E zerosEPi[[#1]]], Cos[Pi zerosEPi[[#1]]], zerosEPi[[#2]] - zerosEPi[[#1]]} & @@@ zeroPairIndices]}, BoxRatios -> {1, 1, 2}, Axes -> True, PlotRange -> {{-1, 1}, {-1, 1}, {0, 8}}, ViewPoint -> {1.37, -3.08, 0.28}, AxesLabel -> {Cos[\[Alpha] Subscript[z, i]], Cos[\[Beta] Subscript[z, i]], Subscript[z, j] - Subscript[z, i]}]

Now, can we find a closed form of this surface? Turns out, we can. For generic and , we will have no algebraic relation between and , so these two expressions are natural independent variables. Assume that is a zero of and that another zero has distance to . Then the sum of the three cosines implies a relation between and and . We can calculate this relation by eliminating unwanted variables from and .

✕ f\[Alpha]\[Beta][x_] := Cos[x] + Cos[\[Alpha] x] + Cos[\[Beta] x]

✕ f\[Alpha]\[Beta][x0 + \[Delta]] // ExpandAll // TrigExpand

✕ c\[Alpha]c\[Beta]\[Delta]Equation = GroebnerBasis[ {f\[Alpha]\[Beta][x0], f\[Alpha]\[Beta][x0 + \[Delta]] // ExpandAll // TrigExpand, (* algebaric identities for Cos[...], Sin[...] *) Cos[x0]^2 + Sin[x0]^2 - 1, Cos[\[Alpha] x0]^2 + Sin[\[Alpha] x0]^2 - 1, Cos[\[Beta] x0]^2 + Sin[\[Beta] x0]^2 - 1}, {}, {Cos[x0], Sin[x0], Sin[\[Alpha] x0], Sin[\[Beta] x0]}, MonomialOrder -> EliminationOrder][[ 1]] /. {Cos[x0 \[Alpha]] -> c\[Alpha], Cos[x0 \[Beta]] -> c\[Beta]} // Factor;

The resulting polynomial (in , ) is pretty large, with more than 2,500 terms.

✕ {Exponent[c\[Alpha]c\[Beta]\[Delta]Equation, {c\[Alpha], c\[Beta]}], Length[c\[Alpha]c\[Beta]\[Delta]Equation]}

Here is a snippet of the resulting equation.

✕ Short[c\[Alpha]c\[Beta]\[Delta]Equation, 8]

The displayed curve and surface are special cases of this equation. Because of its size, we compile the equation for faster plotting.

✕ cf\[Alpha]\[Beta] = Compile[{\[Alpha], \[Beta], c\[Alpha], c\[Beta], \[Delta]}, Evaluate[c\[Alpha]c\[Beta]\[Delta]Equation], CompilationOptions -> {"ExpressionOptimization" -> True}]

We cut part of the surface open to get a better look at the inside of it. (Because of the complicated nature of the surface, we plot over a smaller vertical range compared to the above plot that used points for the zeros.)

✕ Module[{pp = 160, \[CurlyEpsilon] = 10^-1, data}, Monitor[ data = Table[ cf\[Alpha]\[Beta][GoldenRatio, GoldenRatio^2, c\[Alpha], c\[Beta], \[Delta]], {\[Delta], \[CurlyEpsilon], 4, (4 - \[CurlyEpsilon])/pp}, {c\[Beta], -1, 1, 2/pp}, {c\[Alpha], -1, 1, 2/pp}];, N[\[Delta]]]; ListContourPlot3D[data, DataRange -> {{-1, 1}, {-1, 1}, {\[CurlyEpsilon], 4}}, Contours -> {0}, RegionFunction -> (Not[#1 > 0 \[And] #2 < 0] &), MeshFunctions -> {Norm[{#1, #2}] &, #3 &}, BoundaryStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.004]], ViewPoint -> {2.49, -2.22, 0.53}, BoxRatios -> {1, 1, 2}, AxesLabel -> {Cos[\[Alpha] Subscript[z, i]], Cos[\[Beta] Subscript[z, i]], Subscript[z, j] - Subscript[z, i]}]]

In the case of and , we have the above algebraic equation of the Cayley surface at our disposal to remove the term.

✕ auxEq = Cos[x0]^2 + Cos[\[Alpha] x0]^2 + Cos[\[Beta] x0]^2 - (1 + 2 Cos[x0] Cos[\[Alpha] x0] Cos[\[Beta] x0]) /. Cos[x0] -> -Cos[\[Alpha] x0] - Cos[\[Beta] x0] /. {Cos[ x0 \[Alpha]] -> c\[Alpha], Cos[x0 \[Beta]] -> c\[Beta]}

(Note that the last equation is zero only for , , only at zeros and not for all values of . Eliminating the variable cβ gives us an even larger implicit equation for the distances of nonsuccessive zeros, with more than 74,000 terms.)

✕ c\[Alpha]\[Delta]Equation = Resultant[c\[Alpha]c\[Beta]\[Delta]Equation, auxEq, c\[Beta]];

✕ {Exponent[c\[Alpha]\[Delta]Equation, c\[Alpha]], Length[c\[Alpha]\[Delta]Equation]}

Luckily, this large equation factors in about 10 minutes into 4 much smaller equations, each having “only” 300 summands.

✕ (c\[Alpha]\[Delta]EquationFactored = Factor[c\[Alpha]\[Delta]Equation];) // Timing

✕ {Length[c\[Alpha]\[Delta]EquationFactored], Length /@ (List @@ c\[Alpha]\[Delta]EquationFactored)}

Here is a simplified form of the first factor. (For brevity of the resulting expression, we don’t yet substitute and , but will do so in a moment.)

✕ firstFactor = Collect[c\[Alpha]\[Delta]EquationFactored[[1]], c\[Alpha], FullSimplify]

✕ cpc\[Phi] = ContourPlot[ Evaluate[Block[{\[Alpha] = GoldenRatio, \[Beta] = GoldenRatio^2}, (# == 0) & /@ (List @@ c\[Alpha]\[Delta]EquationFactored)]], {c\[Alpha], -1, 1}, {\[Delta], 0.1, 6}, PlotPoints -> 60, AspectRatio -> 3/2] /. Tooltip[x_, _] :> x

We overlay a few thousand randomly selected distances from the 100k zeros of calculated above. The points all fall on the blue curve, meaning the equation firstFactor describes the position of the nonsuccessive zero distances. We indicate the positions of the maxima in the successive zero distances by horizontal gridlines. (As we used generic , in the derivation of the four functions of cαδEquationFactored , we could also use and obtain a similar image.)

✕ c\[Phi]Points = Module[{i, j}, Select[ Table[i = RandomInteger[{1, 99900}]; j = RandomInteger[{1, 10}]; {Cos[GoldenRatio zerosGolden[[i]]], zerosGolden[[i + j]] - zerosGolden[[i]]}, {50000}], #[[2]] < 6 &]];

✕ Show[{cpc\[Phi], Graphics[{Black, Opacity[0.2], PointSize[0.005], Point[c\[Phi]Points]}]}, GridLines -> {{}, {1.8136, 1.0439, 1.4990, 3.0114`} }]

We see that the peak positions of the successive zero distances (gridlines) are horizontal tangents on the curves. Intuitively, this is to be expected: at a horizontal tangent, many zero distances have approximately equal values, and so the singularities in the distribution form. This horizontal tangent observation now gives us a purely algebraic method to determine the peak positions of the zero distances. The condition of horizontal tangents on the given curve described by firstFactor gives two curves whose intersection points determine the peak positions we are looking for. Here are the curves for firstFactor and the curves that represent the conditions of horizontal tangents plotted together. The intersections of the blue and the yellow/brown curves are the relevant points.

✕ ContourPlot[ Evaluate[Block[{\[Alpha] = GoldenRatio, \[Beta] = GoldenRatio^2}, {firstFactor == 0, D[firstFactor == 0 /. \[Delta] -> \[Delta][c\[Alpha]], c\[Alpha] ] /. \[Delta]'[c\[Alpha]] -> 0 /. \[Delta][ c\[Alpha]] -> \[Delta]}]], {c\[Alpha], -1, 1}, {\[Delta], 0, 3.5}, PlotPoints -> 60, GridLines -> {{}, {1.8136, 1.044, 1.5, 3.011}}] /. Tooltip[x_, _] :> x

By eliminating the variable cα , we are left with a univariate equation in the zero spacing . But the resulting intermediate polynomials will be quite large (~877k terms!). So instead, we calculate a numerically high-precision approximation of the values of the horizontal tangents.

✕ firstFactorTangents = FindRoot[Evaluate[ Block[{\[Alpha] = GoldenRatio, \[Beta] = GoldenRatio^2}, {firstFactor == 0, D[firstFactor == 0 /. \[Delta] -> \[Delta][c\[Alpha]], c\[Alpha] ] /. \[Delta]'[c\[Alpha]] -> 0 /. \[Delta][ c\[Alpha]] -> \[Delta]}]], {c\[Alpha], #1}, {\[Delta], #2}, WorkingPrecision -> 50, Method -> "Newton"] & @@@ {{0.3, 1.8}, {-0.6, 1.04}, {0.3, 1.5}, {0.8, 3.01}}

The four peak position values agree perfectly with the previously calculated values.

✕ peakFunctions = Function[x, #1 Cos[x] + #2 Cos[GoldenRatio x] + #3 Cos[ GoldenRatio^2 x]] & @@@ {{1, 1, 1}, {-1, 1, 1}, {1, -1, 1}, {1, 1, -1}};

✕ #1[\[Delta]/2 /. #2] & @@@ Transpose[{peakFunctions, firstFactorTangents}]

As a function of the zero distance , the function firstFactor describes not only the distance between successive zeros, but also between distant zeros. In the following graphic, we connect the points of the zeros and for . The graphic shows how firstFactor does indeed describe the zero distances.

✕ Show[{ Graphics[{Thickness[0.001], Opacity[0.2], Table[{{RGBColor[0.368417, 0.506779, 0.709798], RGBColor[ 0.880722, 0.611041, 0.142051], RGBColor[ 0.560181, 0.691569, 0.194885], RGBColor[ 0.922526, 0.385626, 0.209179], RGBColor[ 0.528488, 0.470624, 0.701351], RGBColor[ 0.772079, 0.431554, 0.102387], RGBColor[ 0.363898, 0.618501, 0.782349], RGBColor[1, 0.75, 0], RGBColor[ 0.647624, 0.37816, 0.614037], RGBColor[ 0.571589, 0.586483, 0.], RGBColor[0.915, 0.3325, 0.2125], RGBColor[0.40082222609352647`, 0.5220066643438841, 0.85]}[[ k - 1]], Line[{ #[[-1]] - #1[[1]], Cos[GoldenRatio #[[1]]]} & /@ Take[Partition[zerosGolden, k, 1], 2000]]}, {k, 2, 7}]}, Axes -> True, PlotRange -> {{0, 10}, {-1, 1}}], ContourPlot[ Evaluate[ Block[{\[Alpha] = (1 + Sqrt[5])/ 2, \[Beta] = ((1 + Sqrt[5])/2)^2}, firstFactor == 0]], {\[Delta], 0.1, 10}, {c\[Alpha], -1, 1}, PlotPoints -> {80, 20}, ContourStyle -> Black] /. Tooltip[x_, _] :> x}, Frame -> True, Axes -> False, PlotRangeClipping -> True, AspectRatio -> 1/4]

The graphic also shows the string correlation between successive zeros that we saw in the previous pair correlation histogram.

✕ Histogram[ VectorAngle[#1 - #2, #3 - #2] & @@@ Partition[{ #[[-1]] - #1[[1]], Cos[GoldenRatio #[[1]]]} & /@ Partition[zerosGolden, 2, 1], 3, 1], 1000, {"Log", "Count"}]

Above, we plotted the zero distances over the complex plane. The expression firstFactor allows us to plot the value of over the complex plane. The next graphic shows the real part of over a part of the first quadrant.

✕ Module[{pp = 60, pts, c\[Alpha]Zeros}, (c\[Alpha]Zeros[\[Delta]_] := c\[Alpha] /. Solve[Block[{\[Alpha] = GoldenRatio, \[Beta] = GoldenRatio^2}, # == 0], c\[Alpha]]) &[firstFactor]; pts = Cases[ Flatten[Table[ N@{\[Delta]x, \[Delta]y, Re[#]} & /@ c\[Alpha]Zeros[N[\[Delta]x + I \[Delta]y]], {\[Delta]y, -0, 2, 2/pp}, {\[Delta]x, 0, 2, 2/pp}], 2], {_Real, _Real, _Real}]; Graphics3D[{RGBColor[0.36, 0.51, 0.71], Sphere[pts, 0.01]}, BoxRatios -> {1, 1, 3/2}, Axes -> True, PlotRange -> {All, All, 6 {-1, 1}}, AxesLabel -> {Re[\[Delta]], Im[\[Delta]], Re[c\[Alpha]]}, ViewPoint -> {1.98, -2.66, 0.65}, Method -> {"SpherePoints" -> 6}]]

Now we have all equations at hand to determine the two remaining peak positions of the function fBrass , which had the approximate values and . We use the implicit equation obeyed by and at zeros of .

✕ c\[Alpha]c\[Beta]Brass[c\[Alpha]_, c\[Beta]_] = Collect[1 + 2 c\[Alpha] - 3 c\[Alpha]^2 - 4 c\[Alpha]^3 + 4 c\[Alpha]^4 + 2 c\[Beta] + 2 c\[Alpha] c\[Beta] - 4 c\[Alpha]^2 c\[Beta] - 3 c\[Beta]^2 - 4 c\[Alpha] c\[Beta]^2 + 8 c\[Alpha]^3 c\[Beta]^2 - 4 c\[Beta]^3 + 8 c\[Alpha]^2 c\[Beta]^3 + 4 c\[Beta]^4, c\[Beta]]

Rather than eliminating the variable symbolically, for each cα , we numerically calculate all possible values for cβ and substitute these values into cαcβδEquation . For faster numerical evaluation, we compile the resulting expression.

✕ c\[Alpha]\[Delta]BrassCompiled = Compile[{{c\[Alpha], _Complex}, {\[Delta], _Complex}}, Evaluate[Block[{\[Alpha] = (1 + Sqrt[3])/ 2, \[Beta] = ((1 + Sqrt[3])/2)^2}, Block[{c\[Beta] = #}, c\[Alpha]c\[Beta]\[Delta]Equation] & /@ (c\[Beta] /. Solve[c\[Alpha]c\[Beta]Brass[c\[Alpha], c\[Beta]] == 0, c\[Beta]])]], CompilationOptions -> {"ExpressionOptimization" -> True}]

Calculating the values of cαcβBrass on a dense set of cα , δ grid points as blue points graphically the represents all possible zero distances as a function of cα . We overlay some of the previously calculated zero distances from above as orange points, the four identified peak positions as dark green lines and the two outstanding peak positions as red lines.

✕ Module[{ppc\[Alpha] = 400, pp\[Delta] = 600, data, cp}, Monitor[data = Table[c\[Alpha]\[Delta]BrassCompiled[ c\[Alpha], \[Delta]], {\[Delta], 0, 4, 4/pp\[Delta]}, {c\[Alpha], -1, 1, 2/ppc\[Alpha]}];, N[{\[Delta], c\[Alpha]}]]; cp = Show[ ListContourPlot[Re[#], Contours -> {0}, ContourShading -> None, DataRange -> {{-1, 1}, {0, 4}}, ContourStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.004]], PlotRange -> {{-1, 1}, {0.1, 4}}] & /@ Transpose[data, {2, 3, 1}]]; Show[{Graphics[{PointSize[0.004], Opacity[0.3], Orange, Point[{Cos[\[Phi]Brass #1], #2} & @@@ RandomSample[Transpose[{Most[zerosBrass], differencesBrass}], 10000]], Darker[Green], Opacity[1], Thickness[0.002], Line[{{-1, #}, {1, #}} & /@ {2.2299, 1.46396, 2.0722, 3.5761}], Red, Line[{{-1, #}, {1, #}} & /@ {1.098, 2.367}]}, AspectRatio -> 3/2, Frame -> True, PlotRangeClipping -> True, PlotRange -> {{-1, 1}, {0.1, 4}}], cp}]]

We clearly see that the two remaining peak positions also occur at points where the zero distance curve has a horizontal tangent. Expressing the condition of a horizontal tangent to the curve cαcβδEquation (with the constraint cαcβBrass ) numerically allows us to calculate high-precision values for the two remaining peak positions. (We also include the differentiated form of the constraint to have four equations in four variables.)

✕ c\[Alpha]c\[Beta]Brass = (-1 + x + 2 y^2)^2 + 4 (-1 + x - 2 x y^2) z^2 + 4 z^4 /. x -> -y - z /. {y -> c\[Alpha], z -> c\[Beta]} // Simplify

✕ horizontalTangentsBrass = Block[{\[Alpha] = \[Phi]Brass, \[Beta] = \[Phi]Brass^2}, {c\[Alpha]c\[Beta]\[Delta]Equation /. { c\[Beta] -> c\[Beta][c\[Alpha]], \[Delta] -> \[Delta][c\[Alpha]]}, D[c\[Alpha]c\[Beta]\[Delta]Equation /. { c\[Beta] -> c\[Beta][c\[Alpha]], \[Delta] -> \[Delta][c\[Alpha]]}, c\[Alpha]], c\[Alpha]c\[Beta]Brass /. c\[Beta] -> c\[Beta][c\[Alpha]], D[c\[Alpha]c\[Beta]Brass /. c\[Beta] -> c\[Beta][c\[Alpha]], c\[Alpha]]} /. {c\[Beta][c\[Alpha]] -> c\[Beta], Derivative[1][c\[Beta]][c\[Alpha]] -> c\[Beta]P, \[Delta][c\[Alpha]] -> \[Delta], \[Delta]'[ c\[Alpha]] -> 0}];

✕ Length[Expand[#]] & /@ horizontalTangentsBrass

✕ FindRoot[Evaluate[horizontalTangentsBrass], {c\[Alpha], -4305/100000}, {\[Delta], 10966/10000}, {c\[Beta], 97/100}, {c\[Beta]P, 78/100}, WorkingPrecision -> 30, PrecisionGoal -> 10, Method -> "Newton", MaxIterations -> 100]

✕ FindRoot[Evaluate[horizontalTangentsBrass], {c\[Alpha], 39/100}, {\[Delta], 2366/1000}, {c\[Beta], -38/100}, {c\[Beta]P, -71/100}, WorkingPrecision -> 30, PrecisionGoal -> 10, Method -> "Newton", MaxIterations -> 100]

And now the calculated peak positions of the zero distances agree perfectly with the numerically observed ones.

✕ Function[c, Histogram[Select[differencesBrass, c - 0.005 < # < c + 0.005 &], 100, GridLines -> {{1.0966482948, 2.3673493597}, {}}, PlotRangeClipping -> True, Method -> {"GridLinesInFront" -> True}]] /@ {1.0977, 2.366}

Power Laws near the Peaks of the Zero Distances?

Now that we have found equations for the positions of the singularities, how does the density of the zero distances behave near such a singularity? To find out, we select bins near the singularities and count the number of distances in these bins.

✕ singData = Table[Module[{sel}, sel = If[j == 2 || j == 3 || j == 4, Select[differencesGolden, Evaluate[peaksGolden[[j]] <= # <= peaksGolden[[j]] + 0.2] &], Select[differencesGolden, Evaluate[peaksGolden[[j]] - 0.2 <= # <= peaksGolden[[j]]] &]]; {Abs[#1 - peaksGolden[[j]]], #2} & @@@ Tally[Round[sel, 0.001]]], {j, 4}];

A log-log plot of the counts versus the distance to the singularities shows straight lines. This encourages us to conjecture that the behavior of the density near a singularity behaves as .

✕ ListLogLogPlot[singData, PlotRange -> {{5 10^-3, 1.8 10^-1}, All}, PlotLegends -> (Row[{"x", "\[TildeTilde]", NumberForm[#, 4]}] & /@ peaksGolden)]

The numerical value of depends on the concrete singularity, and seems to be in the range of . Numerical experiments with other sums of three trigonometric functions show a power law behavior near the singularities in general. The values of the exponents vary.

✕ Coefficient[ Fit[Log[Select[#, 10^-3 <= #[[1]] <= 10^-1 &]], {1, x}, x], x] & /@ singData

We can now also model what would happen if the phases in + + were random (but fulfilling the conditions for the envelopes). Here we do this for the local extrema at with function value –1 at this point.

✕ gGolden[x_, {\[CurlyPhi]1_, \[CurlyPhi]2_, \[CurlyPhi]3_}] := Cos[x + \[CurlyPhi]1] + Cos[GoldenRatio x + \[CurlyPhi]2] + Cos[GoldenRatio^2 x + \[CurlyPhi]3];

We select a random value for and calculate values for and so that , .

✕ phases = DeleteCases[ Table[Check[Block[{\[CurlyPhi]1 = RandomReal[{-Pi, Pi}]}, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3} /. FindRoot[ Evaluate[{gGolden[ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == -1, \ Derivative[1, {0, 0, 0}][gGolden][ 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] == 0}], {\[CurlyPhi]2, RandomReal[{-Pi, Pi}]}, {\[CurlyPhi]3, RandomReal[{-Pi, Pi}]}]], {}], {100000}] // Quiet, {}];

Due to the two constraint equations, the phases are correlated.

✕ Histogram3D[Mod[phases, 2 Pi][[All, {##}]], 100, AxesLabel -> {Subscript["\[CurlyPhi]", #1], Subscript[ "\[CurlyPhi]", #2]}] & @@@ \ \ {{1, 2}, {1, 3}, {2, 3}}

Two of the three remaining envelope zeros are clearly visible. (With the numerical root-finding method, one catches only a few curves that are near the green curve.) While the overall graphics look symmetric with respect to the line, the individual curves are not. In contrast to the behavior of , under the assumption of totally random phases between , and , the sums with random phases quickly take on values smaller than –3/2.

✕ Plot[Evaluate[gGolden[x, #] & /@ RandomSample[phases, 250]], {x, -2, 5}, PlotStyle -> Directive[Thickness[0.001], Opacity[0.2], Gray], Prolog -> {Plot[ Cos[x] - Cos[GoldenRatio x] - Cos[GoldenRatio^2 x], {x, -2, 5}, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Opacity[0.5]]][[1]], Plot[-Cos[x] + Cos[GoldenRatio x] - Cos[GoldenRatio^2 x], {x, -2, 5}, PlotStyle -> Directive[RGBColor[0.88, 0.61, 0.14], Opacity[0.5]]][[1]], Plot[-Cos[x] - Cos[GoldenRatio x] + Cos[GoldenRatio^2 x], {x, -2, 5}, PlotStyle -> Directive[ RGBColor[0.56, 0.69, 0.19], Opacity[0.5]]][[1]]}, GridLines -> {{}, {-1, 3}}, Epilog -> {Directive[Purple, PointSize[0.01]], Point[{#/2, 0} & /@ {1.04398, 1.49906, 3.01144}]}] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]

✕ gGoldenZeros = Quiet[x /. FindRoot[Evaluate[gGolden[x, #]], {x, 1/2}] & /@ phases];

The estimated in obtained from the zeros of these curves with random phases agree remarkably well with the two values from the above histograms. This suggests that, to a certain degree, the assumption of random local phases is justified.

✕ phasesmpp = With[{\[Xi] = peaksGolden[[2]]/2}, Select[gGoldenZeros, \[Xi] < # < \[Xi] + 0.1 &] - \[Xi]]; \[Gamma]2 = Coefficient[ Fit[N[Cases[Log[Tally[Round[phasesmpp, 0.001]]], {_Real, _}]], {1, x}, x], x]; Histogram[phasesmpp, 100, PlotLabel -> Row[{"\[Gamma]", "\[ThinSpace]\[TildeTilde]\[ThinSpace]", NumberForm[\[Gamma]2, 2]}]]

✕ phasespmp = With[{\[Xi] = peaksGolden[[3]]/ 2}, -(Select[gGoldenZeros, \[Xi] - 0.03 < # < \[Xi] &] - \[Xi])]; \[Gamma]3 = Coefficient[ Fit[N[Cases[Log[Tally[Round[phasespmp, 0.001]]], {_Real, _}]], {1, x}, x], x]; Histogram[phasespmp, 100, PlotLabel -> Row[{"\[Gamma]", "\[ThinSpace]\[TildeTilde]\[ThinSpace]", NumberForm[\[Gamma]3, 2]}]]

The appendix contains some calculations based on the computed zero density (assuming uniform distribution of phases) to see if the power law holds exactly or is approximate.

The Inflection Points of f ɸ (x)

For completeness, I repeat the distance investigations that were carried out for the zeros and extrema for inflection points.

The inflection points are the zeros of the second derivative.

✕ findInflectionPoints[f_, n_] := findZeros[f'', n]

✕ inflectionsGolden = Prepend[findInflectionPoints[fGolden, 100000], 0.];

The following plots are the direct equivalent to the case of extrema, so we skip commenting on them individually.

✕ Plot[fGolden[x], {x, 0, 10.3 Pi}, Epilog -> {Darker[Red], Point[{#, fGolden[#]} & /@ Take[ inflectionsGolden, 30]]}]

✕ Graphics[{PointSize[0.001], RGBColor[0.36, 0.51, 0.71], Opacity[0.2], Point[N[{#, fGolden[#]} & /@ inflectionsGolden]]}, AspectRatio -> 1/2, ImageSize -> 400, Frame -> True]

✕ Histogram[Mod[inflectionsGolden, 2 Pi], 200]

✕ Histogram[Differences[inflectionsGolden], 1000, PlotRange -> All]

✕ Histogram[Differences[inflectionsGolden, 2], 1000, PlotRange -> All]

✕ Histogram3D[{#2 - #1, #3 - #2} & @@@ Partition[inflectionsGolden, 3, 1], 100, PlotRange -> All]

✕ summandValuesI = {Cos[#], Cos[GoldenRatio #], Cos[GoldenRatio^2 #]} & /@ inflectionsGolden;

✕ Graphics3D[{PointSize[0.002], RGBColor[0.36, 0.51, 0.71], Point[Union@Round[summandValuesI, 0.01]]}, PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, Axes -> True, AxesLabel -> {Cos[x], Cos[GoldenRatio x], Cos[GoldenRatio^2 x]}]

Maximal Distance between Successive Zeros

Having found the positions of the peaks of the distribution of the distances of the zeros, another natural question to ask about the zero distribution is: what is the largest possible distance between two successive roots? The largest distance will occur in the following situation: starting at a zero, the function will increase or decrease, then have a first extremum, then a second and a third extremum, and then will have another zero. When the middle extremum barely touches the real axis, the distance between the two zeros will be largest. Here are some plots of the function around zeros that are the furthest apart. Note that while the curves look, at first glance, symmetric around x ≈ 1.6, the low maxima on the left side belongs to the curve with the high maxima on the right side and vice versa.

✕ Show[Plot[fGolden[# + t], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ Take[Sort[ Transpose[{differencesGolden, Most[zerosGolden]}]], -500][[All, 2]], PlotRange -> All, Frame -> True, GridLines -> {{0}, {-1, 3}}, ImageSize -> 400] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]

For the random phase model, I calculate the phases that make the middle extrema just touch the real axis.

✕ h[x_] = Cos[x + \[CurlyPhi]1] + Cos[\[Alpha] x + \[CurlyPhi]2] + Cos[\[Beta] x + \[CurlyPhi]3] ; \[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]1[] = (Solve[{h[x] == 0, D[h[x], x] == 0, D[h[x], \[CurlyPhi]1] == 0} /. x -> 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] /. _C -> 0 // Simplify); {Length[\[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]1[]], \[CurlyPhi]2\ \[CurlyPhi]3MaxD\[CurlyPhi]1[][[1]]}

Here is a plot of a solution with the middle extrema touching the real axis.

✕ touchingSolutions\[CurlyPhi]1[x_] = h[x] /. Select[ \[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]1[], Im[h[x] /. # /. {\[Alpha] -> GoldenRatio, \[Beta] -> GoldenRatio^2} /. x -> 5.] == 0 &] /. {\[Alpha] -> GoldenRatio, \[Beta] -> GoldenRatio^2};

Here is one of the solutions.

✕ touchingSolutions\[CurlyPhi]1[x][[1]] // FullSimplify

Here is a plot of a solution with the middle extrema touching the real axis. The four solutions each give the same maximal zero distance.

✕ Plot[touchingSolutions\[CurlyPhi]1[x], {x, -3, 3}]

The so-obtained maximal distance between two zeros agrees well with the observed value. The calculated maximum is slightly larger than the observed one. (The reverse situation would be bad.)

✕ (Min[x /. Solve[touchingSolutions\[CurlyPhi]1[x][[1]] == 0 \[And] 1/10 < x < 4, x, Reals]] - Max[x /. Solve[touchingSolutions\[CurlyPhi]1[x][[1]] == 0 \[And] -4 < x < -1/10, x, Reals]]) // N[#, 5] &

✕ Max[differencesGolden]

Finding the envelopes with respect to and does not give a larger result.

✕ \[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]2[] = (Solve[{h[x] == 0, D[h[x], x] == 0, D[h[x], \[CurlyPhi]2] == 0} /. x -> 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] /. _C -> 0 // Simplify); touchingSolutions\[CurlyPhi]2[x_] = With[{R = {\[Alpha] -> GoldenRatio, \[Beta] -> GoldenRatio^2}}, h[x] /. Select[ \[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]2[], Im[h[x] /. # /. R /. x -> 5.] == 0 &] /. R];

✕ \[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]3[] = (Solve[{h[x] == 0, D[h[x], x] == 0, D[h[x], \[CurlyPhi]3] == 0} /. x -> 0, {\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}] /. _C -> 0 // Simplify); touchingSolutions\[CurlyPhi]3[x_] = With[{R = {\[Alpha] -> GoldenRatio, \[Beta] -> GoldenRatio^2}}, h[x] /. Select[ \[CurlyPhi]2\[CurlyPhi]3MaxD\[CurlyPhi]3[], Im[h[x] /. # /. R /. x -> 5.] == 0 &] /. R] ;

✕ (Quiet[Min[x /. Solve[N[#] == 0 \[And] 0 < x < 4, x, Reals]] - Max[ x /. Solve[N[#] == 0 \[And] -4 < x < 0, x, Reals]]] & /@ #) & /@ {touchingSolutions\[CurlyPhi]1[x], touchingSolutions\[CurlyPhi]2[x], touchingSolutions\[CurlyPhi]3[x]}

Here are the three envelope curve families together with near 97,858.4.

✕ With[{x0 = 97858.38930}, Show[{Plot[fGolden[x0 + x], {x, -3, 3}, PlotStyle -> Directive[Thickness[0.01], Gray, Opacity[0.4]]], Plot[touchingSolutions\[CurlyPhi]1[x], {x, -3, 3}, PlotStyle -> RGBColor[0.88, 0.61, 0.14]], Plot[touchingSolutions\[CurlyPhi]2[x], {x, -3, 3}, PlotStyle -> RGBColor[0.36, 0.51, 0.71] ], Plot[touchingSolutions\[CurlyPhi]3[x], {x, -3, 3}, PlotStyle -> RGBColor[0.56, 0.69, 0.19] ]}, PlotRange -> All]]

Interestingly, the absolute minimum of the distance between two successive zeros in + + is slightly larger.

✕ rootDistance[{\[CurlyPhi]1_Real, \[CurlyPhi]2_Real, \ \[CurlyPhi]3_Real}] := (Min[Select[#, Positive]] - Max[Select[#, Negative]]) &[ N[x /. Quiet[ Solve[Cos[x + \[CurlyPhi]1] + Cos[GoldenRatio x + \[CurlyPhi]2] + Cos[GoldenRatio^2 x + \[CurlyPhi]3] == 0 \[And] -5 < x < 5, x]]]]

If one lets the three phases range over the domains , then one finds a slightly larger maximal distance between zeros.

✕ With[{pp = 24}, Monitor[Table[{{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}, rootDistance[ N@{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}]}, {\[CurlyPhi]1, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}, {\[CurlyPhi]2, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}, {\[CurlyPhi]3, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}];, N@{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}]]

✕ FindMaximum[ rootDistance[{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}], {{\ \[CurlyPhi]1, 23 Pi/12}, {\[CurlyPhi]2, 7 Pi/6}, {\[CurlyPhi]3, Pi/3}}] // Quiet

✕ FindMaximum[ rootDistance[{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}], {{\ \[CurlyPhi]1, 23 Pi/12}, {\[CurlyPhi]2, 7 Pi/6}, {\[CurlyPhi]3, Pi/3}}] // Quiet

This maximum distance is realized when a minimum (maximum) between two zeros barely touches the real axis.

✕ Plot[Cos[x + \[CurlyPhi]1] + Cos[GoldenRatio x + \[CurlyPhi]2] + Cos[GoldenRatio^2 x + \[CurlyPhi]3] /. %[[2]], {x, -3, 3}]

Calculating the first million zeros of does not yield a value larger than the previously calculated value of 3.334. This suggests that the last curve configuration is not realized along the real axis in .

Here is a list of record maxima slightly below the real axis, found by calculating the first few million extrema of . (This list was found using a version of the above function findZeros for the first derivative and only keeping a list of extremal shallow maxima.)

✕ shallowMaxima = {2.4987218706691797180876160929177, 33.704071215892008970482202830654, 98.293712744702474592256555931685, 443.88844400968497878246638591063, 3388.8528289090563871971140906274, 12846.898421437481029976313761124, 55352.647183638537564573877897525, 124704.59412664060098149321964301, 166634.14221987979291743435707392, 304761.83543691954802508678830822, 457972.87856640025996046175960675, 776157.81309371886102983541391071, 1220707.5925697200786171039302735};

The next plot shows how record maxima approach the real axis from below.

✕ ListLogLogPlot[{#, -fGolden[#]} & /@ shallowMaxima, AxesLabel -> {x, -Subscript[f, \[Phi]][x]}]

Here is a plot of at the least shallow value.

✕ {Plot[fGolden[shallowMaxima[[-1]] + x], {x, -2, 2}], Plot[fGolden[shallowMaxima[[-1]] + x], {x, -0.001, 0.001}]}

The root difference for this “near” zero is about 3.326058.

✕ Differences[ N[x /. Solve[ fGolden[Rationalize[shallowMaxima[[-1]], 0] + x] == 0 \[And] -3 < x < 3, x], 15]]

Interpolating the data from the last graphic, one gets for the value of the shallowest maxima up to .

The distances between the nearest roots to the right and left of these maxima seem to approach an upper bound below the 3.3340 value from above.

✕ Differences[ Function[\[CapitalDelta], x /. FindRoot[ Evaluate[fGolden[x] == 0], {x, # + \[CapitalDelta]}, WorkingPrecision -> 100, PrecisionGoal -> 30]] /@ {-2, 2}] & /@ \ shallowMaxima // Flatten // N[#, 10] &

So is it a special property of that the maximal root distance assuming arbitrary phases is not realized, or is it a more general property of all ? Numerically, experiments with various transcendental values of and suggest that generically the maximal possible root distance obtained from the envelope method agrees with the maximum observed.

Unfortunately, this time the algebraic formulation of the distances between the zeros does not provide the ultimate answer. The following graphic shows the four branches of the factored cαδEquation together with the observed points. No special structure is visible in the curves at the maximal-observed zero distances. The two “strands” of points correspond to the two curves shown in the first plot in this section.

✕ ContourPlot[ Evaluate[Block[{\[Alpha] = GoldenRatio, \[Beta] = GoldenRatio^2}, # == 0 & /@ (List @@ c\[Alpha]\[Delta]EquationFactored)]], {c\[Alpha], -1, 1}, {\[Delta], 3.25, 3.4}, PlotPoints -> 50, Epilog -> {Black, PointSize[0.01], Point[{Cos[GoldenRatio #1] , #2 - #1} & @@@ Select[Partition[zerosGolden, 2, 1], #[[2]] - #[[1]] > 3.25 &]]}, GridLines -> {{}, {1.813, 1.044, 1.5, 3.011, 3.334}}] /. Tooltip[x_, _] :> x

The new root that appears at the maximal root distance is in the algebraic formulation visible as vertical tangents in the curve. Plotting the vertical tangents at the positions of the largest zero distances observed and the observed zero distances shows this clearly.

✕ ContourPlot[ Evaluate[Block[{\[Alpha] = (1 + Sqrt[5])/ 2, \[Beta] = ((1 + Sqrt[5])/2)^2}, {firstFactor == 0}]], {c\[Alpha], -1, 1}, {\[Delta], 0, 3.6}, PlotPoints -> 50, Epilog -> {Black, PointSize[0.003], Point[{Cos[GoldenRatio #1] , #2 - #1} & @@@ RandomSample[Partition[zerosGolden, 2, 1], 15000]]}, GridLines -> {Mean /@ Split[Sort[ Cos[GoldenRatio #1] & @@@ Select[Partition[zerosGolden, 2, 1], #[[2]] - #[[1]] > 3.32 &]], Abs[#1 - #2] < 0.4 &], {} }] /. Tooltip[x_, _] :> x

Finding the points of vertical tangents and “lifting” the cα values of these points to the curve near the observed intersections gives the algebraic prediction for maximal root distances.

✕ verticalTangentsc\[Alpha]\[Delta] = FindRoot[Evaluate[ Block[{\[Alpha] = (1 + Sqrt[5])/ 2, \[Beta] = ((1 + Sqrt[5])/2)^2}, {firstFactor == 0, D[firstFactor == 0 /. c\[Alpha] -> c\[Alpha][\[Delta]], \[Delta]] /. c\[Alpha]'[\[Delta]] -> 0 /. c\[Alpha][\[Delta]] -> c\[Alpha]}]], {c\[Alpha], #1}, {\[Delta], #2}, WorkingPrecision -> 50, PrecisionGoal -> 20, MaxIterations -> 200] & @@@ {{-0.15, 1.63}, {0.68, 1.7}}

✕ (FindRoot[ Evaluate[ Block[{\[Alpha] = (1 + Sqrt[5])/ 2, \[Beta] = ((1 + Sqrt[5])/2)^2}, firstFactor == 0 /. #[[1]] ]], {\[Delta], 3.3}, WorkingPrecision -> 40, PrecisionGoal -> 20, MaxIterations -> 200] // N[#, 20] &) & /@ verticalTangentsc\[Alpha]\[Delta]

Both zero distances agree; this means the maximal root distance is about 3.32606.

Summary and Features Still to Look At

Through graphical explorations, statistical tests and numerical checks, we have been able to conjecture the answer to the original MathOverflow question: the positions of the peaks in the distribution of zero distances of are two times the positions of the smallest zeros of .

The concrete function exhibits an interesting mixture of generic and nongeneric properties due to the golden ratio factors.

Many other structures within the zeros, extrema and inflection points of the sum of three trigonometric functions could be investigated, as well as the relations between the zeros, extrema and other special points. Here are a few examples.

Special Points mod 2π

For instance, we could visualize the cosine arguments of the zeros, extrema and inflection points. Here are the argument values modulo .

✕ Graphics3D[{PointSize[0.002], RGBColor[0.36, 0.51, 0.71], Point[Mod[# {1, GoldenRatio, GoldenRatio^2}, 2 Pi] & /@ zerosGolden], RGBColor[0.88, 0.61, 0.14], Point[Mod[# {1, GoldenRatio, GoldenRatio^2}, 2 Pi] & /@ extremasGolden], RGBColor[0.56, 0.69, 0.19], Point[Mod[# {1, GoldenRatio, GoldenRatio^2}, 2 Pi] & /@ inflectionsGolden]}, PlotRange -> {{0, 2 Pi}, {0, 2 Pi}, {0, 2 Pi}}, Axes -> True] // Legended[#, LineLegend[{RGBColor[0.36, 0.51, 0.71], RGBColor[0.88, 0.61, 0.14], RGBColor[0.56, 0.69, 0.19]}, {"zeros", "extrema", "inflection points"}]] &

Correctness of the Random Phase Assumption

Or we could look in more detail at how faithful the zero distances are reproduced if one takes the random phases seriously and looks at all functions of the form + + . For a grid of , , values, we calculate the distance between the smallest positive and largest negative zero.

✕ zeroDistance0[ f_] := (Min[Select[#, Positive]] - Max[Select[#, Negative]]) &[ Table[x /. FindRoot[f, {x, x0}], {x0, -3, 3, 6/17}]] // Quiet

The resulting distribution of the zero distances looks quantitatively different from the zero distance distributions of . But because the peak positions arise from the envelopes, we see the same peak positions as in .

✕ With[{pp = 40}, Monitor[ zeroData = Table[zeroDistance0[ Cos[x + \[CurlyPhi]1] + Cos[GoldenRatio x + \[CurlyPhi]2] + Cos[GoldenRatio^2 x + \[CurlyPhi]3]], {\[CurlyPhi]3, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}, {\[CurlyPhi]2, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}, {\[CurlyPhi]1, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}];, N@{\[CurlyPhi]3, \[CurlyPhi]2, \[CurlyPhi]1}]]

✕ Histogram[Flatten[zeroData], 200, GridLines -> {peaksGolden, {}}]

Function Value Distribution in the Random Phase Assumption

Or we could model the function value distribution in the general case. The function value distribution for generic , is quite different from the one observed above for . Generically (see the examples in the postlude) it has a characteristic flat middle part in the interval . Assuming again that locally around the function looks like + + with and uniformly distributed , and , we can express the probability to see the function value as:

✕ P(y) = \!\( \*UnderoverscriptBox[\(\[Integral]\), \(0\), \(2 \[Pi]\)]\( \*UnderoverscriptBox[\(\[Integral]\), \(0\), \(2 \[Pi]\)]\( \*UnderoverscriptBox[\(\[Integral]\), \(0\), \(2 \[Pi]\)]\* TemplateBox[{RowBox[{"y", "-", RowBox[{"(", RowBox[{ RowBox[{"cos", "(", SubscriptBox["\[CurlyPhi]", "1"], ")"}], "+", RowBox[{"cos", "(", SubscriptBox["\[CurlyPhi]", "2"], ")"}], "+", RowBox[{"cos", "(", SubscriptBox["\[CurlyPhi]", "3"], ")"}]}], ")"}]}]}, "DiracDeltaSeq"] \[DifferentialD] \*SubscriptBox[\(\[CurlyPhi]\), \(1\)] \[DifferentialD] \*SubscriptBox[\(\[CurlyPhi]\), \(2\)] \[DifferentialD] \*SubscriptBox[\(\[CurlyPhi]\), \(3\)]\)\)\)

Integrating out the delta function, we obtain the following.

✕ P(y)~\!\( \*UnderoverscriptBox[\(\[Integral]\), \(0\), \(2 \[Pi]\)]\( \*UnderoverscriptBox[\(\[Integral]\), \(0\), \(2 \[Pi]\)] \*FractionBox[\(\[Theta](1 - \*SuperscriptBox[\((y - cos( \*SubscriptBox[\(\[CurlyPhi]\), \(1\)]) - cos( \*SubscriptBox[\(\[CurlyPhi]\), \(2\)]))\), \(2\)])\), SqrtBox[\(1 - \*SuperscriptBox[\((y - cos( \*SubscriptBox[\(\[CurlyPhi]\), \(1\)]) - cos( \*SubscriptBox[\(\[CurlyPhi]\), \(2\)]))\), \(2\)]\)]] \[DifferentialD] \*SubscriptBox[\(\[CurlyPhi]\), \(1\)] \[DifferentialD] \*SubscriptBox[\(\[CurlyPhi]\), \(2\)]\)\)

Carrying out the integral numerically gives exactly the function value distribution shown below.

✕ Pfv[y_] := NIntegrate[ Piecewise[{{1/(\[Pi] Sqrt[ 1 - (y - Cos[\[CurlyPhi]1] - Cos[\[CurlyPhi]2])^2]), 1 - (y - Cos[\[CurlyPhi]1] - Cos[\[CurlyPhi]2])^2 >= 0}}], {\[CurlyPhi]1, 0, 2 Pi}, {\[CurlyPhi]2, 0, 2 Pi}, PrecisionGoal -> 3]/(2 Pi)^2

✕ lpPfv = ListPlot[Join[Reverse[{-1, 1} # & /@ #], #] &[ Monitor[Table[{y, Pfv[y]}, {y, 0, 3, 0.05}], y]], Joined -> True, Filling -> Axis]

This distribution agrees quite well with the value distribution observed for .

✕ Show[{Histogram[Table[fPi[x], {x, 0, 1000, 0.001}], 100, "PDF"], lpPfv}]

Mean Zero Spacing in the Random Phase Assumption

The observed value of mean spacing between zeros is ≈1.78.

✕ meanGoldenZeroSpacing = Mean[differencesGolden]

Here is a plot showing how the average mean spacing evolves with an increasing number of zeros taken into account. The convergence is approximately proportional to .

✕ {ListPlot[Take[#, 100] &@Transpose[{ Most[zerosGolden], MapIndexed[#1/#2[[1]] &, Accumulate[differencesGolden]]}], PlotRange -> All, GridLines -> {{}, {meanGoldenZeroSpacing}}], Show[{ListLogLogPlot[Transpose[{ Most[zerosGolden], Abs[ MapIndexed[#1/#2[[1]] &, Accumulate[differencesGolden]] - meanGoldenZeroSpacing]}]], LogLogPlot[5 x^-0.88, {x, 1, 2 10^5}, PlotStyle -> Gray]}]}

If one uses the random phase approximation and average over all zero distances with at least one zero in the interval , and using a grid of values for the phases, one gets a value that agrees with the empirically observed spacing to less than 1%.

✕ zeroSpacings[{\[CurlyPhi]1_, \[CurlyPhi]2_, \[CurlyPhi]3_}] := Module[{sols, pos1, pos2}, sols = Quiet[ Sort[x /. Solve[N[Cos[x + \[CurlyPhi]1] + Cos[GoldenRatio x + \[CurlyPhi]2] + Cos[GoldenRatio^2 x + \[CurlyPhi]3]] == 0 \[And] -6 < x < 12, x]]]; pos1 = Max[Position[sols, _?(# < 0 &)]]; pos2 = Min[Position[sols, _?(# > 2 Pi &)]]; Differences[Take[sols, {pos1, pos2}]]]

✕ Module[{pp = 32, sp}, Monitor[ spacingArray = Table[ sp = zeroSpacings[{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}], {\ \[CurlyPhi]1, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}, {\[CurlyPhi]2, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}, {\[CurlyPhi]3, 0, 2 Pi (1 - 1/pp), 2 Pi/pp}];, {N@{\[CurlyPhi]1, \[CurlyPhi]2, \[CurlyPhi]3}, sp}]]

✕ Mean[Mean /@ Flatten[spacingArray, 2]]

Distribution of Zero-Nearest Extrema

We could also look in more detail at the distribution of the distances to the nearest zero from the extrema.

✕ nfZerosGolden = Nearest[zerosGolden]

✕ Histogram[Abs[# - nfZerosGolden[#][[1]]] & /@ extremasGolden, 1000]

Distribution of Areas under the Curve

We can look at the distribution of the (unsigned) areas under a curve between successive zeros.

✕ area[{a_, b_}] = Integrate[fGolden[x], {x, a, b}]

✕ Histogram[Abs[area /@ Partition[zerosGolden, 2, 1]], 500]

Zero Distance Distribution for Sums of Four Cosines

We should check if the natural generalizations of the conjectures’ peak positions hold. For instance, here is a sum of four cosine terms that uses the plastic constant.

✕ P = N@Root[-1 - # + #^3 &, 1] // ToRadicals; fPlastic[x_] = Cos[x] + Cos[P x] + Cos[P^2 x] + Cos[P^3 x];

We again calculate 100k zeros. We also calculate the conjectured peak positions and plot both together. One of the predicted peaks turns out to be an edge; the plastic constant is nongeneric for four cosine terms in the same sense as the golden ratio is for three terms.

✕ zerosPlastic = findZeros[fPlastic, 10^5];

✕ peaksPlastic = (2 x /. FindRoot[{Cos[x], Cos[P x] , Cos[P^2 x], Cos[P^3 x]}.#1, {x, #2}]) & @@@ {{{1, 1, 1, 1}, 1}, {{-1, 1, 1, 1}, 0.8}, {{1, -1, 1, 1}, 0.8}, {{1, 1, -1, 1}, 1}, {{1, 1, 1, -1}, 1.6}}

✕ Histogram[Differences[zerosPlastic], 1000, PlotRange -> All, GridLines -> {peaksPlastic, {}}]

The peak positions again agree perfectly. Plotting the function near the zeros with the peak distances shows similar envelopes as in the three-term case.

✕ With[{aux = Sort@Transpose[{Differences[zerosPlastic], Most[zerosPlastic]}]}, Function[v, Show[Plot[fPlastic[# + t], {t, -2, 0 + 5}, PlotRange -> All, PlotStyle -> Directive[RGBColor[0.36, 0.51, 0.71], Thickness[0.001], Opacity[0.3]]] & /@ Take[Sort[{Abs[#1 - v], #2} & @@@ aux], 100][[All, 2]], PlotLabel -> Row[{"x", "\[ThinSpace]\[Equal]\[ThinSpace]", NumberForm[v, 3]}], PlotRange -> All, Frame -> True, GridLines -> {{0}, {-2, 4}}] /. l_Line :> Mouseover[l, {Red, Thickness[0.002], Opacity[1], l}]] /@ peaksPlastic]

Visibility Graphs of the Extrema

A plot of the extrema invites us to make visibility graphs from the extrema. Remember that a visibility graph can be constructed from time series–like data by connecting all points that are “visible” to each other. The following graphic is mostly self-explanatory. We consider the vertical lines from the axis as visibility blockers (rather than the actual function graph). The function visibleQ determines if the point p1 is visible from the point p2 (and vice versa), where “visible” means that no line from the axis to any point between p1 and p2 blocks sight.

✕ visibleQ[{p1_, p2_}] := True visibleQ[{p1_, middle__, p2_}] := (And @@ (C[{p1, #, p2}] & /@ {middle})) /. C :> tf

✕ tf[{{t1_, v1_}, {s_, u_}, {t2_, v2_}}] := With[{U = v1 + (t1 - s)/(t1 - t2) (v2 - v1)}, If[u < 0, U > 0 || U < u, U > u || U < 0]]

✕ visibilityEdges[pts_] := Monitor[Table[ If[visibleQ[Take[pts, {i, j}]], i \[UndirectedEdge] j, {}], {i, Length[pts] - 1}, {j, i + 1, Length[pts]}], {i, j}] // Flatten

✕ extremasGoldenPoints[n_] := {#, fGolden[#]} & /@ Take[extremasGolden, n]

✕ With[{pts = {#, fGolden[#]} & /@ Take[extremasGolden, 10]}, Show[{Plot[fGolden[x], {x, 0, 11}, PlotStyle -> RGBColor[0.36, 0.51, 0.71]], ListPlot[pts, Filling -> Axis, FillingStyle -> RGBColor[0.88, 0.61, 0.14]], Graphics[{RGBColor[0.36, 0.51, 0.71], PointSize[0.02], MapIndexed[{Point[{#, fGolden[#]}], Black , Text[ #2[[1]], {#, fGolden[#] + 0.2}]} &, Take[ extremasGolden, 10]], Gray, Line[Map[pts[[#]] &, List @@@ visibilityEdges[extremasGoldenPoints[10]], {-1}]]}]}]]

Here is a larger graph. The maxima with are responsible for the far-distance edges (e.g. 1–179).

✕ visGr = Graph[ visibilityEdges[{#, fGolden[#]} & /@ Take[extremasGolden, 200]], VertexLabels -> "Name", EdgeLabels -> Placed["Name", Tooltip]]

The relevant information is contained in the degree distribution of the visibility graph.

✕ ListLogLogPlot[Tally[VertexDegree[visGr]]]

Products of Cosines

In addition to more summands, we could also look at products of cosine functions.

✕ f\[CapitalPi][x_] = Cos[x]*Product[If[IntegerQ[Sqrt[k]], 1, Cos[Sqrt[k] x]], {k, 2, 6}]

✕ Plot[f\[CapitalPi][x], {x, 0, 20}, PlotRange -> All]

✕ cosZeros[\[Alpha]_, n_] := Join[Table[(Pi/2 + 2 \[Pi] k)/\[Alpha], {k, 0, n}], Table[(-Pi/2 + 2 \[Pi] k)/\[Alpha], {k, n}]]

✕ zeros\[CapitalPi] = With[{n = 10000}, Sort[N@Flatten[{cosZeros[1, n], Table[If[IntegerQ[Sqrt[k]], {}, cosZeros[Sqrt[k], n]], {k, 2, 6}]}]]];

✕ Histogram3D[ Partition[Select[Differences[zeros\[CapitalPi]], # < 3 &], 2, 1], 100]

Modulo an increase in the degree of the polynomials involved, such products should also be amenable to the algebraic approach used above for the nonsuccessive zero distances.

Sums of Three Sines

If instead of we had used , the distribution of the zero differences would be much less interesting. (For generic , , there is no substantial difference in the zero-distance distribution between , but , is a nongeneric situation).

Calculating the distribution of spacings gives a mostly uniform distribution. The small visible structures on the right are numerical artifacts, and calculating the zeros with a higher precision makes most of them go away.

✕ fGoldenSin[x_] := Sin[x] + Sin[GoldenRatio x] + Sin[GoldenRatio^2 x] zerosGoldenSin = findZeros[fGoldenSin, 10^5]; Histogram[Differences[zerosGoldenSin], 1000, PlotRange -> All]

Visually, the spacing of + + does not seem to depend on .

✕ ContourPlot[ Cos[x + \[CurlyPhi]] + Cos[Pi x + \[CurlyPhi]] + Cos[Pi^2 x + \[CurlyPhi]] == 0, {x, -4 Pi, 4 Pi}, {\[CurlyPhi], 0, 2 Pi}, PlotPoints -> 120, AspectRatio -> 1/2, GridLines -> {{}, {Pi/2, Pi, 3/2 Pi}}]

Let’s end with the example mentioned in the introduction: + + .

✕ fSqrtSin[x_] := Sin[x] + Sin[Sqrt[2] x] + Sin[Sqrt[3] x] zerosSqrtSin = findZeros[fSqrtSin, 10^5];

The positions of the peaks are described by the formulas conjectured above.

✕ {peak1, peak2, peak3} = {2 x /. Solve[-Cos[x] + Cos[Sqrt[2] x] + Cos[Sqrt[3] x] == 0 \[And] 0 < x < 1, x][[1]], 2 x /. Solve[+Cos[x] + Cos[Sqrt[2] x] + Cos[Sqrt[3] x] == 0 \[And] 1 < x < 4, x][[1]], 2 x /. Solve[+Cos[x] + Cos[Sqrt[2] x] - Cos[Sqrt[3] x] == 0 \[And] 1 < x < 2, x][[1]]}

✕ Histogram[Differences[zerosSqrtSin], 1000, {"Log", "Count"}, PlotRange -> All, GridLines -> {{peak1, peak2, peak3}, {}}]

We note that one of the roots of is very near one of the observed peaks, but it does not describe the peak position.

✕ peak2Alt = x /. Solve[ Sin[x] + Sin[Sqrt[2] x] + Sin[Sqrt[3] x] == 0 \[And] 1 < x < 4, x][[1]]

✕ Histogram[Select[Differences[zerosSqrtSin], 2.2 < # < 2.35 &], 200, PlotRange -> All, GridLines -> {{{peak2, Red}, {peak2Alt, Blue}}, {}}]

We will end our graphical/numerical investigation of sums of three cosines for today.

The Density of the Zeros

A natural question is the following: can one find a symbolic representation of the density of the zero distances along the algebraic curve derived above? Based on the calculated zeros, we have the following empirical density along the curve.

✕ sdkc\[Alpha]\[Delta] = SmoothKernelDistribution[{Cos[GoldenRatio #], #2 - #1} & @@@ Take[Partition[zerosGolden, 2, 1], All], 0.03, PerformanceGoal -> "Quality"]

✕ Plot3D[PDF[ sdkc\[Alpha]\[Delta], {c\[Alpha], \[Delta]}], {c\[Alpha], -1.1, 1.1}, {\[Delta], 0, 4}, PlotRange -> All, Exclusions -> {}, PlotPoints -> 120, MeshFunctions -> {#3 &}, AxesLabel -> {c\[Alpha], \[Delta]}]

Quantitative Classification of the Curve Shapes around the Zeros

How can the behavior of a concrete sum of three cosines be classified around their zeros?

For the sums for which an algebraic relation between the three cosine terms exists, the zeros form curves in the plane. Along these curves, the shape of the function around their zeros changes smoothly. And at the self-intersection of the curve, shapes “split.” The next graphic plots the zeros for the above example . Mouse over the points to see shifted versions of near the zero under consideration.

✕ nfc\[Alpha]\[Delta]Brass = Nearest[c\[Alpha]\[Delta]ListBrass = ({Cos[\[Phi]Brass #], #2 - #1} \ -> #1) & @@@ Take[Partition[zerosBrass, 2, 1], All]];

✕ makeZeroPlotBrass[{c\[Alpha]_, \[Delta]_}] := Module[{pts = nfc\[Alpha]\[Delta]Brass[{c\[Alpha], \[Delta]}, {All, 0.01}]}, Plot[Evaluate[(Cos[# + x] + Cos[\[Phi]Brass (# + x)] + Cos[\[Phi]Brass^2 (# + x)]) & /@ RandomSample[pts, UpTo[10]]], {x, -4, 4}]]

✕ Graphics[{PointSize[0.002], Point[First[#]] & /@ RandomSample[c\[Alpha]\[Delta]ListBrass, 25000]}, Frame -> True, PlotRange -> {{-1, 1}, {0, 2.6}}] /. Point[l_] :> (Tooltip[Point[ l], Dynamic[makeZeroPlotBrass[l]]] )

For generic , with no algebraic relation between them, the zeros do not form curves in the plane. One could display the zeros in space where they form a surface (see the above example of the function ). Even after projection into the plane. While this gives point clouds, it is still instructive to see the possible curve shapes.

✕ zerosSqrt23 = findZeros[fSqrt, 100000];

✕ nfc\[Alpha]\[Delta]Sqrt23 = Nearest[c\[Alpha]\[Delta]ListSqrt23 = ({Cos[ Sqrt[2] #], #2 - #1} -> #1) & @@@ Take[Partition[zerosSqrt23, 2, 1], All]];

✕ makeZeroPlotSqrt23[{c\[Alpha]_, \[Delta]_}] := Module[{pts = nfc\[Alpha]\[Delta]Sqrt23[{c\[Alpha], \[Delta]}, {All, 0.01}]}, Plot[Evaluate[ fSqrt[# + x] & /@ RandomSample[pts, UpTo[10]]], {x, -4, 4}]]

✕ Graphics[{PointSize[0.002], Point[First[#]] & /@ RandomSample[c\[Alpha]\[Delta]ListSqrt23, 25000]}, Frame -> True, PlotRange -> {{-1, 1}, {0, 5}}, AspectRatio -> 2] /. Point[l_] :> (Tooltip[Point[ l], Dynamic[makeZeroPlotSqrt23[l]]] )

Appendix: Modeling the Power Law near the Peaks

Above, based on the observed distribution of the distances between consecutive zeros, a power-law like decay of the distribution was conjectured. In this appendix, we will give some numerical evidence for the power law based on the envelope that defines the smallest zero.

We remind ourselves that we are interested in the smallest zero of + + subject to the constraints , .

Before modeling the power law, let us have a closer look at the data, especially the minima with function values approximately equal to –1 and slope 0 and the following zero. We use the list zerosAndExtremaGolden from above to find such pairs.

✕ minimaAtMinus1 = Select[Cases[ Partition[zerosAndExtremaGolden, 2, 1], {{_, "extrema"}, {_, "zero"}}], (Abs[fGolden[#[[1, 1]]] + 1] < 0.02 \[And] Abs[fGolden'[#[[1, 1]]]] < 0.02) &];

These minima split visibly into two groups.

✕ Show[Plot[Evaluate[fGolden[# + t]], {t, -1, 3}, PlotStyle -> Directive[Thickness[0.001], RGBColor[0.36, 0.51, 0.71]]] & /@ RandomSample[minimaAtMinus1, 100][[All, 1, 1]], PlotRange -> 1 All]

We select the zeros smaller than 0.6.

✕ minimaAtMinus1B = Select[minimaAtMinus1, #[[2, 1]] - #[[1, 1]] < 0.6 &];

✕ minimumAndZeroDistances = {#[[1, 1]], #[[2, 1]] - #[[1, 1]]} & /@ minimaAtMinus1B;

At the minimum , we can locally write + + , which defines three phases , and in = + + .

✕ \[CurlyPhi]1\[CurlyPhi]2\[CurlyPhi]3Triples = ({Mod[-#1, 2 Pi], Mod[-GoldenRatio #, 1 2 Pi], Mod[-GoldenRatio^2 #1, 2 Pi], #2} & @@@ minimumAndZeroDistances) /. x_?(# > 5 &) :> (2 Pi - x);

Plotting and shows an approximate linear relationship.

✕ ListPlot[{{#1, #2} & @@@ \[CurlyPhi]1\[CurlyPhi]2\[CurlyPhi]3Triples, \ {#1, #3} & @@@ \[CurlyPhi]1\[CurlyPhi]2\[CurlyPhi]3Triples}, PlotLegends -> {Subscript["\[CurlyPhi]", 2], Subscript["\[CurlyPhi]", 3]}, AxesLabel -> {Subscript["\[CurlyPhi]", 1]}]

Here are the triples observe