方程式の根の計算（バイセクション法とニュートン法）
パラメータ付きの関数はオブジェクトとして表すと扱いやすい。例えば、べき関数
は、クラス型
class Func:
def __init__(self, a):
self.a = a
def v( self, x ):
return pow(x, self.a)
で表すと、２次関数は
f = Func(2.0)
と表され、ｘにおける値は
f.v(x)
で与えられる。
いま、方程式の根、すなわち、ある値
に対して
となる
の値を求める方法としてバイセクション法を用いる場合について考える。方程式の関数をオブジェクトとして表し、これにバイセクション法を適用する形でプログラムを書いてみた。用意してみる。ここで、バイセクション法とは、区間
内に根があるとき、区間の中点
によって区間を２分（バイセクション）して、
と
のいずれに根が存在するかの判定をして、存在する方の区間を新しい区間として、区間の2分を繰り返すものである（図１）。
図１
バイセクション法のプログラムを以下のように用意した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ check_bisect.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
def bisect( f, c, Lb, Ub ):
s = 1.0 if f.v(Ub) >
f.v(Lb) else -1.0
# s = 1 when f.v(x) is
an increasing functin
while ((abs(Ub - Lb) >
(1.0e-15) * (abs(Ub) + abs(Lb)))
and
(abs(Ub) + abs(Lb) > 1.0e-14)):
m = 0.5 * (Ub + Lb)
if (s * (f.v(m) - c)) > 0.0:
Ub = m
else:
Lb = m
return 0.5 * (Ub + Lb)
class Func:
def __init__(self, a):
self.a = a
def v( self, x ):
return pow(x, self.a)
f = Func(2.0)
# f(x) = x ** 2.0
bs = bisect( f, 4.0, 0.0, 5.0 )
print(bs)
bs = bisect( f, 4.0, -5.0, 0.0 )
print(bs)
f = Func(0.5)
# f(x) = x ** 0.5
bs = bisect( f, 3.0, 0.0, 10.0 )
print(bs)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
関数
def bisect( f, c, Lb, Ub ):
の第１引数に根を求める関数を表すオブジェクトを設定する。第２引数 c は、根に対する関数値である。LbとUbは、根を探索するときの最初の区間の両端の値である。LbとUbがwhileの条件
while ((abs(Ub - Lb) > (1.0e-15) * (abs(Ub) + abs(Lb)))
and
(abs(Ub) + abs(Lb) > 1.0e-14)):
を満たさなくなるまで区間の２分が繰り返される。すなわち、次の条件
(abs(Ub - Lb) > (1.0e-15) * (abs(Ub) + abs(Lb)))
を満たさなくなり精度が１５桁以上になるか、あるいは次の条件
(abs(Ub) + abs(Lb) > 1.0e-14)
を満たさなくなり、区間の両端の値の絶対値が１４桁の精度で０なるまで、根を含む区間の２分が繰り返される。
上のプログラムcheck_bisect.pyを実行すると図２の結果を得る。
図２
f = Func(2.0) # f(x) = x ** 2.0
bs = bisect( f, 4.0, 0.0, 5.0 )
print(bs)
は、
の根bsを、Lb=0.0、Ub=5.0の区間から始めて求めるものであり、バイセクション法における収束条件の精度範囲で2.0が得られていることが分かる。
bs = bisect( f, 4.0, -5.0, 0.0 )
print(bs)
は、
の根bsを、区間Lb=-5.0、Ub=0.0から始めて求めている。負の範囲での根-2.0が収束条件の精度範囲で得られている。
次の部分
f = Func(0.5)
# f(x) = x ** 0.5
bs = bisect( f, 3.0, 0.0, 10.0 )
print(bs)
は、方程式
の解を求めるもので、収束条件として設定されている有効桁数の精度で bs=9 得られている。
バイセクション法より効率の良い方法としてニュートン法がある。これは、関数
の点
における接線
が、
と交わるｘ座標
を次の根の推定値の更新値とするものである（図３）。
図３
ニュートン法のプログラム例がcheck_Newton.pyである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ check_Newton.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
def Newton( f, df, c,
Lb, Ub ):
x1 = Lb + 0.3 * (Ub - Lb)
x0 = Lb + 0.6 * (Ub - Lb)
while ((abs(x1 - x0) >
(1.0e-14) * (abs(x1) + abs(x0)))
and
((abs(x1) + abs(x0)) > 1.0e-14)):
x0 =
x1
x1 =
x0 + (c - f.v(x0)) / df.v(x0)
return x1
class Func:
def __init__(self, a):
self.a = a
def v(self, x):
return pow(x, self.a)
class DFunc:
def __init__(self, a):
self.a = a
def v(self, x):
return self.a * pow(x, self.a - 1.0)
f = Func(2.0)
df = DFunc(2.0)
v = Newton( f, df,
2.0, 0.0, 10.0 )
print (v)
print (math.sqrt(2.0))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
根を求める関数がクラス型Funcで、関数の導関数がクラス型DFuncで表されている。クラス型Funcはべき関数を表すもので、上のバイセクション法におけるものと同じである。クラス型DFuncはその導関数
を表している。ニュートン法を行う関数
Newton( f, df, c, Lb, Ub )
のパラメータは、導関数のパラメータdfが第２パラメータとして追加されていることを除けばバイセクション法の関数bisectと同じである。ニュートン法では区間の両端の値LbとUbは手続きの初期値の算出においてのみ用いられているが、後で示すニュートン法の改良版では積極的に用いられている。
プログラムcheck_Newton.pyの実行結果を図４に示す。
図４
ニュートン法は、かならず根に収束するという保証がない。更新を繰り返しても根に収束しそうにないときは、確実に根に収束するバイセクション法に切り替えることにする。次のプログラムcalcroot.pyでは、ニュートン法で２０回の更新を行っても収束しないときはバイセクション法に切り替えている。また、更新値がLbとUbで設定された範囲を超えるときは、更新値をLbあるいはUbに設定している。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ calcroot.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
def bisect( f, c, Lb,
Ub ):
s = 1.0 if f.v(Ub) >
f.v(Lb) else -1.0
while ((abs(Ub - Lb) >
(1.0e-15) * (abs(Ub) + abs(Lb)))
and
(abs(Ub) + abs(Lb) > 1.0e-14)):
m = 0.5 * (Ub + Lb)
if (s * (f.v(m) - c)) > 0.0:
Ub = m
else:
Lb = m
return 0.5 * (Ub + Lb)
def Newton( f, df, c,
Lb, Ub ):
x = Lb + 0.3 * (Ub - Lb)
x0 = Lb + 0.6 * (Ub - Lb)
niter = 0
while ((abs(x - x0) >
(1.0e-14) * (abs(x) + abs(x0)))
and
((abs(x) + abs(x0)) > 1.0e-14)):
x0 =
x
x =
x0 + (c - f.v(x0)) / df.v(x0)
if x
< Lb:
x = Lb
if x
> Ub:
x = Ub
niter += 1
if
niter > 20:
break
if niter <= 20:
return x
root = bisect( f, c, Lb, Ub
)
return root
class Func:
def __init__(self, a):
self.a = a
def v(self, x):
return pow(x, self.a)
class DFunc:
def __init__(self, a):
self.a = a
def v(self, x):
return self.a * pow(x, self.a - 1.0)
f = Func(2.0)
df = DFunc(2.0)
v = Newton( f, df,
2.0, 0.0, 10.0 )
print (v)
print (math.sqrt(2.0))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムcalcroot.pyの実行結果は、check_Newton.pyと同じである。
このウェブサイトのプログラムは圧縮ファイルRootPrgs.zipにまとめた。ファイル名RootPrgs.zipをクリックしてダウンロードできる。