Mathematics and Programming

Jason I. Brown

Having a computer is one thing; knowing when to use it is another. In this section we explore how Maple and mathematics might be brought together so that the "sum is greater than the parts". There are two principles to keep in mind:

If it's easier to do something by hand, do so! Use Maple in an essential way to help with the computations.

Let the mathematics drive Maple, rather than the other way around. Decide what you want to do mathematically first, and then figure out how to ask Maple for help.

Example 1: Finding Roots

Again, we'll stick to the principle of learning Maple as we go along. Consider the following typical problem.

Find the roots of .

That is, we want to find all values that make . I guess that one way to do this is to first ask Maple to solve this equation. First we'll name the function f(x).

> f:=x->x^7+2*x-5;

All this does is make the polynomial into a function .

Now let's ask Maple to solve .

> solve(f(x)=0,x);

No luck! Let's try thinking mathematically about the function . First, we see that if , then and are both negative, so that . Thus won't have any roots to the left of 0.

A Visual Approach

Let's plot say from to .

> plot(f(x),x=0..5);

Well, that what you'd expect. The dominant term when is large in the polynomial is the leading term , which is positive, so that when is large enough, the function should be positive. In fact, when is at least 2, is at least and is at least , so is at least , and thus positive. Therefore, there can't be any roots to the right of ; all the roots have to be between and . Let's try a plot on this domain.

> plot(f(x),x=0..2);

It certainly looks like there is exactly one root, somewhere between and . You can actually click on where you think the point is, and Maple will display the coordinates, if the plot is in a seperate window. For example, when I clicked, I got an value of about (it was interesting that the value was close to but not quite ). You can continue to find more precise estimates for the root by "zooming" in on the graph. That is, you can keep plotting on smaller and smaller domains that contain the root. Let's zoom a few times.

> plot(f(x),x=1..2);

> plot(f(x),x=1..1.3);

> plot(f(x),x=1.1..1.2);

By clicking (or just by looking) we see that the root is about .

An Analytic Approach

Now can you argue why mathematically there has to be a root around ? Well, a simple calculation shows that and . Now if you reach back into the dark recesses of your mind to your first year calculus class, you'll recall that you learned:

The Intermediate Value Theorem (IVT) : A continuous function on a closed interval takes on all values between and .

Any polynomial is continous (everywhere!), so if we apply the IVT to our function on , as lies between and , there must be a value of between and for which . In fact, I guess we could use the plot above to guess that the answer is about or so. To get a more accurate answer, let's think more mathematically.

The IVT tells us that as is negative and is positive, there is a root of between and . Let's see what happens in the middle:

> f(1.5);

Now we see that is negative at and positive at , so again by the IVT there is a root of between them. A good guess might be half way in between. Let's see what happens at (the half-way point):

> f(1.125);

Aha! It's negative. As is negative at and positive at , the root has to lie between them. Let's continue for a few more via Maple (we'll name the value we're looking at x):

> x:=(1.125+1.5)/2;

> f(x);

> x:=(1.125+1.3125)/2;

> f(x);

> x:=(1.125+1.21875)/2;

> f(x);

Now $f(x)$ is negative at and positive at , so a root of lies between the two values. Thus we can approximate it to 1 decimal place as . As a check , let's try :

> f(1.1);

That's pretty close to , so we're not too far off.

What about if we want to improve the acurracy, say we want an answer that is accurate to within ? It would be a lot of work at the terminal to continue, but it certainly could be done. Maple can automate repetitive tasks, so let's see how Maple can help.

Now we know that an answer is accurate to .00001 if we can determine that the root lies in an interval where . What we want to do is the following: Suppose that we are at a stage where we know that the root is between and . We first check to see if ; if it is, we print say (actually, any point in will do). If not, then since and we then take the midpoint and look at . If , then we know that the root is between and , so we replace by . If on the other hand , then we know that the root is between and , so we replace by . In any event, we repeat until we stop.

I'll write down the Maple lines we need to do this, and then talk about them. It'll be the first time we write a set of lines rather than just a single one.

> a:=1:
b:=2:
for i from 1 to 100 do
if (b-a <= 0.000001) then
print(`The root is approximately`,evalf(a));
break;
else
c:=(b+a)/2:
if f(c) < 0 then
a:=c:
else
b:=c:
fi:
fi:
od;

Now the first two lines initialize and to be the endpoints and respectively. The next line starts a repetetive loop that runs at most the number of times we are going to keep taking midpoints. The following line checks whether the interval is smaller than ; if so, we'll print out the left end point as our answer (and "break" out of the procedure and stop). Otherwise, we set to be the midpoint between and , and look at . If is negative, we replace by (as the IVT tells us that the root lies between and ), and if is positive, we replace by . In any event, we go back and repeat the process. The lines of commands make look complicates, but you'll soon catch on.

The procedure above gives an answer of , and we know that the roots is between and .

[Aside: One thing that may come to mind is whether there can be more than one root of the funtion f between 1 and 2. It seems pretty clear from all the plots that there can't be. Can we justify this mathematically? Well, if you recall Rolle's Theorem, again from that dreaded (but useful!) calculus course you took, between any two roots of a differentiable function there must be a root of its derivative. If we look at the derivative of f,

> D(f);

we see that the derivative of f is always positi ve. Thus by Rolle's Theorem, f can't have two roots (anywhere!).]

You may have notices a couple of problems with the above routine:

1. What would have happened if we landed up exactly on the root? Well, it would have been caught by the "else" statement, and we would have set (the right endpoint) equal to it, and the procedure would have continued. However, at every subsequent midpoint we choose, we see that the function would be negative at it (as the graph showed that the function is negative to the left of the root). Thus our final answer (which would be when the loop runs out) would be the left endpoint of an interval when the root is the right endpoint. Thus our answer would be less accurate than it should be, and in any event Maple might spend a lot more time trying to find the root than it really should (as it had the root for quite a while but didn't know it!).

2. We really don't know how many times to run the loop (we guessed at most 100 times), but for all we know, it may need more.

A general procedure would help, one that we could use and reuse again and again for different function and different intervals (and one that wouldn't necessarily have the function at the left endpoint being the negative one). A procedure is a set of lines that carries out a single process based on whatever items are inputted (the latter are called the parameters of the procedure). Below is the corresponding procedure (which we'll call ivtRoot) that takes into account both of the comments above:

> ivtRoot:=proc(f,a,b,acc)
local i,leftEnd,rightEnd,c:
leftEnd:=a:
rightEnd:=b:
while (rightEnd-leftEnd > acc) do
c:=(rightEnd+leftEnd)/2:
if f(c) = 0 then
RETURN(`The root is exactly`,c);
else
print(c,f(c));
if sign(f(c)) = sign(f(leftEnd)) then
leftEnd:=c:
else
rightEnd:=c:
fi:
fi:
od:
print(`The root is approximately`,evalf(leftEnd));
end;

This routine takes as input a function , the endpoints of the interval and and an accuracy and returns an approximation to a root of within an accuracy of . Here are some more details:

The left hand endpoint, , is renamed leftEnd, and the right hand endpoint, , is renamed rightEnd. You need to do this as the left and right hand endpoints change during the course of the procedure, and you can't change the input parameters and and !

The procedure uses a "while" loop rather than a "for" loop as we don't know how many times to run the loop, and wish to so do as long as the length of the interval is greater than the accuracy we desire.

If we stumble upon the root (i.e. ) then we print out our answer. It is important to break out of the procedure, as otherwise, the loop would continue until the interval is small enough (even though we're done already!). The "RETURN" command takes you out of the procedure you're in and sends back its argument, i.e. what is enclosed in the parentheses (for more details, see you Maple manual or the on-line help).

It's impossible to know in general whether is positive or negative (although we know that will be the opposite, or we wouldn't be searching for the root between and ). The key thing is that and have opposite "signs", and that we always search for the root in an interval for which the values at the endpoints have opposite signs. Luckily, Maple has a built-in function, "sign", for checking the sign of a number.

Let's try out the procedure a few times, first on our original example.

> ivtRoot(f,1,2,0.00001);

Well, that was pretty quick, wasn't it? Let's try another example. Note that the function

> g:=x->cos(x)-x;

is positive at (as ) and negative at (as ); see the figure below as well.

> plot(g(x),x=0..Pi/2);

The IVT (and the plot!) tells us that g has a root between and . Let's use our procedure to find it, say to within .

> ivtRoot(g,0,Pi/2,0.001);

Error, (in ivtRoot) cannot evaluate boolean

Well, what's going on? Why do we get an error message? Well, it seems that Maple has difficulty thinking sometimes of some exact numbers as decimals when Maple has to compare them (as in inequalities) or determine their sign (which is again is an inequality to check). The problem here is . Let's fix the routine by forcing Maple to convert numbers to decimals via the "evalf" command (We'll retype it below, but in actual fact, you can scroll back to where the ivtRoot procedure is and make the changes; don't forget to run it again by hitting the carriage return key afterwards!).

> ivtRoot:=proc(f,a,b,acc)
local i,leftEnd,rightEnd,c:
leftEnd:=evalf(a):
rightEnd:=evalf(b):
while (rightEnd-leftEnd > acc) do
c:=(rightEnd+leftEnd)/2:
if f(c) = 0 then
RETURN(`The root is exactly`,c);
else
if sign(evalf(f(c))) = sign(evalf(f(leftEnd))) then
leftEnd:=c:
else
rightEnd:=c:
fi:
fi:
od:
print(`The root is approximately`,evalf(leftEnd));
end;

> ivtRoot(g,0,Pi/2,0.001);

Well, that's much better (and the answer corresponds well to where the graph of crosses the axis in the previous plot ).

The example of root finding above shows precisely how mathematics and programming (via Maple) can connect to solve a problem. This is typical of the approach we'll take throughout this book. You'll let the mathematics decide what you want to do, and then you'll ask Maple to help solve some of the parts.

[Aside: Of course, Maple has a built -in command to find the solutions to equations:

> fsolve(h=0,x,a..b);

will solve the equation for in the interval . If is a polynomial, it will find all the real roots in the interval (and you can even omit the interval to find all the real roots!), but for general functions h, the command will only return one root in the interval (of course, there may be more).

Is there any reason not use this command? Why did we bother writing our own commands and procedure for something Maple already does? Well, there are a number of reasons. First, we learned how to take a mathematical theorem like the IVT and use it to write a procedure to do something (like find a root). Using a built-in command like the "fsolve" command wouldn't have taught us anything (beyond how to use it). Secondly, we now know exactly how our procedure works, and why it works ; we have to trust Maple that the "fsolve" command really does what its supposed to do (of course, most times we'll do exactly that!). Thirdly, if we want to alter the procedure, we're free to do so, and we know what's going on inside. Finally, we did learn alot about new Maple commands and the Maple programming language, and these will serve us well in the coming sections.]

Example 2: Maximizing and Minimizing a Polynomial on a Closed Interval

We'll go through one more example of a typical mathematical problem.

For example,

What is the maximum and minimum of on the interval ?

Really, the best we can hope for are accurate estimates for the maximum and minimum, so that is what we're after.

Recall from your calculus course:

The Extreme Value Theorem (EVT) : Any continuous function on a closed and bounded interval has a maximum and minimum value.

As any polynomial is continuous everywhere, the problem certainly has a solution. The trick is to convert the mathematics to Maple.

Now the general plan from your calculus course was the following for finding the maximum and minimum of a continuous function on an interval :

Find the critical numbers for f that are in (a,b)

The maximum of f on [a,b] is the largest of the f values at the critical numbers and the endpoints, while the minimum of f on [a,b] is the smallest of the f values at the critical numbers and the endpoints

To find the critical numbers, we differentiate the function and find all the roots of the derivative in the interval. To do this we'll use Maple's built-in "fsolve" command, as it finds all the roots of a polynomial in an interval. We can then use Maple's "max" and "min" commands to pick off the largest and smallest values of the polynomial at endpoints and critical numbers.

First we'll let be the name of the polynomial function.

> p:=x->x^8-x^7-10*x^6+9*x^5+36*x^4-27*x^3-54*x^2+27*x+27;

Next, we'll find the derivative of and use Maple's fsolve command to find all its roots in .

> p1:=D(p);

> fsolve(p1(x),x,-2..3);

There are certainly a number of critical numbers in the sequence above. Let's give it a name, say cS.

> cS:=%;

Now we need apply the polynomial to each value in this list (as well as to the endpoints). We can get at the i-th item in the sequence cS by typing cS[i]. For example,

> cS[2];

Finally, we can calculate the maximum and minimum values of p on the interval:

> max(p(-2),p(3),p(cS[1]),p(cS[2]),p(cS[3]),p(cS[4]),p(cS[5]),p(cS[6]),p(cS[7]));

> min(p(-2),p(3),p(cS[1]),p(cS[2]),p(cS[3]),p(cS[4]),p(cS[5]),p(cS[6]),p(cS[7]));

That does it! The maximum and minimum of on the interval are and respectively. We'll do a quick check of the plot of p on the interval:

> plot(p(x),x=-2..3);

The maximum certainly does seem to happen a little over at the right endpoint, . The minimum is a little harder to see, so why don't we replot the graph from to (as we can certainly see that the function is negative near and positive to the right of (and thus the minimum occurs to the left of ).

> plot(p(x),x=-2..2);

Ah, we see that the minimum seems to occur around -1, and that it is about -8.

Now the problem that we solved was just one of a more general one:

Find the maximum and minimum of a polynomial on a closed interval.

We'll write a general procedure "maxMinPoly" to do the work; the procedure will take as input a polynomial function and the endpoints of the interval, and print out the maximum and minimum values. The fact that the functions we're considering are polynomials ensures that Maple's fsolve command (applied to another polynomial, namely the derivative of ) will find all the roots in the interval

> maxMinPoly:=proc(p,a,b)
local p1,critSeq,maxPoly,minPoly,i:
p1:=D(p):
critSeq:=fsolve(p1(x),x,a..b):
maxPoly:=max(p(a),evalf(p(b))):
minPoly:=min(p(a),evalf(p(b))):
for i from 1 to nops([critSeq]) do
maxPoly:=max(maxPoly,p(critSeq[i])):
minPoly:=min(minPoly,p(critSeq[i])):
od:
print(`The maximum of`,p(x),` on the interval `,[a,b],`is `,maxPoly);
print(`The minimum of`,p(x),` on the interval `,[a,b],`is `,minPoly);
end;

We need to talk about some of the commands. First, the critical numbers are found by the fsolve command (applied to the derivative of ); this sequence is placed into the variable "critSeq". To count the number of items in the sequence critSeq, we need to convert the sequence to a list (simply by enclosing it in square brackets) and using the "nops" command (this picks off the number of "pieces" in a list).

To get the maximum and minimum of the f values at the required points, we use Maple's "max" and "min" commands, which act on sequences of numbers. We'll only use it on pairs of numbers. What we'll do is first calculate the maximum and minimum of the polynomial at the endpoints of the interval. Then we'll suceesively look at the value of the polynomial at each critical number, and update the maximum (and minimum) found so far if the values of the polynomial at that critical number is larger (or smaller) than the maximum (respectively, minimum) found so far.

Let's try out our new maxMinPoly procedure on our particular problem:

> p:=x->x^8-x^7-10*x^6+9*x^5+36*x^4-27*x^3-54*x^2+27*x+27;

> maxMinPoly(p,-2,3);

[Aside: Rather than try out other examples for our maxMinPoly routine, we'll end this section by remarking that we can "tidy" up the maxMinPoly routine a little bit by using an advanced Maple command called "map". What "map" does is to apply a function to a list (or set). For example:

> lst:=[2,-3,5];

> map(ln,lst);

Now max and min operate on sequences, not lists. It is important to know that to remove the brackets from a list L, all you need to do is put "[]" after L, that is write L[].

> %[];

Here is how we might redo the procedure:

> maxMinPoly:=proc(p,a,b)
local p1,lst,maxPoly,minPoly,i:
p1:=D(p):
lst:=map(p,[fsolve(p1(x),x,a..b),a,b])[]:
maxPoly:=max(lst):
minPoly:=min(lst):
print(`The maximum of`,p(x),` on the interval `,[a,b],`is `,maxPoly);
print(`The minimum of`,p(x),` on the interval `,[a,b],`is `,minPoly);
end;

The command first makes a sequence out of the critical numbers and the endpoints and makes them into a list, then applies the polynomial function p to all of them, and then finally converts the list back to a sequence. Pretty complicated, isn't it?

Trying out the particular example again, we get

> maxMinPoly(p,-2,3);

That worked pretty well. If you're not happy with the complexity and brevity of the procedure above, you're certainly not alone. Feel free to write procedures that are easy to understand; writing "tighter" routines will come as you become more and more at home with Maple. For right now, stick with simple but understandable!]

Exercises

1. For each of the following functions and intervals , find a root of in the interval to within an accuracy of 0.00000001.

(a)

(b)

(c)

(d)

(e)

2. For each of the problems in question 1, explain clearly why they have exactly one root in the given interval.

3. For each of the following functions, find all the roots of the given function in the interval (use your ivtRoot procedure and some mathematics; don't use Maple "fsolve" command!):

(a)

(b)

(c)

(d)

(e)

4. Alter the procedure ivtRoot so that it checks the value of function at the endpoints. If there is a root at an endpoint, return the value(s), and if the function does not have opposite signs at the endpoints, and otherwise print s that the procedure cannot be used.

5. For each of the following polynomial functions and intervals , find accurate estimates for the maximum and minimum of on the interval.

(a)

(b)

(c)

6. For each of the following functions, find accurate estimates for the maximum and minimum on the given interval (hint: use a mix of both mathematics and Maple!).

(a)

(b)

(c)

7. Write a Maple procedure "rollesProc" that takes as input a continuous function and numbers and (with ) and prints "At most one root in the interval" if the derivative of has no roots in the interval, and prints "Unsure of the number of roots" otherwise (recall that Rolle's Theorem states that if the derivative of a function has roots in an interval , then can have no more than roots in ). Afterwards, test your procedure on some examples.

8. [Harder] Suppose we are considering the maximum and minimum of a continuous function (which is not necessarily a polynomial) on an interval , and we have determined (perhaps through Maple plots) that f has exactly one root in each of the subintervals ,É, . Write a Maple procedure maxMinFunc that takes as input , where intList = [a_1,a_2,É,a_2n], and returns (accurate estimates for) the maximum and minimum of on the interval . Afterwards, test your procedure on the problems in exercise 6.