###################################################### # The bisection method # Jim Carlson # http://www.math.utah.edu/~carlson) # December 14, 1995 # root( f, a, b, n) tries to find an approximate root # of f in the interval [a,b] using n # iterations of the bisection method. You should # verify that the sign of f changes from one # endpoint to the other. The output of root # is a list consisting of # # xa = the left endpoint the last interval # xb = the right endpoint the last interval # # If you started with a legitimate value for a # and b then xa and xb will be # reasonably good approximations to the root. # If the are then # # ya = f(xa), yb = f(xb) # # should be small. # # Put the comment symbol (#) in front of the # print statement if you don't want to watch # root do its work. # # Example: # # f := x -> x^2 - 2.0; # root( f, 0, 2, 10 ); # # This computes an approximation to the square root of two. root := proc( f, a, b, n ) local i, xa, xm, xb, ya, yb, ym; xa := evalf( a ); xb := evalf( b ); ya := f(a); yb := f(b); for i from 1 to n do xm := (xa+xb)/2; ym := f(xm); print( [ xa, xb ] ); if ya*ym < 0 then xb := xm; yb := f(xb); else xa := xm; ya := f(xa); fi; od; [ xa, xb ]; end; ###################################################### # root2( f, a, b, epsilon ) is like root( f, a, b, n ) # except that one gives a maximum error range epsilon # instead of a number of iterations. Note the two # kinds of loop used: for-loop above, while-loop below. # # Example: # # > f := x -> x^2 - 2.0; # > root2( f, 0, 2, 0.01 ); # # This computes an approximation to the square root of two. root2 := proc( f, a, b, epsilon ) local xa, xm, xb, ya, yb, ym; xa := evalf( a ); xb := evalf( b ); ya := f(a); yb := f(b); while ya*yb < 0 and abs(xa-xb) > epsilon do xm := (xa+xb)/2; ym := f(xm); print( [ xa, xb ] ); if ( ya*ym < 0 ) then xb := xm; yb := f(xb); else xa := xm; ya := f(xa); fi; od; [xa, xb ]; end; ######################################################