63 lines
1.0 KiB
Python
63 lines
1.0 KiB
Python
import numpy as np
|
|
|
|
def bisec(f, a, b, eps, nmax):
|
|
"""
|
|
Takes f:[a,b] -> |R and tries to compute the null point
|
|
using bisection.
|
|
"""
|
|
|
|
if(f(a)*f(b) >= 0):
|
|
raise ValueError("f(a)*f(b) >= 0")
|
|
|
|
x_minus = a
|
|
x_plus = b
|
|
|
|
if(a > b):
|
|
x_minus = b
|
|
x_plus = a
|
|
|
|
for i in range(nmax):
|
|
x_minus, x_plus, err = bisec_one(x_minus, x_plus, f, eps)
|
|
if(err < eps):
|
|
break
|
|
if(err < eps):
|
|
return x_minus, err
|
|
raise ValueError("bisection hit nmax")
|
|
|
|
|
|
def bisec_one(x_minus, x_plus, f, eps):
|
|
a = x_minus
|
|
b = x_plus
|
|
|
|
x_m = (b + a) / 2
|
|
y_m = f(x_m)
|
|
|
|
if(y_m < 0):
|
|
if(y_m > -1*eps):
|
|
return x_m, x_m, -y_m
|
|
|
|
if(f(x_m)*f(b) >= 0):
|
|
return a, x_m, -y_m
|
|
return x_m, b, -y_m
|
|
|
|
if(y_m < eps):
|
|
return x_m, x_m, y_m
|
|
|
|
if(f(a)*f(x_m) >= 0):
|
|
return x_m, b, y_m
|
|
return a, x_m, y_m
|
|
|
|
|
|
if( __name__ == "__main__"):
|
|
f1 = lambda x: x
|
|
f2 = lambda x: x**3
|
|
f3 = lambda x: -x + 1
|
|
f4 = lambda x: -x**3 + 1
|
|
f4 = lambda x: (x - 2)*np.exp(-x**2)
|
|
|
|
fs = [f1, f2, f3, f4]
|
|
for f in fs:
|
|
print(bisec(f, -12, 10, 0.0000001, 100))
|
|
|
|
print(bisec(f4, 1.2, 2.4, 0.001, 100))
|