2020年4月24日

Rでマクローリン展開する ((1:n), cumprod, plot)

今回は、関数の値を近似するマクローリン展開を紹介します。マクローリン展開をまったく知らない人も、Rを使えばイメージできるようになります。やってみましょう。

目次
  1. 指数関数のマクローリン展開
  2. (1:n) 、cumprod
  3. マクローリン級数による値の近似

*  *  *

1 指数関数のマクローリン展開

指数関数とは、$y=2^x$ のように、$x$ が(べき)指数にある関数です。ネイピア数 $e$ を底とする指数関数は

$$y=e^x$$

と書けます。この関数をマクローリン展開すると下式のようになります。(導出の詳細は次回の補論で説明する予定です。)

$$e^x=1+\frac{1}{1!}x^1+\frac{1}{2!}x^2+\frac{1}{3!}x^3+\frac{1}{4!}x^4+\frac{1}{5!}x^5+...$$

式中の!は階乗です。たとえば、$3!=3\times 2\times 1=6$ です。階乗の計算をして式をきれいにすると

$$e^x=1+x+\frac{1}{2}x^2+\frac{1}{6}x^3+\frac{1}{24}x^4+\frac{1}{120}x^5+...$$

はじめてみる人は「なにこれ?」という感じだと思います。この多項式は関数の値を近似します。$x=1$ なら、左辺は $e^1=e$ となります。$e$ は2.71828くらいの値をとるネイピア数です。右辺の多項式がこの値とおなじになるというのが、マクローリン級数の意味するところです。

*  *  *

2 (1:n) 、cumprod

コードを効率よく書き、短時間で実行するのに便利なのがベクトル型の命令文です。今回新たに登場するのは (1:n)  cumprod です。

(1:n) は、1から n まで整数を並べる命令文です(njk などほかの文字でも構いません)。たとえば、n=5 のとき (1:n) は

(1:5)=1, 2, 3, 4, 5

となります。cumprod  は並んだ数を順に掛ける関数です。たとえば

 cumprod((1:5))

であれば、{1,2,3,4,5} という数の並びを順に掛けます。すなわち

1, 1*2, 1*2*3, 1*2*3*4, 1*2*3*4*5

それぞれ計算すると

cumprod((1:5))=1, 2, 6, 24, 120

となります。興味深いことに、cumprod((1:5)) は、1から5の階乗を並べたものになります。すなわち

cumprod((1:5))=1!, 2!, 3!, 4!, 5!

ベクトル型の命令文は、コードをすっきりさせます。

*  *  *

3 マクローリン級数による値の近似

ネイピア数はマクローリン級数で近似できるか、Rで検証します。下の画像をクリックするとコードが大きく表示されますので、それをみて入力してください。入力が終わりましたら「➡️Run」を押してみましょう。








うまくいけば、左下と右下に次のように表示されます。










左に示されている [1] 2.731266058e-08 はネイピア数とマクローリン級数の値の差です。指数で $2.73...\times 10^{-8}$ と表されます。これは、差がほとんどないことを意味します。下のグラフは、右のグラフを拡大したものです。マクローリン級数の値を示す○は、級数の項が増えるにしたがい、ネイピア数の値を表す点線に近づきます。


*  *  *

以下、準備、処理、出力の順にコードを説明します。まず準備からです。2行目と3行目はいつものおまじないです。Rのコードはこの2行からはじめましょう。

5行目の x=1 x の値を1に指定しています。このコードでは、x の値が1でない場合についても検証できます。6行目の n=10 は、マクローリン級数の項数の上限を指定しています。n=10 のとき x を含む10の項と第1項の1、あわせて11項を上限として級数を計算します。







*  *  *

つづいて処理です。8行目の vx は多項式の分子をまとめて計算するためのものです。n が5なら、(1:n)=1,2,3,4,5 ですので

vx=x^1, x^2, x^3, x^4, x^5

となります。ここでは n=10 ですので、x^10 まで数がならびます。9行目の vf は多項式の分母をまとめて計算するためのものです。上で説明したように、n が5なら

cumprod((1:5))=1!, 2!, 3!, 4!, 5!

となります。ここでは n=10 ですので、10! まで数がならびます。10行目の mac は、マクローリン級数の第2項から第11項を表しています。vxvf で要素ごとに割り算しますので

mac=x^1/1!, x^2/2!, x^3/3!, x^4/4!, x^5/5!, ...

となります。11行目は、mac の先頭に1をつけてできた数列の要素を順に足し上げて mac を再定義しています。言葉で説明するより数列を示した方が理解しやすいかもしれません。次のような数列が11項までつづきます。

mac
=1, 1+x^1/1!, 1+x^1/1!+x^2/2!, 1+x^1/1!+x^2/2!+x^3/3!, ...

わずか4行ですが、たくさんの工程がたたみ込まれています。




*  *  *

さいごに出力です。前々回、curve というグラフを描く関数を紹介しました。今回は plot を使います。13行目にある plot の第1要素 (0:n) は横軸の数字です。0から10までをとります。第2要素 mac は、マクローリン級数の数列です。14行目の abline は、級数の計算結果のグラフに、ネイピア数の値(h=exp(x))を点線で(lty=3)描き加える命令文です(ここでは x=1 です)。

15行目は、ネイピア数と第11項までのマクローリン級数の値の差を計算して印字する命令文です。先ほどみましたように、差はほぼ0です。




*  *  *

わずか15行のコードですが、説明すると盛りだくさんですね… Rはそれだけ効率のよい言語だともいえます。効率のよいコードがうまく動くと爽快です。是非トライしてみてください。