非線型方程式の解法

方程式 ｆ（ｘ）＝０ を解くには関数 ｙ＝ｆ（ｘ） の逆関数 ｘ＝ｆ^−１（ｙ） において ｙ＝０ とし，
      ｘ＝ｆ^−１（0） 
を求めればいいことであるが，一般に逆関数は容易には求まらない。そこで数値的解法では以下のような試行錯誤を行う。
　　 ｙ＝ｆ（ｘ） の値が　２点　ｘ＝ａ，ｘ＝ｂ　において異符号ならばその間において少なくとも１個の解を持つはずだ（もちろんｆ（ｘ）はこの区間で連続である）。ｘ＝ａ から ｘ＝ｂ までの区間をｎ等分し
　　ｈ＝（ｂ−ａ）／ｎ　とする。
ｘ＝ａ＋ｉｈ　（ｉ＝０，１，２，・・・ｎ）　における関数値 ｆ（ｘ） を計算し ｆ（ｘ）〜０ となる ｘ が得られればそれは解と言える。 ｆ（ｘ）〜０ とはεを任意に小さな正数とし｜ ｆ（ｘ）｜＜ε を満足するときである。しかし ｘ を０．００１の精度で求めるには ｈ〜０．００１すなわち ｎ〜１０００ としなければならず計算回数は膨大なものになり非能率である。そこで ｘ＝ａ から ｘ＝ｂ までしらみつぶしではなく重点的に調べることにしよう。
　　解の近傍のある点（ｘ１，ｆ（ｘ１）） において接線を引き, ｘ軸との交点を ｘ２ とすると　　　　　
　　　　　ｘ２＝ｘ１−ｆ（ｘ１）／ｆ’（ｘ１）
であることは容易にわかる。この式より ｘ１ を与えると ｘ２ が求まりその値を ｘ１ として再度 ｘ２ を計算する・・・・という繰り返しを行っていくと漸近的に解に近づいていくと考えられる。実際，数回の繰り返しで ｆ（ｘ２）〜０．００００１　となる。よってこのときのｘ２が解である。これがＮｅｗｔｏｎ法といわれる古くから有名なアルゴリズムである。ｆ’（ｘ１）＝０ の近傍ではこの式は使えないが，別の ｘ１ でやり直せばよい.。ｘ１を初期値，ε を収束判定定数と言う。

［プログラム］  主要演算部のみ

VBでは
Do　　
　　If Abs(f(x1)) < eps Then
　　　　Exit Do　　　　　　　　　　　　　　　　　　 　　→  ｘ１　を解として書き出し　
　　Else
　　　　x2 = x1 - f(x1) / g(x1)　　　　　　　　　　　　　　ｇ（ｘ）＝ｆ’（ｘ）
　　　　fmx1 = Format(x1, "##0.00000")
　　　　fmx2 = Format(f(x1), "##0.00000")
　　　　List1.AddItem fmx1 & " " & fmx2　　　　　→  ｘ１，ｘ２を書き出し
　　　　x1 = x2　　　　　　　　　　　　　　　　　　　　　→　次のｘ１ｔとして代入　繰り返し
　　　End If
Loop

Cでは
 while (fabs(f(x1)) > eps ){
　　　x2　=　x1　-　f(x1)/g(x1);
　　　i　=　i++;
　　　printf ("%3d  %10.6f  %10.6f\n", i , x2 , f(x2));
　　　x1　=　x2;
　　}

［サンプル］　ビジュアルなVBプログラムはここからダウンロード，安価なHTMLファイルはこちらへ
　　　　Ｃソースプログラム　　繰り返し　　階乗　　関数値　　曜日干支*　　Newton法*

［入門］　次のプログラムを作れ。

　�@　ｙ＝ｘ^２の値を−１≦ｘ≦１の範囲で求める。
　�A　ｙ_ｋ＝ sin(x)　（ｋ＝１）　　の値を−１≦ｘ≦１の範囲で求める。　
　　　　　＝　ｘ　　　（ｋ＝２）　　                  地味なＣプログラムのソースはここ
　�B　閏年か平年かを判定する。
　�C　ｙｙ年mm月dd日の曜日を求める。　　通日の公式は次ページ参照。　　　
　�D 惑星の軌道要素(理科年表参照)を読み込んで，公転周期・近日点・遠日点距離を計算する。 配列 構造体
　�E n! を求める関数を作り e の近似値を求める。

［問題］　次の方程式を解け。

　�@　ｘ^２＝５
　�A　x＝ｅ sin(x)＋M     　　 Ｋｅｐｌｅｒ　の方程式 （eとMは定数）
　�B　ｘｅｘｐ（−ｘ）＋２＝０
　�C　ｘ^ｘ＝１０
　�D　Planckの公式を微分してWienの法則を導け。