y=x のフーリエ級数展開を可視化する【Sage】【Python】

はじめに

小出昭一郎『量子力学（Ⅰ）』（p.60）より周期 2l の関数 f(x) は以下のように展開できます。

複素フーリエ級数

$f(x)=\sqrt{\frac{\pi}{2l^2}}\Sigma_{n}\exp(in\pi x/l)F_n$

複素フーリエ係数

$F_n=\frac{1}{2\pi}\int_{-l}^{l}\exp(-in\pi x/l)f(x)dx$

$l=1$ とし、 $f(x)=x$ の場合についてシミュレーションを行います。計算や可視化はSageShellを使います。

プログラム

$y=x$ について、自明ですが、一応グラフを書いてみます。

x=var('x')
sage: plot(x)
Launched png viewer for Graphics object consisting of 1 graphics primitive

f:id:riverta1992:20201228115738p:plain

関数を級数で表すためには係数Fnの計算が大切ですが、無限個の係数を計算するわけにもいきません。項数をNとし、与えられた関数をN個の関数の和として返す関数fourierを下記のように定義します。

def fourier(f,N,l):
    F = {}
    ns = [0]
    for i in range(1,N+1):
        ns.append(i)
        ns.append(-i)

    for n in ns:
        F[n] = ( numerical_integral(cos(-n*pi*x/l)*f(x),-l,l)[0] + I*numerical_integral(sin(-n*pi*x/l)*f(x),-l,l)[0] ) / sqrt(2*pi)

    return sqrt(pi/2/l^2) * sum([exp(I*n*pi*x/l)*F[n] for n in ns])

N=5 の場合についてグラフをプロットしてみます。

sage: plot(fourier(x,5,1).real(),xmin=-1.1,xmax=1.1,ymin=-1.1,ymax=1.1)

f:id:riverta1992:20201228125011p:plain

Nを増やしていくとどうなるか見てみます。

sage: l = [plot(x,color="gray",xmin=-1.1,xmax=1.1,ymin=-1.1,ymax=1.1)+plot(fourier(x,i,1).real(),xmin=-1.1,xmax=1.1,ymin=-1.1,ymax=1.1)+text("N={}".format(i),(-0.5,0.5),bounding_box={'boxstyle':'round','fc':'w'},fontsize=20) for i in range(1,31)]
sage: animate(l)
Launched gif viewer for Animation with 30 frames