# hplgit/scipro-primer

Switch branches/tags
Nothing to show
Fetching contributors…
Cannot retrieve contributors at this time
138 lines (113 sloc) 3.86 KB
 def Newton_basic(f, dfdx, x, eps=1E-7): n = 0 # iteration counter while abs(f(x)) > eps: x = x - f(x)/dfdx(x) n += 1 return x, f(x), n def test_Newton_basic_precomputed(): from math import sin, cos, pi def f(x): return sin(x) def dfdx(x): return cos(x) x_ref = 0.000769691024206 f_x_ref = 0.000769690948209 n_ref = 3 x, f_x, n = Newton_basic(f, dfdx, x=-pi/3, eps=1E-2) tol = 1E-15 # tolerance for comparing real numbers assert abs(x_ref - x) < tol # is x correct? assert abs(f_x_ref - f_x) < tol # is f_x correct? assert n == 3 # is n correct? def test_Newton_basic_linear(): """Test that a linear func. is handled in one iteration.""" f = lambda x: a*x + b dfdx = lambda x: a a = 0.25; b = -4 x_exact = 16 eps = 1E-5 x, f_x, n = Newton_basic(f, dfdx, -100, eps) tol = 1E-15 # tolerance for comparing real numbers assert abs(x - 16) < tol, 'wrong root x=%g != 16' % x assert abs(f_x) < eps, '|f(root)|=%g > %g' % (f_x, eps) assert n == 1, 'n=%d, but linear f should have n=1' % n def Newton(f, dfdx, x, eps=1E-7, maxit=100): if not callable(f): raise TypeError( 'f is %s, should be function or class with __call__' % type(f)) if not callable(dfdx): raise TypeError( 'dfdx is %s, should be function or class with __call__' % type(dfdx)) if not isinstance(maxit, int): raise TypeError('maxit is %s, must be int' % type(maxit)) if maxit <= 0: raise ValueError('maxit=%d <= 0, must be > 0' % maxit) n = 0 # iteration counter while abs(f(x)) > eps and n < maxit: try: x = x - f(x)/float(dfdx(x)) except ZeroDivisionError: raise ZeroDivisionError( 'dfdx(%g)=%g - cannot divide by zero' % (x, dfdx(x))) n += 1 return x, f(x), n def test_Newton_numerics(): """Test that linear function is handled in one iteration.""" f = lambda x: a*x + b dfdx = lambda x: a a = 0.25; b = -4 x_exact = 16 eps = 1E-5 x, f_x, n = Newton(f, dfdx, -100, eps) tol = 1E-15 # tolerance for comparing real numbers assert abs(x - 16) < tol, 'wrong root x=%g != 16' % x assert abs(f_x) < eps, '|f(root)|=%g > %g' % (f_x, eps) assert n == 1, 'n=%d, but linear f should have n=1' % n def test_Newton_divergence(): from math import tanh f = tanh dfdx = lambda x: 10./(1 + x**2) x, f_x, n = Newton(f, dfdx, 20, eps=1E-4, maxit=12) assert n == 12 assert x > 1E+50 def test_Newton_div_by_zero1(): from math import sin, cos f = cos dfdx = lambda x: -sin(x) success = False try: x, f_x, n = Newton(f, dfdx, 0, eps=1E-4, maxit=1) except ZeroDivisionError: success = True assert success import nose.tools as nt def test_Newton_div_by_zero2(): from math import sin, cos f = cos dfdx = lambda x: -sin(x) nt.assert_raises( ZeroDivisionError, Newton, f, dfdx, 0, eps=1E-4, maxit=1) def test_Newton_f_is_not_callable(): success = False try: Newton(4.2, 'string', 1.2, eps=1E-7, maxit=100) except TypeError as e: if "f is " in e.message: success = True def test_Newton_dfdx_is_not_callable(): nt.assert_raises_regexp( TypeError, "dfdx is ", Newton, lambda x: x**2, 'string', 1.2, eps=1E-7, maxit=100) def test_Newton_maxit_is_not_int(): nt.assert_raises_regexp( TypeError, "maxit is ", Newton, lambda x: x**2, lambda x: 2*x, 1.2, eps=1E-7, maxit=1.2) def test_Newton_maxit_is_neg(): nt.assert_raises_regexp( ValueError, "maxit=-2 <= 0", Newton, lambda x: x**2, lambda x: 2*x, 1.2, eps=1E-7, maxit=-2) if __name__ == '__main__': test_Newton_basic_precomputed()