from numpy import * from matplotlib.pyplot import * def newton(f, fp, x, tolerance=1.0e-12): global error error = [] corr = 10 while abs(corr) > tolerance: corr = -f(x)/fp(x) 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: corr = -f(x1)*(x1-x0)/(f(x1)-f(x0)) x0 = x1 x1 += corr error.append(abs(corr)) return x1 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 cos(x) def f2p(x): return -sin(x) r1 = newton(f1, f1p, 3) print "newton =",r1 r1_err = array(error) r2 = bisect(f1, 0, 3) print "bisect =",r2 r2_err = array(error) r3 = secant(f1, 0, 3) print "secant =",r3 r3_err = array(error) figure('x**2-1') semilogy(r1_err, 'o-', label='newton') semilogy(r2_err, 'o-', label='bisection') semilogy(r3_err, 'o-', label='secant') legend() r1 = newton(f2, f2p, 2) print "newton =",r1 r1_err = array(error) r2 = bisect(f2, 0, 2) print "bisect =",r2 r2_err = array(error) r3 = secant(f2, 0, 2) print "secant =",r3 r3_err = array(error) figure('cos(x)') semilogy(r1_err, 'o-', label='newton') semilogy(r2_err, 'o-', label='bisection') semilogy(r3_err, 'o-', label='secant') legend() # 7.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)" show()