もう数日で大人の科学27号8ビットマイコンが発売ですねー。 発売後は しばらくこちらで遊ぶ予定なので、その前に区切りのところまでメモ。
こちらの記事のpsolve.pyモジュールを使って、三角関数のSin(サイン)を近似してみます。
出来上がりはコチラ
sinカーブの上に、近似した点が ちゃんと重なってるのが見えます。
今回は三角関数のsin(x) を、xの多項式によって近似したいのですけれども
sin(x) = a0*x**0 + a1*x**1 + a2*x**2+ a3*x**3 + ...
私の苦手なsinのテーラー展開をSymPyで調べると
>>> from sympy import *
>>> x=symbols("x")
>>> sin(x).series(x)
''x - x**3/6 + x**5/120 + O(x**6)''
>>> N(_)
x - 0.166666666666667*x**3 + 0.00833333333333333*x**5 + O(x**6)
x**yは xのｙ乗を表す
教科書だと階乗を使って
x-x**3/3! + x**5/5! - x**7/7! +...
てな感じに書かれているかと思います。
関連
でも、前回作成したpsolve.pyは 一次連立方程式の近似解を求めるモジュールなので、そのままではn次多項式の解を求められません。。。。
なので
x0=x**0 x1=x**1 x2=x**2 。。。
として、xのn乗をあらかじめ計算して渡してやります。
sin(x) = a0*x**0 + a1*x**1 + a2*x**2+ a3*x**3 + ... これを sin(x) = a0*x0 + a1*x1 + a2*x2 +a3*x3 + ... という一次方程式だと考える。 -3.14＜x＜3.14 の範囲でxを変化させ、適宜 一次連立方程式をいくつか作成して解く
SVMのカーネルトリックをヒントにしました。
では、上記のsinのテーラー展開のような結果が出ることを期待しつつ、Pythonのプログラムにしてみます。
from __future__ import division
from psolve import psolve
from math import sin,pi,cos
import Image,ImageDraw
w,h=640,512
im=Image.new("RGB",(w,h))
draw=ImageDraw.ImageDraw(im)
oy=h//2
lastx,lasty=0,oy+200
slist=[]
M=32
alist=[]
blist=[]
M2=M//2
draw.rectangle((0,0,w,h),fill=(255,255,255))
for t in xrange(M):
rad=pi*(t-M2)/M
#x=rad=pi*t/M/2 #正の領域だけで近似する場合はこちらの式
sn=sin(rad)
px=w*t/M
py=oy-sn*200
draw.line((lastx,lasty,px,py),fill=(255,0,0))
lastx,lasty=px,py
x=rad
#alist.append([1,x,x**2,x**3,x**4,x**5,x**6,x**7])
alist.append([1,x,x**2,x**3])
blist.append(sn)
x=1
param=psolve(alist,blist,z=[1,x,x**2,x**3],loop=100)
print param
lastx,lasty=3,oy+200
for t in xrange(M):
x=rad=pi*(t-M2)/M
#x=rad=pi*t/M/2 #正の領域だけで近似する場合はこちらの式
sn=sum([param[i]*(x**i) for i in xrange(4)])
px=w*t/M
py=oy-sn*200
draw.line((lastx+3,lasty,px+3,py),fill=(0,0,255))
draw.ellipse((px-4,py-4,px+4,py+4),fill=(0,255,0))
lastx,lasty=px,py
im.save("sin.png")
実行してみると
sin (x) ≒ 1 * -1.9938157693729221e-005+ #ほぼゼロ x * 0.99157263793573691 + (x**2)*8.8301243756742161e-005+ #これも ほぼゼロ (x**3)*-0.14566809542315673
SymPyで求めたのが
x - 0.166666666666667*x**3 +O
でしたから、かなり近い、、、*1
このパラメータで -3.14 ＜x ＜ 3.14 の範囲をグラフにして、sin関数と重ねたのが、記事冒頭の画像。 これぐらいの精度であれば、パっと見、かなり綺麗に近似できてるみたいです。
今回は xの3乗 (3次式)までで 近似しました。あまり次数を増やしても、かえって うまく近似できないようです。謎。*2
*1
と言っていいのか微妙ですけれども
*2
フーリエ級数が 波形のつなぎ目で暴れるのと同じような現象かも？と思ったけど、違うかもしれません