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

6.4 TheBisectionMethod 199 asingle functioncall (f(x1)) is required ineachiterationsincef(x0)becomesthe “old”f(x1)andmaysimplybecopiedasf_x0 = f_x1 (theexception is thevery first iterationwhere two functionevaluationsareneeded). Runningsecant_method.py,gives the followingprintouton thescreen: Number of function calls: 19 A solution is: 3.000000 AswiththefunctionNewton,weplacesecantinthefilenonlinear_solvers. py foreasy importanduse later. 6.4 TheBisectionMethod NeitherNewton’smethodnor the secantmethodcanguarantee that anexisting so- lutionwill be found (seeExercises 6.1 and 6.2). The bisectionmethod, however, does that. However, if thereare several solutionspresent, itfindsonlyoneof them, just asNewton’smethod and the secantmethod. The bisectionmethod is slower than theother twomethods, so reliability comeswithacostof speed. To solvex2 9 D 0, x 2 Œ0;1000 , with the bisectionmethod,we reason as follows. Thefirstkey ideais that iff.x/Dx2 9 iscontinuousonthe intervaland the functionvalues for the interval endpoints (xL D 0,xR D 1000) haveopposite signs, f.x/ must cross thex axis at least once on the interval. That is, we know there is at least onesolution. The second key idea comes fromdividing the interval in two equal parts, one to the left and one to the right of themidpointxM D 500. By evaluating the sign off.xM/, wewill immediately knowwhether a solutionmust exist to the left or right ofxM. This is so, since iff.xM/ 0,weknowthatf.x/has to cross thex axis betweenxL andxM at least once (using the sameargumentas for theoriginal interval). Likewise, if insteadf.xM/ 0, we know thatf.x/ has to cross thex axisbetweenxM andxR at least once. In any case, wemay proceedwith half the interval only. The exception is if f.xM/ 0, in which case a solution is found. Such interval halving can be continued until a solution is found. A “solution” in this case, is when jf.xM/j is sufficiently close to zero,more precisely (as before): jf.xM/j < , where is a small numberspecifiedby theuser. The sketched strategy seems reasonable, so let uswrite a reusable function that can solveageneralalgebraicequationf.x/D0 (bisection_method.py): def bisection(f, x_L, x_R, eps, return_x_list=False): f_L = f(x_L) if f_L*f(x_R) > 0: print "Error! Function does not have opposite \ signs at interval endpoints!" sys.exit(1) x_M = float(x_L + x_R)/2.0 f_M = f(x_M) iteration_counter = 1