【プログラミング】解けない方程式を解きたい！【数値解析：ニュートン法①】

本日も、りけいのりからお届けします。

今回は，

解けない方程式を解く！

をテーマにしていきます．まずは，解けない方程式ってどんなものかを説明し，具体的に解く手順を説明していきます．

そして，次回以降，実際にC言語でプログラムを書き，解いていきます．

f:id:ReK2Science:20201001230316p:plain

解けない方程式
ニュートン法
計算例
おわりに
- 参考文献

解けない方程式

皆さん，以下の式を解いてみましょう！

$\displaystyle{ 8x + 4 = 36 }$

中学校一年生の一次方程式の問題ですね．答えは次のようになりますね．

$\displaystyle{ x = 4 }$

これは，線形方程式(linear equation)に分類されます．これに対して，非線形方程式(nonlinear equation)もあります．この，線形と非線形の言葉については別の機会に説明しますが，簡単に紹介すると，厳密には違いますが，なんとなく線形は「直線的」というイメージでいいと思います．話を戻します．非線形方程式の例は以下の通りです．

$\displaystyle{ x^4 - 7x^3 + 4x^2 - x + 5 = 0}$

$\displaystyle{ \frac{ 4x - 5 }{ x^2 - 5x +5 } + \frac{ x + 2 }{ x^2 - 1} + 3x = 0}$

$\displaystyle{ ( x - 3 )^{\frac{5}{2}} + ( x^2 - 5 )^{\frac{4}{3} }= 0 ) }$

$\displaystyle{ e^x - sin x - 2 = 0}$

もちろん，非線形方程式のすべてが解くことができないわけではありません．解くことができる方程式もたくさんあります．しかし，これから，紹介する方法を使えば，だいたいの方程式を解くーその方程式を満たす未知数の値を導くーことができます．もちろん，場合によっては解けない場合がありますが，普通は解けます．

ニュートン法

それでは，解けない方程式を解いていく手順を紹介します．まず，今回はニュートン法の手順を説明しますが，ほかにも，反復法と呼ばれる手順もあります．

まず，未知数を $\displaystyle{ x }$ として，その方程式の項を左辺にまとめ，右辺をゼロにしたものを

$\displaystyle{ f(x) = 0}$

と表現します．この方程式を満たす， $\displaystyle{ x }$ を近似的に求めていきます．

ここで，このニュートン法の理解のコツとして，図を用いて説明していきます．この $\displaystyle{ f(x) }$ を関数として，

$\displaystyle{ y = f(x) }$

と考えると，求める解は以下の図の $\displaystyle{ y = f(x) = 0}$ を満たす点，つまり，軸との交点ということがわかります．

f:id:ReK2Science:20201003165241p:plain

では，具体的な手順を見ていきます．

初期値 $\displaystyle{ x^{(0)} }$ の点で垂線を引きます．この $\displaystyle{ x^{(k)} }$ は，微分や累乗を表しているのでなく，ｋ回目の計算で求めた解という意味で，第ｋ次近似解とか言われます．この初期値 $\displaystyle{ x^{(0)} }$ は，適当に決めます．大体の場合は，どんな値にしても最終的にはうまくいきます．
曲線 $\displaystyle{ f(x) }$ との交点をAとします．
A点での接線の方程式を立てます．
$\displaystyle{ y - f(x^{(0)} )= f'(x^{(0)}) ( x - x^{(0)}) }$

より

$\displaystyle{ y = f'(x^{(0)}) ( x - x^{(0)}) + f(x^{(0)} ) }$
その接線と水平軸との交点をBとし，座標値 $\displaystyle{ x^{(1)} }$ とおきます

上の接線の式より次のように求まります．

$\displaystyle{ 0 = f'(x^{(0)}) ( x - x^{(0)}) + f(x^{(0)} ) }$
線の方程式を立てます．
$\displaystyle{ x^{(1)} = x^{(0)} - \frac{f(x^{(0)})}{f'(x^{(0)})} }$
この操作を繰り返す．
一般に， $\displaystyle{ x^{(k+1)} }$ は次のように表せます．
$\displaystyle{ x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} }$
操作を繰り返すにつれて， $\displaystyle{ x^{(k+1)} }$ は限りなく $\displaystyle{ f(x) = 0}$ の解に近づいていきます．

以上です！こんな手順です！

手順はわかったけど，まだなんかしっくりこないなあ，，，って人が多いと思います．

それは当然で，実際にこれを利用した具体例がないとイメージがわかないと思います．

計算例

具体例を，計算してみましょう．

$\displaystyle{ x^{2} = 3 }$ をニュートン法で解け．

この問題は，別に数値解析でなくても解けますね．この解は，

$\displaystyle{ x = \pm{\sqrt{3}} }$

ですね！だいたい，1.7320508・・・（ひとなみにおごれや）ですね．

では，これを，ニュートン法を用いて解いたとき，どうなるか見ていきましょう．

まず，次のように変形します．

$\displaystyle{ f(x) = x^{2} - 3 = 0 }$

この時，以下のように計算できます．

$\displaystyle{ f'(x) = 2 x }$

これらを用いて，以下のようにニュートン法の計算を行います．

$\displaystyle{ x^{(k+1)} = x^{(k)} - \frac{ (x^{(k)})^{2} - 3 }{ 2 x^{(k)} } }$

となります．ちなみにですが，私の場合以下のように変形します．計算ミスを減らすためです．累乗の計算をしなくて済みます．

$\displaystyle{ x^{(k+1)} = \frac{ 1 }{ 2 } ( x^{(k)} + \frac{ 3 }{ x^{(k)} } ) }$

では，さっそく計算してみましょう．ここでは，初期値を $\displaystyle{ x^{(0)} = 2.0 }$ とします．

\begin{eqnarray} x^{(1)} &=& \frac{ 1 }{ 2 } ( x^{(0)} + \frac{ 3 }{ x^{(0)} } ) \\ &=& \frac{ 1 }{ 2 } ( 2 + \frac{ 3 }{ 2 } ) \\ &=& 1.75 \end{eqnarray}

次も行きます！

\begin{eqnarray} x^{(2)} &=& \frac{ 1 }{ 2 } ( x^{(1)} + \frac{ 3 }{ x^{(1)} } ) \\ &=& \frac{ 1 }{ 2 } ( 1.75 + \frac{ 3 }{ 1.75 } ) \\ &=& 1.732142857 \end{eqnarray}

\begin{eqnarray} x^{(3)} &=& \frac{ 1 }{ 2 } ( x^{(2)} + \frac{ 3 }{ x^{(2)} } ) \\ &=& \frac{ 1 }{ 2 } ( 1.732142857 + \frac{ 3 }{ 1.732142857 } ) \\ &=& 1.73205081\end{eqnarray}

以上です！だいぶ， $\displaystyle{ \sqrt{3} }$ に近い値が計算で求まってきましたね！これをずっと繰り返すと，さらに精度が高まってきます！