Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
2  
FourierTransform? –  Verbeia May 21 '12 at 6:20
4  
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
3  
You might find this paper of some use. Symmetry, Anti-Symmetry, Periodicity, are all discussed. cs.berkeley.edu/~fateman/papers/symmetry.pdf –  Richard Fateman May 6 at 22:13

4 Answers 4

up vote 25 down vote accepted

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`*"]
share|improve this answer
6  
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. –  J. M. May 21 '12 at 8:39
3  
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$. –  J. M. 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
1  
@Anixx: nothing wrong with that... it's when you're dealing with products of trigonometric functions that things get slightly hairy. –  J. M. May 21 '12 at 14:59
1  
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]

enter image description here

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]

enter image description here

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

share|improve this answer
6  
@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. –  J. M. 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. –  J. M. 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

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.

share|improve this answer

My earlier post, directing people to a paper, seems to have been relegated to a comment on some other marginally-related post, so I will try to explain in-line enough of why the answer to the question is related to symmetry, and what one might in fact use a computer algebra system to do about it. I suppose it is possible to copy the whole paper into stackexchange, but is that really useful?

Here's an except from the paper..

Consider symmetry as invariance under a group of transformations. Thus, given a transformation $T$ that maps some domain $R$ to itself, we say that $f: ~R\rightarrow R$ is symmetric under $T$ if and only if $f = T[f]$. Analyzing the set of all possible $T$ in full generality is neither practical nor particularly helpful. Instead, let us focus on transformations $T$ that can be expressed in a particular, but still interesting, form. Given $f(x)$, define $T[f](x)$ as a combination of linear transforms around the function $f$ as follows: $$ T_{a, b, c, d}[f](x) = a * f( b * x + c) + d; $$ Such a family of transformations is parameterized by 4 constants. For our current discussion you may assume they are real numbers $\{a,~b,~c,~d\}$, though this could be generalized. The family of transformations is fairly powerful; in particular, it allows us to express the forms previously mentioned plus other traditional notions of symmetry. Thus the notion that a function $f$ is even simply means that it is invariant under $T_{1, -1, 0, 0}$. Similarly, the notion that a function is odd simply means that it is invariant under $T_{-1, -1, 0, 0}$. This family of transformations also allows us to express periodicity, since a function periodic with period $p$ will be invariant under $T_{1, 0, p, 0}$.*

..........

A solution to this problem can be expressed as a call to Reduce in Mathematica, though exactly what Reduce will do depends on the version of Mathematica, past, present, and future. The paper was written in 2004. http://www.cs.berkeley.edu/~fateman/papers/symmetry.pdf

and no, it does not mention Fourier series, which is an alternative but rather unrelated approach to determining periodicity.

share|improve this answer
2  
Hello Richard, I was the one that "relegated" your previous post to a comment, but not to a "marginally-related post"... rather, I converted it to a comment on the very post that you answered (so surely must be related, right?)! I was responding to flags from other knowledgable users that your post was insufficient to stand alone as an answer as per the site's standards and was more appropriate as a comment. The excerpt is a good start, but as Artes suggested earlier, some Mathematica code demonstrating it would be helpful to the user and the reader. –  rm -rf May 7 at 21:13
1  
I'll admit, I know nothing of this topic to judge whether your statement that "A solution to this problem can be expressed as a call to Reduce in Mathematica,..." is sufficient for the reader, but perhaps someone else can weigh in on that. –  rm -rf May 7 at 21:14
1  
I suppose that one could examine the (rather complicated) output of Reduce[Sin[x + a] == Sin[x], x] and discover an a==2Pi*C[1] near the end. It's really incumbent upon Richard to provide the example he suggests. –  Mark McClure May 7 at 21:26
3  
Generally, it's fine to provide pointers to papers but, if that's all and "answer" does, then a comment is really more appropriate. Also, the comment was placed under the original question. If that's genuinely a "marginally-related post", then why are you posting an answer at all? –  Mark McClure May 7 at 21:40
1  
@RichardFateman There is an FAQ. Among other things, you'll learn that you can create new tags after you have 300 rep and how to format your posts. –  Mark McClure May 9 at 9:54

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.