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) print "a=%lf c=%lf b=%lf" % (a, m, b) 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 f2(x): return np.cos(x) def f2p(x): return -np.sin(x) 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')) 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')) # 6.2.1 def f3(x): return x**3 - 13*13 def f3p(x): return 3.0*x**2 print "13**(2/3)=", newton(f3, f3p, 13), " (via newton)" print "13**(2/3)=", 13**(2./3), " (via python)" def root(x, k): f = lambda y: y**k-x fp = lambda y: k*y**(k-1.) return newton(f, fp, x) print "13**(1/3)=", root(13,3), " (via root)" print "13**(1/3)=", 13**(1./3), " (via python)" plt.show()