주기적으로 변화하는 함수의 주기를 알고 싶습니다.

계산의 결과가 주기 함수적인 행동을 할 수 있고, 그 주기의 구체적인 값을 구하고 싶은 경우가 있습니다. 그 때문에, 임의의 계산 결과로부터 주기를 구하는 알고리즘을 자작해 보았습니다. 적절한 데이터 간격이면 좋은 정밀도로 올바른 값을 구할 수 있을 것 같습니다.

주기를 찾는 아이디어 설명

주기 함수 $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)$ 부근에서 충분한 데이터 수가 필요한 것 같습니다.

요약

주기 함수의 주기를 구하는 프로그램을 제작해 보았습니다. 충분한 데이터 점이 있으면 좋은 정밀도가 요구될 것 같습니다. 좋은 보간법과 조합하는 것으로 적은 데이터점으로부터 주기를 구할 수도 있을지도 모르겠네요. 그것은, 또 떠올랐을 때에…….