フーリエ級数を用いて正多角形を描く（その２）

前回の続き、フーリエ級数を用いた正多角形を描くアルゴリズムについて。

正n角形を描く：
正n角形を描くには、

y = Σ(1/C² * sin(C * θ))
θ = 0〜2π

Cは無限級数であり、例えば正方形（n=4）を描くには、

C = [..., -3*n+1, -2*n+1, -1*n+1, 0*n+1, 1*n+1, 2*n+1, ...]
C = [..., -11, -7, -3, 1, 5, 9, ...]

になり、1を基準としてn=4ずつ加算、あるいは負の方向であれば4ずつ減算した級数。
あるいは、k%n==1　(k=-k, ..., -3, -2, -1, 0, 1, 2, 3, ... k)のときのkの値。
仮に級数を6つに限定すれば、

C = [-11, -7, -3, 1, 5, 9]

となります。そして各波形を

y = 1/C² * sin(C * θ)

とすれば以下のようになります。

横軸はθ、縦軸の値はCの値。C=9やC=-11のときはほぼ直線に見えますが、僅かながら振幅しています。

さらに、これらの波形を足し合わせると、

y = Σ(1/C² * sin(C * θ)),　　C = [-11, -7, -3, 1, 5, 9]

以下のような波形になります。

この波形を横軸x、縦軸yに置き換えると、

x = Σ(1/C² * cos(C * θ)) / Σ(1/C²), C = [-11, -7, -3, 1, 5, 9]

y = Σ(1/C² * sin(C * θ)) / Σ(1/C²), C = [-11, -7, -3, 1, 5, 9]

正方形に近似した軌跡になります。まだ正方形としては不完全なかたちですが、Cの個数を増やすとより正方形に近似していきます。

仮にCの個数を20個にすると、

C = [-39, -35, -31, -27, -23, -19, -15, -11, -7, -3, 1, 5, 9, 13, 17, 21, 25, 29, 33, 37]

になります。

sin()、cos()の代わりにexp()をつかって複素数であらわせば、

V = Σ(1/C² * exp(i * C * θ)) / Σ(1/C²),

になります。iは虚数、Vは複素数でVの実部がX座標、Vの虚部がY座標になります。

これをコーディングすると以下。

正方形の場合（n=4）：

import numpy as np
import matplotlib.pyplot as plt

n = 4
p = n * 10
theta = np.linspace(0, 2*np.pi, p+1)
F = 20
C = np.array([(1 + f * n) for f in range(-F//2, F//2)])
print('C =', C)
V = sum([1 / c**2 * np.exp(1j * c * theta) for c in C]) / sum(1 / C**2)

X = V.real
Y = V.imag

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.axis('equal')
ax.plot(X, Y)
ax.scatter(X, Y, s=10, c='red')

nは正n角形、pはtheta：0〜2πにおける分解能。Pythonでは虚数iは1j。
Cは無限フーリエ級数の係数ですが、20個（F=20）にしておきます。
V.realはVの実部でX座標、V.imagはVの虚部でY座標。

n=4（正方形）、　p=n*10=40、F=20で描くと以下。

赤い点は設定した分解能pの数だけありますが、より正確な図形を描く場合はFの数を増やします。F=20であればかなり正方形に近似します。

アニメーションで表現：

n=4（正方形）、　p=n*10=40、F=20

コード：

中心角thetaを分解能pで分けアニメーションにしています。

nを変えれば、正三角形（n=3）、正五角形（n=5）などになります。また、pの分解能をあげればより細かいフレームレートになります。あるいはFで級数の数を増やせば、より正確な図形の軌跡を描きます。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML

n = 4
p = n * 10
theta = np.linspace(0, 2*np.pi, p+1)
F = 20
C = np.array([(1 + f * n) for f in range(-F//2, F//2)])
print('C =', C)

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.axis('equal')

ims = []
X = []
Y = []
for t in theta:
    V = sum([1 / c**2 * np.exp(1j * c * t) for c in C]) / sum(1 / C**2)
    x = V.real
    y = V.imag
    X.append(x)
    Y.append(y)
    P1 = ax.plot(X, Y, c='tab:blue')
    P2 = ax.scatter(x, y, s=40, c='red')
    ims.append(P1 + [P2])
    
ani = animation.ArtistAnimation(fig, ims, interval=150)
ani.save('ani0.gif', writer="imagemagick")
#ani.save('ani0.mp4', writer="ffmpeg")
plt.close()
HTML(ani.to_jshtml())

Epicycleで表現：

正n角形の軌跡は複数の波形の重ね合わせで生成されているので、それぞれの波形の動きがわかるようにEpicycleで表現すると以下。

複数の円があり、その中で半径の直線が回転しています。F=20という設定なので20種類の波形があり、正方形の中心から外側に向けて円も20個あります（表示上小さすぎて見えないですが）。各円の中で独自の係数Cに対応した周期で回転しています。

C = [-39, -35, -31, -27, -23, -19, -15, -11, -7, -3, 1, 5, 9, 13, 17, 21, 25, 29, 33, 37]

ですが、絶対値順に並び替えて、

C = [1, -3, 5, -7, 9, -11, 13, -15, 17, -19, 21, -23, 25, -27, 29, -31, 33, -35, 37, -39]

にしています。

円の半径はCの要素をcとすれば、

r = abs(1/c² * exp(i * c * θ) / sum(1/C²))

　= abs(1/c² / sum(1/C²))

で求めているため、正方形の中心から外側にいくほど小さな円になります。

c=1なら反時計回りに周期1、c=-3なら時計回りに周期3となります（cが正なら反時計回り、負なら時計回り）。

上記アニメーションのコード：

p = n * 20
theta = np.linspace(0, 2*np.pi, p+1)
F = 20
C = np.array([(1 + f * n) for f in range(-F//2, F//2)])
AS = np.argsort(abs(C))
C = C[AS]
print('C =', C)

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
ax.grid()
ax.axis('equal')

ims = []
cir = np.linspace(0, 2*np.pi, 65)
LN = np.array([])
for t in theta:
    VR = np.array([])
    V = np.array([0])
    v = 0
    for c in C:
        vr = 1 / c**2 * np.exp(1j * c * t) / sum(1 / C**2)
        VR = np.append(VR, vr)
        v += vr
        V = np.append(V, v)
     
    P1 = ax.plot(V.real, V.imag, c='blue')
    
    cr = abs(VR) * np.exp(1j * C * np.vstack(cir))
    P2 = ax.plot(V[:-1].real + cr.real, V[:-1].imag + cr.imag, c='gray', alpha=0.5)
    
    LN = np.append(LN, V[-1])
    P3 = ax.plot(LN.real, LN.imag, c='tab:blue')
    P4 = ax.scatter(V[-1].real, V[-1].imag, s=40, c='red')

    ims.append(P1 + P2 + P3 + [P4])
    
ani = animation.ArtistAnimation(fig, ims, interval=100)
ani.save('epi.gif', writer="imagemagick")
#ani.save('epi.mp4', writer="ffmpeg")
plt.close()
HTML(ani.to_jshtml())

フーリエ級数の数が少ない場合：

F=2、C = [1, -3]であれば、二つの波形の合成（二つの円による軌跡）になり、フーリエ級数が少ないため正方形の形は以下のようにやや丸くなってしまいます（不完全な正方形）。

C = [1, -3]であるため、大きい円は周期1で反時計回りで回っているのに対して、小さな円は周期3で時計回りに回っているのがわかります。

それぞれの円の半径は、

r(1) = (1/1)²/ ((1/1)² + (1/(-3))²) = 0.9

r(-3) = (1/(-3))² / ((1/1)² + (1/(-3))²) = 0.1

になります。