置き去りになったメモ

数値計算、代数、関数型プログラミングとか?

Haskellで求根アルゴリズム

数値計算 Haskell

こんにちは、やきつか(@sy_t_)です。
最近はめっきり寒くなってきたのに、いつも半袖でいる知り合いがいて怖いです。
友人がHaskellで数値計算の記事を書いていて、全く同じ内容を書くのも面白いかな、といのうがこの記事の動機です。(友人の記事→Haskellで数値計算やってみた - Qiita)

目次

はじめに

タイトルにある通り、今回は求根アルゴリズムの解説とHaskellによる実装を紹介したいと思います。
求根アルゴリズムとは、

f:\mathbb{R}\rightarrow\mathbb{R}
という関数 fにおいて、
f(x)=0
を満たす x\in \mathbb{R}を求める問題で、数値計算では二分法やニュートン法がよく知られています。

二分法

アルゴリズム

二分法は、その名の通り区間を二分していくことで解を求める手法です。
具体的には、

f(x_1)f(x_2)<0
を満たす、つまり関数の値が異符号となる x_1, x_2について区間
[x_1,x_2]
を初期区間として、中点
x_M=\frac{x_1+x_2}{2}
を考え、 x_Mが x_1と同符号なら区間を [x_M,x_2]として、 x_2と同符号なら区間を [x_1,x_M]として更新し、さらに更新された区間で中点を求めて...を繰り返して
|f(x_M)|<\epsilon
が満たされた時、 x_Mを解として出力します。ここで、 \epsilonは0に近い正の値で、この値と比べて絶対値が小さい時、 f(x)=0が満たされたとします。(この方法では厳密に0になることはないので)

実装

それでは、アルゴリズムに則って実装していきましょう。

ε :: Double
ε = 1E-15 --0に限りなく近い正の定数

bisection :: (Double -> Double) --関数
               -> Double --区間下限
               -> Double --区間上限
               -> Double --解
bisection f x1 x2
    | abs (f x') < ε = x' -- -ε<f(x)<εならば中点を解として返す
    | f x1 * f x' > 0 = bisection f x' x2 --同符号ならばx1を中点と交換
    | otherwise = bisection f x1 x' --それ以外ならばx2を中点と交換する
    where
        x' = (x1+x2)/2 --中点

実装終わり。解説していきます。
\epsilonは、0とみなす値域を決定する定数で、0に近ければ近いほど解の精度が良くなります。
実際にその処理を実装しているのが次の行です。

abs (f x') < ε = x' -- -ε<f(x)<εならば中点を解として返す

さらにその次の行は、 x_1と同符号ならば x_Mに置き換え自身に引数として渡し、 x_2の場合も同様の処理をしています。(otherwiseとなっているのは、 x_1と異符号ならば x_2と同符号なので)

    | f x1 * f x' > 0 = bisectionMethod f x' x2 --同符号ならばx1を中点と交換
    | otherwise = bisectionMethod f x1 x' --それ以外ならばx2を中点と交換する

ニュートン法

アルゴリズム

ニュートン法は、アイザック・ニュートンに由来する手法で、微分を用いて解を求める手法です。
この手法は次の漸化式によって、端的に表されます。

x_{n+1}=x_n - \frac{f(x_n)}{f'(x_n)}

f:id:yakituka:20191109120801p:plain
図で書くとこんな感じ。なぜ、上記のような漸化式になるかというと、
x_nでの接線の方程式が

y=f'(x_n)(x-x_n)+f(x_n)
なので、 y=0となる x=x_{n+1}を求めると、
0=f'(x_n)(x_{n+1}-x_n)+f(x_n)
\Leftrightarrow x_{n+1}=x_n - \frac{f(x_n)}{f'(x_n)}, \ where \ f'(x_n)\neq 0
と変形できるからです。
つまり、 x_nでの接線を引いて、 x軸と接線が交わった点を次の点とします。
これを繰り返していくことにより、解に近づいていきます。
解の判定は二分法と同様に、εを用いて行います

実装

それでは、実装していきます。

ε :: Double
ε = 1E-15 --0に限りなく近い正の定数

newton :: (Double -> Double) --関数
             -> (Double -> Double) --導関数
             -> Double --初期のx座標
             -> Double --解
newton f f' x
    | abs (f x) < ε = x -- -ε<f(x)<εならば中点を返す
    | otherwise = newton f f' x' --次の点に進めて渡す
    where
        x' = x - f x / f' x --次のx座標

実装終わり。二分法との違いは、引数に導関数が増え、1つのx座標のみを引数にしているところですね。
ただ、この方法では導関数が明確な場合しか使えないので、適当に差分で置き換えて使うことが多いと思います。
ということで、微分を差分で置き換えて実装してみるとこんな感じ

ε :: Double
ε = 1E-15 --0に限りなく近い正の定数

newton' :: (Double -> Double) --関数
             -> Double --初期のx座標
             -> Double --解
newton' f x
    | abs (f x) < ε = x -- -ε<f(x)<εならば中点を返す
    | otherwise = newton' f x' --次の点に進めて渡す
    where
        dx = 0.00001 --格子幅
        f' x = ( f x - f (x-dx) ) / dx --微分ぽいもの
        x' = x - f x / f' x --次のx座標

実行結果

f:id:yakituka:20191109131538p:plain
上から、方程式、真値、二分法、ニュートン法、ニュートン法(差分)の順です。
εは 10^{-8}にして、解に差が表れるようにしています。
ニュートン法の精度が良いですね。

感想

今回はゆるい数値計算の記事を書いてみました、懐かしい気持ちになれました。
ニュートン法(差分)はより精度の良いスキームを使って改良できそうですね、そういうのも一興だと思います。
もしここまで読んでくれた方がいたら、本当にありがとう!!