周期的に変化する関数の周期を知りたい

　計算の結果が周期関数的なふるまいをすることがあり、その周期の具体的な値を求めたいことがあります。そのため、任意の計算結果から周期を求めるアルゴリズムを自作してみました。適切なデータ間隔であれば良い精度で正しい値を求められそうです。

周期を求めるためのアイデアの説明

　周期関数$y = f(x)$が閾値$y = h$の上下で振動しているとします。基本的には、$h$を下から上、または、上から下に横切る時から$h = f(x_h)$となる$x_h$を3つ以上記録し、それらから周期を求めることがゴールです。しかし、不連続なデータの集まりである計算結果から正確な$x_h$をどのように求めたらよいでしょうか。ここで$h
\simeq f(x)$の$f(x)$が直線的なふるまいとなると仮定して、三角形の相似を用いることとしました。(下図)

このように考えれば、$x_h$はデータ点$(x_i,y_i), (x_{i+1},y_{i+1})$から、


x_h = x_i + (x_{i+1}-x_i) \times \frac{|y_i - h|}{|y_i - h|+|y_{i+1}-h|}

と求めることができるでしょう。

プログラムの作製

　プログラムは以下のようになりました。$x_h$が3つ未満であればエラー文を返すのと、それ以上であれば各点から周期を計算して平均化し返すようにしました。使用したパッケージはnumpyだけです。

'''
周期的なふるまいの計算結果list_xとlist_yと
閾値thresholdを与えることで周期を返す関数
'''
import numpy as np

def period_checker(list_x,list_y,threshold):
    #xhを格納するためのリスト
    list_x_zero = np.array([])

    #xhを計算するためのデータ点を選別。thresholdから±0.5のyの値のみ。
    list_y_check = list_y[abs(list_y - threshold) < 0.5]-threshold
    list_x_check = list_x[abs(list_y - threshold) < 0.5]

    #選別したデータから、thresholdを横切る時タイミングを確認しxhを計算・記録。
    for i in range(list_y_check.shape[0]-1):
        if list_y_check[i]*list_y_check[i+1]<0:
            ratio = abs(list_y_check[i])/(abs(list_y_check[i])+abs(list_y_check[i+1]))
            list_x_zero = np.append(list_x_zero, list_x_check[i]+ratio*(list_x_check[i+1]-list_x_check[i]))

    #記録されたxhが3未満であればエラー文を返す。
    if list_x_zero.shape[0] < 3:
        return print('please input valid interval \n')
    #xhが3以上であれば、それらより周期を求め平均化したのち、求めたxhと周期を表示
    else:
        periodT = 0
        Num_of_period = 0
        for j in range(list_x_zero.shape[0]-2):
            periodT += list_x_zero[j+2] - list_x_zero[j]
            Num_of_period += 1
        periodT = periodT/Num_of_period
        return print('x(y = {:.4f}) \n'.format(threhold),list_x_zero,'\n period : {:.18f} \n'.format(periodT))

プログラムのテスト

　周期関数の代表であるsin関数とcos関数でプログラムをテストしてみました。これらの周期の厳密解は$2 \pi$なので比較することでプログラムの正確性が測れるでしょう。
　例えば、$x$の範囲を$0$~$17$とし、その間を$5000$分割した時のsin関数とcos関数の周期を以下のように計算してみました。ついでにnumpyの組み込み関数を用いて$2 \pi$の値も表示しました。

In:
h = 0
list_angle = np.linspace(0,17,5000)
list_sin = np.sin(list_angle)+h
list_cos = np.cos(list_angle)+h

period_checker(list_angle,list_sin,h)
period_checker(list_angle,list_cos,h)

print(2*np.pi)

out:
x(y = 0.0000) 
 [ 3.14159265  6.28318531  9.42477796 12.56637061 15.70796327] 
 period : 6.283185307836515854 

x(y = 0.0000) 
 [ 1.57079633  4.71238898  7.85398163 10.99557429 14.13716694] 
 period : 6.283185307875393200 

6.283185307179586

　$x_h$がそれぞれ5つずつ記録されており、計算された周期が出力されました。また、組み込み関数の値と周期の値を比べると、それぞれの計算結果は$10^{-9}$の桁数まで一致していることが分かります。より分割数を増やすことでより高い精度の答えを求めることもできます。正確な結果を求めるには$h
\simeq f(x)$付近での十分なデータ数が必要なようです。

#まとめ

　周期関数の周期を求めるプログラムを作製してみました。十分なデータ点があればよい精度が求められそうです。よい補間法と組み合わせることで少ないデータ点から周期を求めることもできるかもしれませんね。それは、また思いついたときに……。