How to find the period of an arbitrary mathematical function?

Question

Is there a function to find the period of an arbitrary (possibly complex) function in Mathematica?

Answer 1

You can check out this one. I don't know how well it works

Periodic`PeriodicFunctionPeriod[E^(I 2 Pi t) + Cos[3/9 Pi  t], t]

6

Perhaps you are also interested in the other functions in that context. Check

Names["Periodic`*"]

Answer 2

Here is a numeric approximation method that can be useful when no analytic information is known. I will illustrate with the function WeierstrassPPrime[t, {2, 3}] that was mentioned in a comment to one response.

We begin by taking random steps, and sampling the function at those steps (I'll explain the random step size presently). We then plot the coordinates thus obtained, to do a visual check that we have at least a few periods.

SeedRandom[111222333];
xvals = Accumulate[RandomReal[.01, 3000]];
xyvals = Chop[Table[{t, WeierstrassPPrime[t, {2, 3}]}, {t, xvals}]];

ListPlot[xyvals]

So far so good. We now order these by their y coordinates, and subtract x values that correspond to neighboring y values.

sortd = xyvals[[Ordering[xyvals[[All, 2]]]]];
xdiffs = Differences[sortd[[All, 1]]];

Of course the neighboring y values might not be terribly close. I will use the following heuristic to "order" them in terms of which corresponding x differences I want to use. We take each y difference, and divide by the sqrt of sum of squares of the two y values from which it arises. We regard x differences as being more useful when that ratio is smaller. The idea here is that if the y difference is zero then, for this function (which is monotonic on a period), we are "close" to an integral number of periods apart in corresponding x values.

ydiffs = Abs[Differences[sortd[[All, 2]]]];
ynorms = Map[Norm, Partition[sortd[[All, 2]], 2, 1]];
weighteddiffs = ydiffs/ynorms;

So we now sort according to the heuristic described above.

sortdxdiffs = Abs[xdiffs[[Ordering[weighteddiffs]]]];

Next we discard those x differences that are (approximately) zero periods apart.

usediffs1 = 
  Select[sortdxdiffs[[1 ;; 500]], # >= Max[sortdxdiffs]/10 &];

From these we can cull the differences that are approximately one (as opposed to more) period apart.

usediffs = Select[usediffs1, # <= 1.5*Min[usediffs1] &];

Finally we average these. This gives a fairly good approximation to a period.

In[237]:= Mean[Abs[usediffs]]

Out[237]= 2.39443

So we have an approximate period of 2.3944 or so.

I will now mention why we did not take regular spacing: had we done so, we would require at least an order of magnitude more points to get an approximation this good; best we might get from points spaced at intervals of .001 is an approximation to that many places. While I do not claim our averaging approximation to be perfect, I will surmise it is good to at least three or so places.

One might use the differences that are multiple periods apart to get a possibly better estimate via weighted averaging (the idea being that more samples giving a lower variance). I'm just showing a simple computation requiring minimal work to get an average.

There are (at least) two problems with this approach. One is that we used real data, whereas the request was to handle complex periodic functions. It is not too difficult to extend to complex. One might work with real and imaginary parts separately, for example, and either average or use the lesser if one is roughly an integral multiple of the other.

A bigger problem is that we relied on monotonicity within periods in order to assert that coordinates with approximately equal y values had x values that differ by roughly an integer number of periods. Provided there are no flat parts (e.g. from a Haar wavelet), one can often deduce which x differences correspond to period separations by clumping them. The ones that are an integral number of periods apart will tend to be more common than differences that come from close y values where corresponding x pairs were not separated by an integral number of periods. This idea at least works for functions that are, to crude approximation, roughly sinusoidal. I'm sure there are more pathological examples where it will fail. In such cases some plotting and human intervention might be needed in order to figure out which x differences correspond to actual period separations and which come from inadvertent closeness of corresponding y values.

This process can be much harder for experimental data. For an example of such, see this cepheid demonstration or this related pulsar star period demonstration How to handle such a problem is another topic for another day.

--- edit ---

Here is an example of how one might deal with nonmonotonicity. I will use a function mentioned in the comments.

func[x_] := WeierstrassP[x, {2, 3}]

Again we will take a large number of samples. Some experimenting indicated what would be a good set. I'll show the plot as well.

SeedRandom[111222333];
xvals = Accumulate[RandomReal[.001, 10000]];
xyvals = Chop[Table[{t, Chop[func[t]]}, {t, xvals}]];

ListPlot[xyvals]

Again we create differences of x coordinates that correspond to neighboring y values. This time I'll be less fussy and use all of them, excluding ones that are quite small

sortd = xyvals[[Ordering[xyvals[[All, 2]]]]];
xdiffs = Differences[sortd[[All, 1]]];
sortdxdiffs = Sort[Select[Abs[xdiffs], # >= .1 &]];

Here is where the clumping happens. We will split this list, grouping pairs that are within .001 of one another. I'll explain why in a moment.

collected = Split[sortdxdiffs, #2 - #1 <= .001 &];


In[211]:= collected // Length

Out[211]= 1606

Now we look at the lengths of these groupings. One will stand out as being huge relative to the rest. This is the "main clump" (apologies for use of technical jargon).

lens = Map[Length, collected];
Table[RankedMax[lens, k], {k, 1, 5}]

Out[231]= {2800, 169, 45, 30, 29}

We take this clump and average its members. This gives our estimate for the period.

mainclump = First[Select[collected, Length[#] > 2000 &]];
Mean[mainclump]

Out[233]= 2.39443

So why does this grouping tactic work? The idea is that for a curve like this, it is roughly a 50-50 proposition whether y neighbors will belong to x neighbors that are roughly a period apart vs. on the flip side of the same period. In the latter case the x distances are spread out over a wide range, whereas the (approximately) period-separated x distances will aggregate near the period. As nearly half fall into that class, we go for that large set. Hence the terminology of "clumping".

The 50-50 above is loose talk. Since y neighbors can come from actual x neighbors, and since we do not have as many full periods as mirrored pairs of half periods, it is more like a bit under one in three that we get full period separation of x values for neighboring y pairs. But so what? We still get our clump. Had we used more periods, tehre would be other clumps of integral period separation, and we could use those as well in the averaging process, of course weighting by the appropriate integer for each.

--- end edit ---

Answer 3

I found another method:

Period[f_, t_] := 
 D[First[a /. 
    Solve[Refine[
      Reduce[{Function[t, f][x + a] == Function[t, f][x]}, a], 
      C[1] \[Element] Integers], a]], C[1]]
Period[Cos[3 x], x]

But it seems it works only for basic functions.