bitterharvest’s diary

A Bitter Harvestは小説の題名。作者は豪州のPeter Yeldham。苦闘の末に勝ちえた偏見からの解放は命との引換になったという悲しい物語

ライプニッツの公式から円周率を求める

1.円周率を求める

ライプニッツは一七世紀の哲学者、数学者、科学者であり、政治家、外交官でもあった。彼は、ドイツ東部のライプチヒの出身である。

Haskellを学んでいる人であれば、モナドという言葉を知っていることと思うが、哲学用語でのモナドという概念を打ち出した学者である。ニュートンとは独立に微積分法を発見し、二進法を研究したことでも知られている。

ここで取り上げるのは、円周率を求めるライプニッツの公式である。彼によれば、円周率は次の式により求められる。

\[ \pi = 4 - \frac{4}{3} + \frac{4}{5} - \frac{4}{7} + \frac{4}{9} - \frac{4}{11} + ... \]

Haskellで円周率を求めるために、この式を分解してみる。

分子は、4, -4, 4, -4 を繰り返している。

分母は、1, 3, 5, 7, 9と1から始まる奇数の数列である。

Haskellでは、分子は(iterate negate 4)で表すことができ、分母は[1,3..]で表すことができる。

ライプニッツの公式での項の列は、Haskellを用いると

zipWith (/) (iterate negate 4) [1,3..]

となる。ここで、zipWithは、二つの数列の対応する項に演算(/)を施して新たな数列を作る。例えば、3番目までの数列を取り出すためには、

take 3 $ zipWith (/) (iterate negate 4) [1,3..]

2000項までの数列を作り、その総和は

sum $ take 2000 $ zipWith (/) (iterate negate 4) [1,3..]

で求めることができる。実際に実行すると、3.1410926536210413となる。本当の値は、3.14159265359であるので、小数点3桁まで等しくなっていることが分かる。

1.1 問題

誤差が0.01以下になるためには、何項目まで求めたらよいか?適当な関数を作成して、答えなさい。

2.フィボナッチ数を求める

フィボナッチ数は

\begin{align*} fib_1 & = 1, \\ fib_2 & = 1, \\ fib_n & = fib_{n-1} + fib_{n-2} \end{align*}

で表すことができる。即ち、フィボナッチ数は、その前の二つの数の和となっている。フィボナッチ数を数列(フィボナッチ数列と呼ばれる)で表すと、

[1,1,2,3,5,8,13,..]

となる。ライプニッツの公式から円周率を求めるときに、zipWithを学んだので、フィボナッチ数列をzipWithで表すことを考える。

円周率の計算の時と同様に、フィボナッチ数列を二つの数列に分解してみる。以下のようになる。

          1    1    1    1    2    3    5    8

    +               1    2    3    5    8   13
――――――――――――――――――――――――――――――――――――――――――――――――――
          1    1    2    3    5    8   13   21

フィボナッチ数列をfibと表すと、加えられる方の数列はhaskellでは1:1:fibで表すことができる。また、加える方の数列はフィボナッチ数列の先頭を取り除いた数列であるのでhaskellではtail fibと表される。(Haskellでは、数列のように、データが並んでいるものをリストと呼んでいるが、リストの表し方には幾通りかある。例えば、1から5までの数字の列は、1:2:3:4:5:[]とあらわすことが基本であるが、見やすくするために、[1..5], [1,2..5],1:2:[3..5]などと表すことができる。)

加えた結果は、1:1の後に、数列fibとtail fibを項ごとに加えて得た数列zipWith (+) fib (tail fib)となる。

従って、フィボナッチ数列

fib = 1:1:zipWith (+) fib (tail fib)

となる。最初の10項を取り出したければ、

take 10 fib

とすればよい。

2.1 問題

ペル数は、\begin{align*} P_1 & = 1, \\ P_2 & = 2, \\ P_n & = 2P_{n-1} + P_{n-2} \end{align*}で表される。これを求める関数を定義しなさい。

3. 無限級数

ライプニッツの公式で無限級数の計算の仕方を学んだので、\( \sin x \) の無限級数について考える。これは次の式により求めることができる。

\[ \sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + ...\].

3.1 分母の数列を得る

上記の無限級数で、分母は、1,3,5,7,9...のそれぞれの階乗を取ったものである。これは次のように求めることができる。

map (\n -> product $ enumFromTo 1 n)  [1,3..]


上の関数で、enumFromTo 1 nはHaskellでの[1..n]と同じで、1からnまでの数からなるリスト(数列)、即ち、[1,2,3,4,..,n]を求めるものである。productは与えられたリストのすべての数の積である。従って、\( n! \)を与える。mapは与えられたリストの要素それぞれに対して、第一引数で定義された関数、この場合には、(\n -> product $ enumFromTo 1 n) 、を施す。即ち、[1,3..]の要素のそれぞれに対して、階乗の計算を行う。

例えば、最初の四つの項を取り出すために、

take  4 $ map (\n -> product $ enumFromTo 1 n)  [1,3..]


とすると、[1,6,120,5040]を得る。

3.2 分子の数列を得る

分子の方は、

(\x -> iterate (*(x*(-x))) x)

となる。ここで、次の項は、前で得た項に\( - x^2 \)をかけたもの、即ち、(*(x*(-x)))である。

ちなみに、x=2の場合は

take 4 $ (\x -> iterate (*(x*(-x)))  x) 2

とすると、

[2,-8,32,-128]

を得る。

3.3 \( \sin \)関数の級数を得る

ここでzipWithを利用し、分子と分母のそれぞれの数列の対応する項を除算(/)し、その結果を数列として返す関数sequenceは次のようになる。

sequence y = zipWith (/) ((\x -> iterate (*(x*(-x)))  x) y) (map (\n -> product $ enumFromTo 1 n) [1,3..])

従って、近似として最初の10項を求め、それらの総和を\( \sin y \)とすると、次のように表せる。

sin y = sum $ take 10 $sequence y

例えば、\( \sin 1.5706 (= \frac{\pi}{2}) \)

を求めると、

0.9999999807278945

を得る。

3.4 問題:余弦関数\( \cos x \)

\( \cos x \) は次の無限級数で定義される。

\[ \cos x = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!} + ...\].

関数\( \cos x \) を定義しなさい。

3.5 問題:ネイピア数

ネイピア数\( e \) は次の無限級数で定義される。

\[ e = \frac{1}{0!} + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + ...\].

この級数を表す関数napierを定義しなさい。