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 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 halley(f, fp, fpp, x, tolerance=1.0e-12): global error error = [] corr = 10 if f(x)==0.0: return f(x) while abs(corr) > tolerance: fx = f(x) fpx = fp(x) fppx = fpp(x) corr = (-2.0*fx) / ( fpx + np.sign(fpx) * np.sqrt(fpx**2.0 - 2.0*fx*fppx) ) x += corr error.append(abs(corr)) return x def secant(f, x0, x1, tolerance=1.0e-12): global error error = [] corr = 10 while abs(corr) > tolerance: fx0 = f(x0) fx1 = f(x1) fpx = (fx1 - fx0)/(x1 - x0) corr = -fx1/fpx x0 = x1 x1 += corr error.append(abs(corr)) return x1 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) r4 = halley(f1, f1p, f1pp, 3) print "halley =",r4 r4_err = np.array(error) r5 = secant(f1, 0, 3) print "secant =",r5 r5_err = np.array(error) plt.figure() plt.semilogy(r1_err, 'o-') plt.semilogy(r2_err, 'o-') plt.semilogy(r3_err, 'o-') plt.semilogy(r4_err, 'o-') plt.semilogy(r5_err, 'o-') plt.legend(('newton', 'bisection', 'regula falsi', 'halley', 'secant')) # 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) r4 = halley(f2, f2p, f2pp, 2) print "halley =",r4 r4_err = np.array(error) r5 = secant(f2, 0, 2) print "secant =",r5 r5_err = np.array(error) plt.figure() plt.semilogy(r1_err, 'o-') plt.semilogy(r2_err, 'o-') plt.semilogy(r3_err, 'o-') plt.semilogy(r4_err, 'o-') plt.semilogy(r5_err, 'o-') plt.legend(('newton', 'bisection', 'regula falsi', 'halley', 'secant')) # 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) r4 = halley(f3, f3p, f3pp, 3) print "halley =",r4 r4_err = np.array(error) r5 = secant(f3, 0.5, 3) print "secant =",r5 r5_err = np.array(error) plt.figure() plt.semilogy(r1_err, 'o-') plt.semilogy(r2_err, 'o-') plt.semilogy(r3_err, 'o-') plt.semilogy(r4_err, 'o-') plt.semilogy(r5_err, 'o-') plt.legend(('newton', 'bisection', 'regula falsi', 'halley', 'secant')) plt.show()