setup diary

_ Rで多項式やPade近似

Rで係数cnから多項式の値を求めるためには、例えば以下のようにすれば良い。

cn<- 1/gamma(1:11)
x<- -100:100/10
f<-0
for(i in 0:10)f<-f+cn[i+1]*x^i
plot(x,f,type="l")

しかし、値を求めるためにゼロで初期化したりforループを書いたりが面倒である。そこで、xのべき乗の行列を作って、行列の積を使うとスマートに書ける。

f<-outer(x,0:10,"^")%*%cn

これを使って、Pade近似をいじってみた。Pade近似というのは、関数を２つの多項式の比として近似する方法であり、n次の多項式をm次の多項式で割る場合に[n,m]次のPade近似という。例えば、Taylor展開は[n,0]次のPade近似ということもできる。sin(x)の10次のTaylor展開の係数cnは以下で表される。

cn<-c(0,1,0,-1)/gamma(1:11)

これを[4,6]次のPade近似で表してみよう。pracmaというlibraryを使うと係数が簡単に求められる。

library(pracma)
pq<-pade(rev(cn),,4,6)
pn<-rev(pq$r1)
qn<-rev(pq$r2)
f<-(outer(x,0:4,"^")%*%pn)/(outer(x,0:6,"^")%*%qn)

pracmaでは多項式の係数を逆順にしないといけないことに注意が必要である。libraryを使うのは反則かも知れないので、library無しで、係数を求めるプログラムも書いてみた。

mat<-cbind(diag(1,10,4),-sapply(0:5,function(i)c(rep(0,i),cn)[1:10] ))
pqn<-solve(mat)%*%cn[-1]
pn<-c(cn[1],pqn[1:4,1])
qn<-c(1,pqn[5:10,1])

関数によっては、Taylor展開よりもPade近似の方が、良く合う場合があり、 Rで簡単に使えるようになったので、今後は使ってみよう。