次のルートを決定するために二分法を実行する Python プログラムを作成したいと考えています。
f(x) = -26 + 85x - 91x2 +44x3 -8x4 + x5
二分法は、多項式 f(x) の根を推定するための数値的方法です。
答えを教えてくれる利用可能な疑似コード、アルゴリズム、またはライブラリはありますか?
次のルートを決定するために二分法を実行する Python プログラムを作成したいと考えています。
f(x) = -26 + 85x - 91x2 +44x3 -8x4 + x5
二分法は、多項式 f(x) の根を推定するための数値的方法です。
答えを教えてくれる利用可能な疑似コード、アルゴリズム、またはライブラリはありますか?
基本的なテクニックを示すコードは次のとおりです。
>>> def samesign(a, b):
return a * b > 0
>>> def bisect(func, low, high):
'Find root of continuous function where f(low) and f(high) have opposite signs'
assert not samesign(func(low), func(high))
for i in range(54):
midpoint = (low + high) / 2.0
if samesign(func(low), func(midpoint)):
low = midpoint
else:
high = midpoint
return midpoint
>>> def f(x):
return -26 + 85*x - 91*x**2 +44*x**3 -8*x**4 + x**5
>>> x = bisect(f, 0, 1)
>>> print(x, f(x))
0.557025516287 3.74700270811e-16
所定の許容値が達成されたときに早期に終了するには、ループの最後にテストを追加します。
def bisect(func, low, high, tolerance=None):
assert not samesign(func(low), func(high))
for i in range(54):
midpoint = (low + high) / 2.0
if samesign(func(low), func(midpoint)):
low = midpoint
else:
high = midpoint
if tolerance is not None and abs(high - low) < tolerance:
break
return midpoint
scipy.optimize.bisectを使用する以前のStackOverflowの質問で解決策を見ることができます。または、学習を目的としている場合は、前の質問のコメント投稿者が示唆しているように、二分法に関するWikipediaエントリの擬似コードはPythonで独自の実装を行うための優れたガイドです。
公差あり:
# there is only one root
def fn(x):
return x**3 + 5*x - 9
# define bisection method
def bisection( eq, segment, app = 0.3 ):
a, b = segment['a'], segment['b']
Fa, Fb = eq(a), eq(b)
if Fa * Fb > 0:
raise Exception('No change of sign - bisection not possible')
while( b - a > app ):
x = ( a + b ) / 2.0
f = eq(x)
if f * Fa > 0: a = x
else: b = x
return x
#test it
print bisection(fn,{'a':0,'b':5}, 0.00003) # => 1.32974624634
ライブ: http://repl.it/k6q
私はいくつかの変更を加えようとしています:
これは私のコードです:
import numpy as np
def bisection3(f,x0,x1,tol,max_iter):
c = (x0+x1)/2.0
x0 = c
Xs = 0.3604217029603
err_Abs = np.linalg.norm(x0-Xs)
itr = 0
f_x0 = f(x0)
f_x1 = f(x1)
xp = [] # sucesion que converge a la raiz
errores_Abs = []
errores_Aprox = []
fs = [] # esta sucecion debe converger a cero
while(((x1-x0)/2.0 > tol) and (itr< max_iter)):
if f(c) == 0:
return c
elif f(x0)*f(c) < 0:
x1 = c
else :
x0 = c
c = (x0+x1)/2.0
itr +=1
err_Abs = np.linalg.norm(x0-Xs)
err_Aprox = np.linalg.norm(x1 - x0)
fs.append(f(c))
xp.append(c)
errores_Abs.append(err_Abs)
errores_Aprox.append(err_Aprox)
return x0,errores_Abs, errores_Aprox,fs,xp
そして、私は実行の例を持っています:
f = lambda x : 3*x + np.sin(x) - np.exp(x)
X0_r1 , err_Abs_r1,err_Aprox_r1, fs_r1 , xp_r1 = bisection3(f,0,0.5,1e-5,100)
上記の Bisection-method は、許容範囲をストッパーとして変更することができます。
def samesign(a, b):
return a*b > 0
def bisect(func, low, high):
assert not samesign(func(low), func(high))
n = 20
e = 0.01 # the epsilon or error we justify for tolerance
for i in range(n):
if abs(func(low)-func(high)) > e:
midpoint = (low + high) / 2
print(i, midpoint, f(midpoint))
if samesign(func(low), func(midpoint)):
low = midpoint
else:
high = midpoint
else:
return round(midpoint, 2)
return round(midpoint, 2)