1 Root Finding In mathematics, a root finding algorithm is an algorithm for finding zeros of continuous functions. At the same time these algorithms could be used to solve for f(x) = g(x) by finding the roots of f(x) g(x). These methods typically take an initial estimate of the zero or range for the zero and successively approximate better estimates or tighter bounds of the range until it can determine a root of the function. Consider the polynomial x 516.8x 4+22.4x 3+700.8x 22102.4x3456 in the range [5, 5] as seen if Fig. 1. In many cases, including for this function shown here, there are no explicit ways to solve for the zeroes of this function as opposed to finding the zeros of a quadratic function or a trigonometric function. Estimates for the zeroes of this function, on the other hand, can be determined by using root finding algorithms to find each root of the function. Two methods of finding the roots of the function include the Bracketing Method and NewtonRaphson Method. The former involves finding tighter bounds for the zero by successively halving the interval until the range where the zero is located is sufficiently small. Using the derivative of the function, the latter method involves following the curvature of the function towards a 0 crossover. Neither method can find an exact solution but, depending on the threshold allowed, a really good estimate can be found. 2 Bracketing Method The Bracketing Method involves finding tighter bounds for a root of a function. As an input, it take an initial bound for the zero - one point must be a positive value and the other must be a negative value. The value for the midpoint is found; if its negative, a new bound is made with the midpoint and the original positive value, and if its positive, a new bound is made with the midpoint and original negative value. This process is then repeated until it converges on a zero. A flowchart for this method could be seen in Fig. 2. 2.1 Matlab Code 1 function [xval,yval,iteration] = roots bisection(func, interval, tolerance) 2 %*********************************************************** 3 % * 4 % FUNCTION ROOTS BISECTION * 5 % * 6 % Purpose: Calculates a root of a function using the * 7 % Bisection Method for Roots Finding * 8 % * 9 % Function call: [xval,yval,iteration] = roots bisection * 10 % (func, interval, tolerance) * 11 % * 12 % Input: func = the function handle of the function in * 13 % question * 14 % interval = the upper and lower bounds of interest * 1 Choose an Initial Estimate on the Bounds Determine the Midpoint of the Bounds Calculate f(midpoint) Replace Upper Bound with Midpoint Replace Lower Bound with Midpoint Calculate Error Take Midpoint as Root If f(midpoint) > 0 If f(midpoint) < 0 If Error < Threshold If Error > Threshold Figure 1: Graph of x 5 16.8x 4 + 22.4x 3 + 700.8x 2 2102.4x 3456 15 % tolerance = the error threshold before function is * 16 % considered converged * 17 % * 18 % Outputs: xval = x value of the estimated zero * 19 % yval = y value of the estimated zero * 20 % iteration = how many iterations needed until convergence * 21 % * 22 % David Cinciruk * 23 % 1/16/2012 * 24 %*********************************************************** 25 26 % Checks to see if the interval is between two points with different signs 27 if sign(func(interval(1))) == sign(func(interval(2))) 28 % If they have the same sign, it sets up for an error in the output 29 % using NaN 30 error = 0; 31 mid = NaN; 32 fmid = NaN; 33 else 2 Figure 2: Flowchart of the Bracketing Method for Root Finding 34 % If they don't, forces the error to infinite to enter the while loop 35 error = Inf; 36 end 37 38 % Initializes the iteration loop 39 iteration = 0; 40 41 % Initializes the first start and end points 42 lpx = interval(1); 43 upx = interval(2); 44 45 while tolerance < error 46 % Increases the number of iterations 47 iteration = iteration + 1; 48 49 % Finds the value of the function at the lower point of the interval 50 lpy = func(lpx); 51 52 % Finds the midpoint of the interval and the value at it 53 mid = (lpx + upx)/2; 54 fmid = func(mid); 55 56 % If the midpoint has the same sign as the lower point, it sets the 57 % new lower point of the interval equal to the midpoint. Otherwise it 58 % sets the new upper point of the interval to the midpoint 59 if sign(fmid) == sign(lpy) 60 lpx = mid; 61 else 62 upx = mid; 63 end 64 65 % When there's more than one iteration done (so there's at least two 66 % checked data points) error checking is done to see if the algorithm 67 % converged yet 68 if iteration > 1 69 error = abs(prev mid)/abs(prev); 70 end 71 72 % Sets the current midpoint value to the old midpoint value 73 prev = mid; 74 end 75 3 76 % Outputs the x value and y value 77 xval = mid; 78 yval = fmid; 2.2 Solving the Example Problem This function was modified slightly to be a script with extra print commands. If we were to find a zero of the function x 5 16.8x 4 + 22.4x 3 + 700.8x 2 2102.4x 3456 in the range of [1, 2] with a threshold of 106 , our code will run through the following steps: Current Iteration is 1 Current Midpoint is -1.500000 with a value of 1.106156e+03 Error is Inf Interval is now -1.500000 -1.000000 Current Iteration is 2 Current Midpoint is -1.250000 with a value of 1.791826e+02 Error is 1.666667e-01 Interval is now -1.250000 -1.000000 Current Iteration is 3 Current Midpoint is -1.125000 with a value of -2.644561e+02 Error is 1.000000e-01 Interval is now -1.250000 -1.125000 Current Iteration is 4 Current Midpoint is -1.187500 with a value of -4.444153e+01 Error is 5.555556e-02 Interval is now -1.250000 -1.187500 Current Iteration is 5 Current Midpoint is -1.218750 with a value of 6.693006e+01 Error is 2.631579e-02 Interval is now -1.218750 -1.187500 Current Iteration is 6 Current Midpoint is -1.203125 with a value of 1.113279e+01 Error is 1.282051e-02 Interval is now -1.203125 -1.187500 Current Iteration is 7 Current Midpoint is -1.195313 with a value of -1.668241e+01 Error is 6.493506e-03 Interval is now -1.203125 -1.195313 Current Iteration is 8 Current Midpoint is -1.199219 with a value of -2.781801e+00 Error is 3.267974e-03 Interval is now -1.203125 -1.199219 Current Iteration is 9 Current Midpoint is -1.201172 with a value of 4.173748e+00 Error is 1.628664e-03 Interval is now -1.201172 -1.199219 4 Current Iteration is 10 Current Midpoint is -1.200195 with a value of 6.955375e-01 Error is 8.130081e-04 Interval is now -1.200195 -1.199219 Current Iteration is 11 Current Midpoint is -1.199707 with a value of -1.043241e+00 Error is 4.068348e-04 Interval is now -1.200195 -1.199707 Current Iteration is 12 Current Midpoint is -1.199951 with a value of -1.738789e-01 Error is 2.035002e-04 Interval is now -1.200195 -1.199951 Current Iteration is 13 Current Midpoint is -1.200073 with a value of 2.608225e-01 Error is 1.017294e-04 Interval is now -1.200073 -1.199951 Current Iteration is 14 Current Midpoint is -1.200012 with a value of 4.347007e-02 Error is 5.085953e-05 Interval is now -1.200012 -1.199951 Current Iteration is 15 Current Midpoint is -1.199982 with a value of -6.520485e-02 Error is 2.543106e-05 Interval is now -1.200012 -1.199982 Current Iteration is 16 Current Midpoint is -1.199997 with a value of -1.086750e-02 Error is 1.271585e-05 Interval is now -1.200012 -1.199997 Current Iteration is 17 Current Midpoint is -1.200005 with a value of 1.630126e-02 Error is 6.357845e-06 Interval is now -1.200005 -1.199997 Current Iteration is 18 Current Midpoint is -1.200001 with a value of 2.716875e-03 Error is 3.178902e-06 Interval is now -1.200001 -1.199997 Current Iteration is 19 Current Midpoint is -1.199999 with a value of -4.075312e-03 Error is 1.589456e-06 Interval is now -1.200001 -1.199999 Current Iteration is 20 Current Midpoint is -1.200000 with a value of -6.792187e-04 5 Error is 7.947294e-07 Interval is now -1.200001 -1.200000 3 Newton-Raphson Method The Newton-Raphson Method uses the derivative of the function and a single estimate of the zero in order to find a better estimate for a zero. This method finds the zero intercept value of the tangent line at the value of interest. This zero intercept value is then set to the value of interest and this process is repeated until convergence. To find this zero-intercept value, the following equation can be used: xi+1 = xi f(xi) f 0(xi) (1) A flowchart for this method could be seen in Fig. 3. Choose an Initial Estimate of the Zero Calculate x_new = x_old - f(x_old)/f (x_old) Calculate Error Take x_new as root If Error > Threshold Figure 3: Flowchart of the Newton-Raphson Method for Root Finding 3.1 Matlab Code 1 function [xval,yval,iteration] = roots newton raphson(func, deriv, guess, tolerance) 2 3 %*********************************************************** 4 % * 5 % FUNCTION ROOTS NEWTON RAPHSON * 6 6 % * 7 % Purpose: Calculates a root of a function using the * 8 % NewtonRaphson Method for Roots Finding * 9 % * 10 % Function call: [xval,yval,iteration] = * 11 % roots newton raphson(func, deriv, guess, tolerance) * 12 % * 13 % Input: func = the function handle of the function in * 14 % question * 15 % deriv = the function handle of the derivative of the * 16 % function in question * 17 % guess = the initial guess of the zero value * 18 % tolerance = the error threshold before function is * 19 % considered converged * 20 % * 21 % Outputs: xval = x value of the estimated zero * 22 % yval = y value of the estimated zero * 23 % iteration = how many iterations needed until convergence * 24 % * 25 % David Cinciruk * 26 % 1/16/2012 * 27 %*********************************************************** 28 29 % Sets the initial error to infinity to force it into the while loop 30 error = Inf; 31 32 % Initializes the iteration variable 33 iteration = 0; 34 35 % Sets the old estimate to the input guess 36 x old = guess; 37 38 while tolerance < error 39 % Increases the iteration count 40 iteration = iteration + 1; 41 42 % Calculates the new estimate from the old estimate, the function, and 43 % the derivative 44 x new = x old func(x old)/deriv(x old); 45 46 % Calculates the error of the new estimate 47 error = abs((x new x old)/x new); 48 49 % Rewrites the old estimate with the new estimate 50 x old = x new; 51 52 % If it has gone on for 1000 iterations, it breaks out of the while 53 % loop 54 if iteration > 1000 55 break 56 end 57 end 58 59 % Sets the output variables 60 xval = x old; 61 yval = func(x old); 3.2 Solving the Example Problem This function was modified slightly to be a script with extra print commands. If we were to find a zero of the function x 5 16.8x 4 + 22.4x 3 + 700.8x 2 2102.4x 3456 with an initial estimate of -1.5 with a threshold of 106 , our code will run through the following steps: Old Guess is -1.500000e+00. New Guess is -1.209020e+00 7 Iteration is 1 Old Guess is -1.209020e+00. New Guess is -1.200010e+00 Iteration is 2 Old Guess is -1.200010e+00. New Guess is -1.200000e+00 Iteration is 3 Old Guess is -1.200000e+00. New Guess is -1.200000e+00 Iteration is 4 4 Python Lambda Functions As can be seen in the Matlab code examples, functions were passed as arguments. These functions, typically anonyous functions are defined like func = @(x)(x^5-16.8*x^4) in Matlab. Python has similar capacities with Lambda functions. These are anonymous functions that exist within the realm of a single program. One can call a lambda function like func = lambda x: x**5 - 16.8*x**4. This lambda function can then be used as an argument in another function like one would a variable. 5 Exercises Using the information provided in this handout, implement both the Bracketing Method and the Newton-Raphson Method in Python and solve the following problems using both methods. 1. Solve for the zero of f(x) = x 5 2. Solve for all the zeros of f(x) = e 2x4 x. 3. Solve for the zero of f(x) = x 2 4. Solve for the zero of f(x) = xex 2 . For Newton-Raphson, use an initial estimate larger than sqrt2/2