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

188 6 SolvingNonlinearAlgebraicEquations x = linspace(0, 4, 10001) y = f(x) root = None # Initialization for i in range(len(x)-1): if y[i]*y[i+1] < 0: root = x[i] - (x[i+1] - x[i])/(y[i+1] - y[i])*y[i] break # Jump out of loop if root is None: print ’Could not find any root in [%g, %g]’ % (x[0], x[-1]) else: print ’Find (the first) root as x=%g’ % root (See thefilebrute_force_root_finder_flat.py.) Note the nice use of setting root to None: we can simply test if root is None to see ifwe founda root andoverwrote theNonevalue, or ifwedidnotfind anyrootamong the testedpoints. Running thisprogramwith somefunction,sayf.x/D e x2 cos.4x/ (whichhas a solution at x D 8 ), gives the root 0.392701, which has an error of 1:9 10 6. Increasing the number of pointswith a factor of ten gives a rootwith an error of 2:4 10 8. After such aquick“flat” implementationof an algorithm,we should always try tooffer thealgorithmasaPythonfunction,applicable toaswideaproblemdomain aspossible. Thefunctionshould takef andanassociated interval Œa;b as input,as well asanumberofpoints (n), and returna list ofall the roots in Œa;b .Here isour candidate foragood implementationof thebrute force rootingfindingalgorithm: def brute_force_root_finder(f, a, b, n): from numpy import linspace x = linspace(a, b, n) y = f(x) roots = [] for i in range(n-1): if y[i]*y[i+1] < 0: root = x[i] - (x[i+1] - x[i])/(y[i+1] - y[i])*y[i] roots.append(root) return roots (See thefilebrute_force_root_finder_function.py.) This timeweuseanotherelegant techniqueto indicate if rootswerefoundornot: roots is anempty list if the rootfindingwasunsuccessful,otherwise it containsall the roots. Applicationof the function to thepreviousexamplecanbecodedas def demo(): from numpy import exp, cos roots = brute_force_root_finder( lambda x: exp(-x**2)*cos(4*x), 0, 4, 1001) if roots: print roots else: print ’Could not find any roots’