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()
確かにkを大きくすると元ののこぎり波に近づいているように見える。ただしn=1000の時でも不連続点で飛び出している。これはギブズ現象と言って仕方のないことらしい。