matsulibの日記

Ingredients as Code

Sageでのこぎり波のフーリエ級数展開

 そろそろ本格化する就活に向けて(?)数学の基礎を復習中、と言っても動画を見て暇をつぶしてるだけだが、Youtubeにある慶應大の講義は非常に分かりやすい。こういうOCWはもっと広まってくれると良いな。今日はこれを見た。

慶應大学講義 物理情報数学C 第六回 複素フーリエ級数からフーリエ変換へ - YouTube

 さて、のこぎり波のフーリエ級数展開は分かったけど、部分積分とか面倒なのでSageで試してみた。下で示すようにSageはリスト内包表記などPythonのコードが普通に使える数式処理システムで、リソースが限定されたオンライン版は誰でもすぐに使えるからPythonに慣れた僕的にはWolfram|Alphaより好き。ipythonと同様、hoge?と入力すればhogeのドキュメントが参照できるし、TABキー補完もある。
 またPythonの関数定義とは別に f(x) = 式 の形式で定義されるシンボリック関数というのもある。f(x)のフーリエ級数による近似をfk(x)とするとこう書ける。

n, t = var('n t')  # 変数の宣言
f(t) = t/pi  # 0<t<pi
b(n) = (2/pi)*integral(f(t)*sin(n*t), t, 0, pi)  # フーリエ正弦係数
fk(t) = sum(b(n)*sin(n*t) for n in range(1, 6))  # 最初の5項までの近似
print fk(t)

 結果はこんな感じ。

2/5*sin(5*t)/pi - 1/2*sin(4*t)/pi + 2/3*sin(3*t)/pi - sin(2*t)/pi + 2*sin(t)/pi


 これだけ見てものこぎり波になってるのかよく分からないので、この項数kを増やした場合も計算し、それぞれプロットして比較したい。ところでシンボリック関数とpython関数の使い分けはよく分からないが、少なくともこういう書き方はできなかった。

fk(t, k) = sum(b(n)*sin(n*t) for n in range(1, k+1))

 python関数と違って、宣言した時点で展開できなければならないらしい。なのでこの場合はpython関数を使って書く。複数のグラフを重ね合わせるには、plot()が返すオブジェクトを足しあわせて最後にshow()メソッドで表示する。

def fk(t, k):
    return sum(b(n)*sin(n*t) for n in range(1, k+1))
    
p = plot(fk(t, 2),-10,10,color='blue',legend_label='n=2')
p += plot(fk(t, 5),-10,10,color='red',legend_label='n=5')
p += plot(fk(t, 10),-10,10,color='green',legend_label='n=10')
p += plot(fk(t, 1000),-10,10,color='purple',legend_label='n=1000')
p.show()

f:id:matsulib:20131120011513p:plain:w500
 確かにkを大きくすると元ののこぎり波に近づいているように見える。ただしn=1000の時でも不連続点で飛び出している。これはギブズ現象と言って仕方のないことらしい。