Page - 198 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

198 6 SolvingNonlinearAlgebraicEquations Comparing (6.3) to the graph in Fig. 6.2, we see how two chosen starting points (x0 D 1000, x1 D 700, and corresponding function values) are used to compute x2. Oncewehavex2,wesimilarlyusex1 andx2 to computex3. AswithNewton’s method, the procedure is repeated until f.xn/ is below some chosen limit value, or some limit on the number of iterations has been reached. We use an iteration counterhere too,basedon thesame thinkingas in the implementationofNewton’s method. We can store the approximations xn in an array, but as in Newton’s method, wenotice that the computationofxnC1 onlyneeds knowledgeofxn andxn 1, not “older”approximations. Therefore,wecanmakeuseofonly threevariables: x for xnC1,x1 forxn, andx0 forxn 1. Note thatx0 andx1mustbegiven (guessed) for thealgorithmtostart. Aprogramsecant_method.pythatsolvesourexampleproblemmaybewritten as: def secant(f, x0, x1, eps): f_x0 = f(x0) f_x1 = f(x1) iteration_counter = 0 while abs(f_x1) > eps and iteration_counter < 100: try: denominator = float(f_x1 - f_x0)/(x1 - x0) x = x1 - float(f_x1)/denominator except ZeroDivisionError: print "Error! - denominator zero for x = ", x sys.exit(1) # Abort with error x0 = x1 x1 = x f_x0 = f_x1 f_x1 = f(x1) iteration_counter += 1 # Here, either a solution is found, or too many iterations if abs(f_x1) > eps: iteration_counter = -1 return x, iteration_counter def f(x): return x**2 - 9 x0 = 1000; x1 = x0 - 1 solution, no_iterations = secant(f, x0, x1, eps=1.0e-6) if no_iterations > 0: # Solution found print "Number of function calls: %d" % (2 + no_iterations) print "A solution is: %f" % (solution) else: print "Solution not found!" Thenumberof functioncalls isnowrelated tono_iterations, i.e., thenumber of iterations, as2 + no_iterations, sinceweneed twofunctioncalls beforeen- tering thewhile loop,and thenonefunctioncall per loop iteration.Note that, even thoughwe need two points on the graph to compute each updated estimate, only