import numpy as np import matplotlib.pylab as plt def newton(f, fp, x, tolerance=1.0e-12): global error error = [] corr = 10 while abs(corr) > tolerance: fx = f(x) fpx = fp(x) corr = -fx/fpx oldx = x x += corr error.append(abs(corr)) return x def bisect(f, a, b, tolerance=1.0e-12): global error error = [] fa = f(a) fb = f(b) if fa*fb > 0: raise Exception("Probably no root in interval [%lf, %lf]" % (a, b)) corr = 10 while abs(corr) > tolerance: m = 0.5*(a+b) fm = f(m) corr = a-m if fa*fm < 0: b = m fb = fm else: a = m fa = fm error.append(abs(corr)) return 0.5*(a+b) def regula_falsi(f, a, b, tolerance=1.0e-12): global error error = [] fa = f(a) fb = f(b) if fa*fb > 0: raise Exception("Probably no root in interval [%lf, %lf]" % (a, b)) corr = 10 while abs(corr) > tolerance: m = (a*fb-b*fa)/(fb-fa) fm = f(m) if fa*fm < 0: corr = b-m b = m fb = fm else: corr = a-m a = m fa = fm error.append(abs(corr)) return m def f1(x): return x**2.0 - 1.0 def f1p(x): return 2.0*x def f1pp(x): return 2.0 def f2(x): return np.cos(x) def f2p(x): return -np.sin(x) def f2pp(x): return -np.cos(x) def f3(x): return x**3.0 - 8.0*x**2.0 + 7.0*x def f3p(x): return 3.0*x**2.0 - 16.0*x + 7.0 def f3pp(x): return 6.0*x -16.0 # f1 r1 = newton(f1, f1p, 3) print "newton =",r1 r1_err = np.array(error) r2 = bisect(f1, 0, 3) print "bisect =",r2 r2_err = np.array(error) r3 = regula_falsi(f1, 0, 3) print "regula_falsi =",r3 r3_err = np.array(error) plt.figure() plt.semilogy(r1_err, 'o-') plt.semilogy(r2_err, 'o-') plt.semilogy(r3_err, 'o-') plt.legend(('newton', 'bisection', 'regula falsi')) # f2 r1 = newton(f2, f2p, 2) print "newton =",r1 r1_err = np.array(error) r2 = bisect(f2, 0, 2) print "bisect =",r2 r2_err = np.array(error) r3 = regula_falsi(f2, 0, 2) print "regula_falsi =",r3 r3_err = np.array(error) plt.figure() plt.semilogy(r1_err, 'o-') plt.semilogy(r2_err, 'o-') plt.semilogy(r3_err, 'o-') plt.legend(('newton', 'bisection', 'regula falsi')) # f3 r1 = newton(f3, f3p, 3) print "newton =",r1 r1_err = np.array(error) r2 = bisect(f3, 0.5, 3) print "bisect =",r2 r2_err = np.array(error) r3 = regula_falsi(f3, 0.5, 3) print "regula_falsi =",r3 r3_err = np.array(error) plt.figure() plt.semilogy(r1_err, 'o-') plt.semilogy(r2_err, 'o-') plt.semilogy(r3_err, 'o-') plt.legend(('newton', 'bisection', 'regula falsi')) plt.show()