#============================================================================
def Bisect(Func, a, b):
#----------------------------------------------------------------------------
# Determines a real root x of function Func isolated in interval [a,b] by
# the bisection method
#----------------------------------------------------------------------------
eps = 1e-10 # relative precision of root
itmax = 100 # max. no. of iterations
x = a; fa = Func(x) # is a the root?
if (abs(fa) == 0e0): return (x)
x = b; fb = Func(x) # is b the root?
if (abs(fb) == 0e0): return (x)
if (fa*fb > 0):
print("[a,b] does not contain a root!")
return (x) # [a,b] does not contain a root
# or contains several roots
for it in range(1,itmax+1):
x = 0.5e0 * (a + b) # new approximation
fx = Func(x)
if (fa*fx > 0): a = x # choose new bounding interval
else: b = x
if (((b-a) <= eps*abs(x)) or (abs(fx) <= eps)): return (x)
print("Bisect: max. no. of iterations exceeded !"); return (x)
#============================================================================
def FalsPos(Func, a, b):
#----------------------------------------------------------------------------
# Determines a real root x of function Func isolated in interval [a,b] by
# the false position method
#----------------------------------------------------------------------------
eps = 1e-10 # relative precision of root
itmax = 100 # max. no. of iterations
x = a; fa = Func(x) # is a the root?
if (abs(fa) == 0e0): return (x)
x = b; fb = Func(x) # is b the root?
if (abs(fb) == 0e0): return (x)
if (fa*fb > 0):
print("[a,b] does not contain a root!")
return (x) # [a,b] does not contain a root
# or contains several roots
for it in range(0,itmax):
x = (a*fb - b*fa)/(fb - fa) # new approximation
fx = Func(x)
if (fa*fx > 0): # choose new bounding interval
dx = x - a; a = x; fa = fx
else:
dx = b - x; b = x; fb = fx
if ((abs(dx) <= eps*abs(x)) or (abs(fx) <= eps)): return (x)
print("FalsPos: max. no. of iterations exceeded !"); return (x)
#============================================================================
def NewtonNumDrv(Func, a, b, x):
#----------------------------------------------------------------------------
# Determines a real root x of function Func isolated in interval [a,b] by
# the Newton-Raphson method using the numerical derivative. x contains on
# input an initial approximation.
#----------------------------------------------------------------------------
eps = 1e-10 # relative precision of root
itmax = 100 # max. no. of iterations
for it in range(0,itmax):
f = Func(x)
dx = eps*abs(x) if x else eps # derivation step
df = (Func(x+dx)-f)/dx # numerical derivative
dx = -f/df if abs(df) > eps else -f # root correction
x += dx # new approximation
if ((x < a) or (x > b)):
print("[a,b] does not contain a root!")
return (x) # [a,b] does not contain a root
if (abs(dx) <= eps*abs(x)): return (x) # check convergence
print("NewtonNumDrv: max. no. of iterations exceeded !"); return (x)
#============================================================================
def Secant(Func, a, b, x):
#----------------------------------------------------------------------------
# Determines a real root x of function Func isolated in interval [a,b] by
# the secant method. x contains on entry an initial approximation.
# Error code: 0 - normal execution
# 1 - interval does not contain a root
# 2 - max. number of iterations exceeded
#----------------------------------------------------------------------------
eps = 1e-10 # relative precision of root
itmax = 1000 # max. no. of iterations
x0 = x; f0 = Func(x0)
x = x0 - f0 # first approximation
for it in range(0,itmax):
f = Func(x)
df = (f-f0)/(x-x0) # approximate derivative
x0 = x; f0 = f # store abscissa and function value
dx = -f/df if abs(df) > eps else -f # root correction
x += dx # new approximation
if ((x < a) or (x > b)):
print("[a,b] does not contain a root!")
return (x) # [a,b] does not contain a root
if (abs(dx) <= eps*abs(x)): return (x) # check convergence
print("Secant: max. no. of iterations exceeded !"); return (x)
def f(x):
return (x**3-x-3)
print 'Bisection method'
print Bisect(f, 0.0, 2.0)
print 'False position method'
print FalsPos(f, 0.0, 2.0)
print 'Newton-Raphson method using the numerical derivative'
print NewtonNumDrv(f,0.0, 2.0,1.5)
print 'Secant method'
print Secant(f,0.0, 2.0,1.5)
