Calculus I
Lesson 22: Solving Equations Numerically with Bisection and Newton's Method
In this worksheet, we will investigate methods for finding approximate solutions to equations of the form . We will look at two important methods: repeated bisection, which is based on the Intermediate Value Theorem, and a method based on linear approximation known as Newton's method.
Repeated Bisection
The idea behind the repeated bisection method is very simple: if a continuous function is negative somewhere and positive somewhere else, it must have a zero in between. The traditional first example of the method is to approximate by solving the equation .
> restart;
> f := x -> x^2 - 2;
> plot(f, 0..4);
It seems clear from the picture that has a zero somewhere in the interval (1,2). We can quickly confirm this by evaluating at these two points:
> f(1); f(2);
The function changes sign on the interval [1,2]; since it is a polynomial, and therefore continuous, it must have a zero in the open interval (1,2). We can obviously locate this zero more accurately by testing values of at the endpoints of shorter intervals, but the trick is to do this in the most efficient manner possible. It turns out that it is not possible, in general, to do any better than halving the length of the interval at each stage, so we bisect the interval in which we have found the zero, and evaluate at the midpoint. Since there was a sign-change over the whole interval, there must be one over one of the half-intervals, and we have the zero located in a shorter interval than before. We can repeat this process as often as necesary to obtain the required accuracy.
> m1 := (1 +2)/2; f(m1);
The sign change is between 1 and 3/2, so we know that 1 < < 1.5 . We now repeat the procedure, starting withe the interval [1, 3/2].
> m2 := (1 + m1)/2; f(m2);
> m3 := (m2 + m1)/2; f(m3);
> m4 := (m3 + m1)/2; f(m4);
> evalf([m3,m4]);
> m5 := (m3 + m4)/2; f(m5);
> m6 := (m5 + m4)/2; f(m6);
> evalf([m5,m6]);
Since changes sign on the interval [m5, m6], and is continuous, we now know that = 1.4 to 2 significant figures.
(Guaranteed!) It may seem as though we have done a lot of work for a small result, but you might remember the following points:
1) We are illustrating a very general method, which can be used to solve much more complicated equations.
2) In principle, the method can be used to get any desired accuracy.
3) Most importantly, the accuracy is guaranteed: we not only get an approximation to the zero, but we get upper and lower bounds. That is, the zero is guaranteed to lie within a certain interval.
You should also observe that the method is highly repetitious, and is therefore a good candidate for programming. Here is a crude implementation, which takes a function and an initial interval, checks to see whether the function changes sign on the interval, and then iterates 10 times.
> bisect10 := proc(f,a,b) local a1,b1,m,j:
> if (evalf(f(a)*f(b)) >= 0) then
> RETURN(`No sign change detected on ` (a,b)):
> else
> a1 := a: b1 := b:
> for j from 1 to 10 do
> m := (a1 + b1)/2:
> if (f(a1)*f(m) < 0) then b1 := m:
> else a1 := m:
> fi:
> od:
> print(`There is a zero in ` (evalf(a1), evalf(b1)) ):
> end:
> bisect10(f,2,3);
> bisect10(f,1,2);
Can you explain the next result?
> bisect10(sin,-1,4);
>
Newton's Method
Despite the advantages of repeated bisection listed towards the end of the last section, it can be slow if you need a highly accurate approximation, especially if you are doing the computations by hand. If we assume that our function is not only continuous, but differentiable, we can use the idea of linear approximation: as you can see from the picture below, if we have somehow managed to find an approximation to a zero of , the intersection of the tangent line at that point with the -axis might be expected to be a better approximation. (In the picture, the original approximation is .)
> x0 := 1;
> t1 := x -> f(x0) + D(f)(x0)*(x-x0);
> plot({f,t1}, 0..3);
Now, of course, we repeat: we find where the tangent line crosses the -axis. This is our new approximation to the zero, so we find the tangent line at this new point and see where it crosses the axis. Hopefully, by repeating the procedure as often as necessary, we can obtain any accuracy we desire.
> x1 := solve(t1(x) = 0, x);
> t2 := x -> f(x1) + D(f)(x1)*(x-x1);
> plot({f,t2}, 0..3);
Note that even after one iteration, we will need to zoom in to see the difference between the true zero of and the intersection of the tangent line with the -axis.
> x2 := solve(t2(x) = 0, x); evalf(%);
Repeat if necessary:
> t3 := x -> f(x2) + D(f)(x2)*(x - x2);
> x3 := solve(t3(x) = 0, x); evalf(%);
As a comparison, here is Maple's evaluation of (to 10 digits):
> evalf(sqrt(2));
It is clear that Newton's method is extermely fast---at least in this example. Can you think of any disadvantages it may have?