りけいのり

かがくをやさしくおもしろく

MENU

【プログラミング】ニュートン法を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)} }]は,微分や累乗を表しているのでなく,k回目の計算で求めた解という意味で,第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 }]とすると,分母がゼロとなり,計算ができないこと.

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

おわりに

いかがでしたか?

 

もし,どうしても解けない方程式があなたの目の前に現れても,あっても今日の方法を使えば,何とか太刀打ちできます.

しかし,このニュートン法には繰り返しの処理が必ず必要です.今回の計算例では,すぐに大体の値に答えが収束しましたが,場合によっては50回,100回繰り返す必要が出てきます.このため,コンピュータを使っていきます!

次回は,ニュートン法での方程式の解法を,本格的にC言語でコンピュータでの計算を詳しく行っていきましょう!

www.rek2u.com

本日も、りけいのりがお届けしました。

参考文献

1)高橋麻奈『やさしいC 第4版』 風工舎

www.amazon.co.jp

 2)安田仁彦 『数値解析基礎』 コロナ社 

www.amazon.co.jp