【プログラミング】ニュートン法をC言語で実行してみよう！【数値解析：ニュートン法②】

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

今回も，

解けない方程式を解く！

をテーマにしていきます．

まず，ニュートン法の復習です！

その方程式を以下のように変形します．

未知数を $\displaystyle{ x }$ とした，方程式を以下のように変形します．

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

次に初期値 $\displaystyle{ x^{(0)} }$ を定め，以下の式で繰り返していきます．

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

$\displaystyle{ x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} }$

以上です．

ここで，このニュートン法には繰り返しの処理が必ず必要です．場合によっては50回，100回繰り返す必要が出てきます．このため，コンピュータを使っていきます！

今回は，ニュートン法での方程式の解法を，C言語でコンピュータでの計算を詳しく行っていきます．

C言語のプログラムの実行などのついては，過去記事を参照してください！

f:id:ReK2Science:20200919023135j:plain

問題設定
ニュートン法
計算例
おわりに
- 参考文献

問題設定

それでは，さっそくやっていきましょう．

ソースコード

#include<stdio.h>
#include<math.h>
#pragma warning(diaable:4996)
#define PI 3.141592

int main(void)
{

 double a = 26600;
 double e = 0.75;
 double u = 398600.5;
 double t0 = 0;
 double t[25];
 double M0 = 0;
 double M[25];
 double E[25];
 double n = 0;

 int i, j= 0;

 n = sqrt(u / pow(a, 3));
 printf("%lf\n", n);

 for (i = 0; i <= 24; i++)
 {
 t[i] = 1800 * i; M[i] = n * (t[i] - t0) + M0;
 E[i] = 0;
 for (j = 0; j < 50; j++)
 E[i] = E[i] - (E[i] - e * sin(E[i]) - M[i]) / (1 - e * cos(E[i]));
 }

 //for (i = 0; i <= 24; i++)
 // printf("t-t0 = %lf\n   M = %lf\n   E = %lf\n", t[i]/60, M[i]*(180/PI), E[i]*(180/PI));

 for (i = 0; i <= 24; i++) printf("%lf\n", M[i] * (180 / PI));
 for (i = 0; i <= 24; i++) printf("%lf\n", E[i] * (180 / PI));

 return 0;
}

\\\

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

[tex:\displaystyle{ 8x + 4 = 36 }]

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

[tex:\displaystyle{ x = 4 }]

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

[tex:\displaystyle{ x^4 - 7x^3 + 4x^2 - x + 5 = 0}]

[tex:\displaystyle{ \frac{ 4x - 5 }{ x^2 - 5x +5 } + \frac{ x + 2 }{ x^2 - 1} + 3x = 0}]

[tex:\displaystyle{ ( x - 3 )^{\frac{5}{2}} + ( x^2 - 5 )^{\frac{4}{3} }= 0 ) }]

[tex:\displaystyle{ e^x - sin x - 2 = 0}]

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

ニュートン法

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

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

[tex:\displaystyle{ f(x) = 0}]

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

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

[tex:\displaystyle{ y = f(x) }]

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

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

ソースコード

#include <stdio.h>
#pragma warning(disable:4996)

int main(void)
{
printf("hello world\n");
return 0;
}

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

より

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

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

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

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

計算例

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

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

今回は，次の具体例を，計算してみましょう．

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

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

[tex:\displaystyle{ x = \pm{\sqrt{3}} }]

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

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

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

[tex:\displaystyle{ f(x) = x^{2} - 3 = 0 }]

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

[tex:\displaystyle{ f'(x) = 2 x }]

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

[tex:\displaystyle{ x^{(k+1)} = x^{(k)} - \frac{ (x^{(k)})^{2} - 3 }{ 2 x^{(k)} } }]

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

[tex:\displaystyle{ x^{(k+1)} = \frac{ 1 }{ 2 } ( x^{(k)} + \frac{ 3 }{ x^{(k)} } ) }]

では，さっそく計算してみましょう．ここでは，初期値を[tex:\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}

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

ですが，ここで，気づいてほしい問題が2つあります．

1つ目は，もう一つの解の[tex:\displaystyle{ - \sqrt{3} }]の結果が得られなかったこと.ちなみに今回も含め，大体の場合，初期値の与え方で，ほかの解も得られます．

2つ目は，初期値を[tex:\displaystyle{ x^{(0)} = 0.0 }]とすると，分母がゼロとなり，計算ができないこと．

今回はそこまで詳しい考察はしませんが，ニュートン法を利用される際は，正解が複数個あるときの初期値の与え方や，初期値の与え方による計算の変化，について検討が必要です．