How to find the period of an arbitrary mathematical function?

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

-
–  Verbeia May 21 '12 at 6:20
Please see here: stats.stackexchange.com/a/16234/4764 An important question is: do you have the analytic form of the function, or do you have have it sampled in a number of points (i.e. you have some numerical data)? If you have numerical data, are the data points equally spaces in time or not? –  Szabolcs May 21 '12 at 7:35

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

PeriodicPeriodicFunctionPeriod[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*"]

-
Works nicely for Sech[t], t-Floor[t], (-1)^Floor[t], and (surprisingly) JacobiDN[t, 1/3] and EllipticTheta[3, t, 1/2], but fails for WeierstrassPPrime[t, {2, 3}]. Clearly the thing is still under development, but it's getting there. –  Ｊ. Ｍ. May 21 '12 at 8:39
Another possible weakness of the function is that it does not always return the minimal period. Sin[3t] Sin[5t] is $\pi$-periodic, but the function returns a period of $2\pi$. –  Ｊ. Ｍ. May 21 '12 at 9:30
Wolfram Alpha has the function "period". To what it is rendered when Mathematica kernel is called? –  Anixx May 21 '12 at 13:22
@Anixx: nothing wrong with that... it's when you're dealing with products of trigonometric functions that things get slightly hairy. –  Ｊ. Ｍ. May 21 '12 at 14:59
If you first apply any of TrigReduce, TrigFactor, or TrigExpand to Sin[3t]Sin[5t], you get the correct result. Strangely, if you apply TrigToExp, you get $-\pi$. –  Mark McClure May 21 '12 at 17:57

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 ---

-
@Mr.Wizard Real programmers don't define all their variables. –  Daniel Lichtblau May 21 '12 at 20:54
Nice. Comparing the result of the approximate method with the result of N[2 WeierstrassHalfPeriods[{2, 3}]][[1]] // Chop shows that the algorithm returns pretty good estimates. I wonder if this can be extended to find the other period? However, if we look at the function WeierstrassP instead of WeierstrassPPrime, the method here does not do so well, because as already said, the current method relies on monotonicity within a single period. –  Ｊ. Ｍ. May 22 '12 at 2:20
For another example of a function where the method in this answer does not perform so well, try the Dixon elliptic function 6 WeierstrassP[x, {0, 1/27}]/(1 - 3 WeierstrassPPrime[x, {0, 1/27}]), whose real period is Beta[1/3, 1/3]/3; it is monotonic within a single period, but has a flat section. –  Ｊ. Ｍ. May 22 '12 at 2:51
I tried that example. It appears to have a period of Beta[1/3, 1/3], or around 5.29992. My method gave 5.2998. If I make use of more of the data I get 5.999. These are, suffice it to say, quite close to the mark. –  Daniel Lichtblau May 22 '12 at 15:33
@J.M. WeierstrassP is a harder case. I'll show how to handle that in an edit. It will show what I had meant by "clumping". –  Daniel Lichtblau May 22 '12 at 15:51
show 1 more comment

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.

-