SAMPLE EXERCISES USING MAPLE

N. P. STRICKLAND

Polynomials Exercise 1. Expand out the following products: (x 1)(1 + x) (x 1)(1 + x + x2 ) (x 1)(1 + x + x2 + x3 ) (x 1)(1 + x + x2 + x3 + x4 )

What would be the next two equations in this sequence? What is the general rule? (Try to write down your answer as a complete, self-contained statement that could be understood by someone who had not seen the calculation before.) Exercise 2. Consider the following polynomials 1 (x) = x 1 2 (x) = x + 1 3 (x) = x2 + x + 1 4 (x) = x2 + 1 6 (x) = x2 x + 1 12 (x) = x4 x2 + 1. Expand out the following products: 1 (x)2 (x) 1 (x)3 (x) 1 (x)2 (x)4 (x) 1 (x)2 (x)3 (x)6 (x). Can you see the pattern? Can you guess what is the corresponding equation involving 12 (x)? Check your guess using Maple. Remark 3. These functions are called cyclotomic polynomials; they are important in Number Theory, and have many interesting properties. There are various dierent ways to dene them, one of which is based on the pattern that you should have observed. Exercise 4. Factor the polynomials x2 1, x4 1, x8 1 and x16 1 using Maples factor() command. What is the pattern? What is the next thing in the sequence? Check your answer with Maple.
Date: November 21, 2003.
1

N. P. STRICKLAND

Exercise 5. Put u = a + b and v = ab. Expand out the following expressions: u2 2v u3 3uv u4 4u2 v + 2v 2 u5 5u3 v + 5uv 2 u6 6u4 v + 9u2 v 2 2v 3 . You should do the rst three by hand, and get Maple to do the remaining two. You should see an obvious pattern in your answers. (There is a pattern in the questions as well, but it is a lot harder to describe.) Graphs of powers Exercise 6. Plot the function y = x3 for 1 x 1, like this: > plot(x^3,x=-1..1); Now plot the functions y = x, x2 , ..., x6 together, like this: > plot({x,x^2,x^3,x^4,x^5,x^6},x=-1..1); Describe the pattern in words. What do you think that the graph of x100 looks like? How about x101 ? Check your answer like this: > plot({x,x^2,x^3,x^4,x^5,x^6,x^100,x^101},x=-1..1); Exercise 7. Plot the functions x5 and x1/5 together. In this case, it is convenient to ensure that the graph is square, with the same units on the two axes; for this we use the option scaling=constrained. > plot({x^5,x^(1/5)},x=0..1,scaling=constrained); What is the relationship between the two graphs? Try it again with the 5 replaced by 6. Then just plot the graph of x7 , and draw by hand the graph of x1/7 (on paper or in your imagination). Decaying oscillations Consider the function y = exp(t) sin(20t). We rst plot it to get a feel for what it looks like: > plot(exp(-t) * sin(20*t),t=0..10); The function oscillates rapidly, but also gets smaller and smaller, becoming very small when t gets large. This is the sort of response that you get if you pluck a guitar string; it vibrates, but the vibrations die away. You get similar behaviour with water waves, oscillations in electrical circuits, and almost any other physical system that you can think of. Of course, the exact formula will be a bit dierent in each case, depending on the size and frequency of the oscillations and the rate of decay. The most general type of function that we will consider is y = a exp(bt) sin(ct), where a, b and c are constants. What dierence do these constants make? To investigate, enter the following in Maple. > a := 2; b := 1.5; c := 20; > plot(a * exp(-b*t) * sin(c*t),t=0..5,y = -2..2);

SAMPLE EXERCISES USING MAPLE

Now go back to the rst line, change some of the numbers, press RETURN, then go to the plot line and press RETURN again to plot the new function. First change a several times, then change b several times, then c. You should conclude that a just aects the size of the oscillations (which is easy to see anyway). Increasing b makes the oscillations decay more quickly, and increasing c makes the oscillations more rapid. From now on, we want to treat a, b and c as symbols, not as particular numbers: > unassign(a,b,c); We next calculate the derivatives of y. We can write y in the form y = uv, where u = a exp(bt) and v = sin(ct). You should remember that d du = a exp(bt) = ab exp(bt) dt dt dv d = sin(ct) = c cos(ct). dt dt We also have the rule for dierentiating a product of two things: dv du d (uv) = u +v . dt dt dt This gives dy = ac exp(bt) cos(ct) ab exp(bt) sin(ct). dt You can check this in Maple as follows: y = > y := a*exp(-b*t)*sin(c*t); > diff(y,t); (As quite often happens, Maple gives the terms in a slightly unnatural order, but nonetheless you can check that tey match up with what we wrote above.) Next, you can compute the second derivative of y: diff(y,t,t); The answer contains three terms, all a bit like the ones we have seen before. There is in fact a precise relationship between y, y and y , given by the equation y + 2by + (b2 + c2 )y = 0. You can check this in Maple: expand(diff(y,t,t) + 2*b*diff(y,t) + (b^2+c^2) * y); This identity is a dierential equation for y. Physically speaking, our whole discussion is really the wrong way around. The physical theory gives the dierential equation that we ended with, and one has to solve the equation to nd the formula that we started with. Gaussians In this example, we consider variants of the function f (x) = exp(x2 ). We rst get the general picture: > plot(exp(-x^2),x=-5..5); We have a hump of height one, centred at x = 0. The function is nicely symmetrical. It seems to be zero when x is bigger than 3 or so, and is never negative. The curve is reasonably close to a triangle with vertices at x = 2 on the x-axis, and at y = 1 on the y-axis. The area of a triangle 1 is half the base length times the height, which in our case is 2 4 1 = 2, so the area under the curve is roughly 2. At x = 0 the curve is horizontal, meaning that the derivative f (0) is zero. The curve is at when x is large, so again it looks like f (x) = 0 when x is large. It is important to be able to make observations like these from the graph, but it is equally important to be able to prove them mathematically using formulae.

N. P. STRICKLAND

(1) Where is the hump, and how big is it? Clearly f (0) = exp(0) = e0 = 1. For any x = 0 2 we have x2 < 0 and so exp(x2 ) = 1/ex < 1. Thus, the maximum value is 1, and this value occurs when x = 0. Another approach is to say that the maximum is found where the derivative is zero; we will discuss the derivative later. (2) To say that the function is symmetrical around the y-axis just means that f (x) = f (x). This is easy to verify, because (x)2 = x2 and so (x)2 = x2 and so f (x) = exp((x)2 ) = exp(x2 ) = f (x). (3) Is it really true that f (x) = 0 when x 3? We have f (3) = e9 = 1/e9 . The number e9 is large (about 8000) but not innite, so 1/e9 cannot be zero. Take another look at the graph, concentrating on the region where 3 x 4: > plot(exp(-x^2),x=3..4); The values are small, but strictly positive. 2 (4) Is it true that f (x) is never negative? As e > 0 we have ex > 0 for all x and so 2 f (x) = 1/ex > 0 as claimed. (5) What is the area under the curve? Recall that the area under a curve is just the denition (or at least one of the possible denitions) of integration. The area that we want is given by:

f (x)dx

What is the integral of f (x)? Is it perhaps something like exp(x2 )/(2x)? If you think you know an integral but you are not sure, there is an easy way to check: just dierentiate, and see if you get back to the function you started with. We can do this by hand, or ask Maple: > diff(exp(-x^2)/(-2*x),x); The answer is not equal to f (x), so our integration was incorrect. We can just ask Maple to do the integration instead: > int(exp(-x^2),x); The answer involves something called erf(x), which you probably have not seen before. This is not so surprising; Maple knows a long list of strange functions, many of them unknown even to most professional mathematicians (although erf itself is only moderately obscure.) If you are curious, try > help(erf); In any case, we do not need the general value of the integral, just the integral from to , and it turns out that this can be given in familiar terms: > int(exp(-x^2),x=-infinity..infinity); The answer is really quite amazing: where on earth did the come from? Let us check the numerical value: > evalf(sqrt(Pi));

SAMPLE EXERCISES USING MAPLE

The answer is about 1.77; recall that our crude estimate with a triangle was 2, which is not great, but not too bad either. (6) We next consider the derivative of our function y = f (x). This can be calculated with the chain rule: we put v = x2 , so y = exp(v) and so dy = exp(v) dv dv = 2x dx dy dy dv f (x) = = = 2x exp(v) dx dv dx = 2x exp(x2 ). As exp(x2 ) is never zero, we see that f (x) is only zero when x = 0. This is another way to see that the maximum is at x = 0. You may be aware that the function f (x) (or related functions such as a exp((x b)2 /c)) are very important in statistics, as the probability distributions for random variables. It is less well-known that these functions also occur in the theory of temperature: if you heat a long bar in the middle and then let the heat spread out for a while, then the temperature prole will be given by a function like f (x). Of course, the prole will change over time: as the heat spreads out, the hump will get lower and broader. If we take the time dependence into account, the temperature is actually given by a function like u(x, t) = exp(x2 /4t)/ t. Let us plot this function, to check that the formula seems reasonable: > u := ((x,t) -> exp(-x^2/(4*t))/sqrt(t)); > plots[animate](u(x,t),x=-50..50,t=0.5..100); You will see a picture of the hump as it is at time t = 0.5. If you click on the plot, a box will appear around it and some new controls will appear at the top of your Maple window, like the controls on a CD player. Click on the play button to see how the plot changes over time; it spreads downwards and outwards, as expected. Just as in the case of damped oscillations, the formula for u really comes from solving a differential equation. As this is a course in mathematics rather than physics, we allow ourselves to work backwards and just check the dierential equation as a mathematical exercise. We can dierentiate u just as before, by putting v = x2 /4t (so u = exp(v)/ t) and using the chain rule. We have dv d x2 x 2x = = . = dx dx 4t 4t 2t du du dv exp(v) x x exp(x2 /4t) = . = = dx dv dx 2t3/2 t 2t Here we have just carried the t through treating it as though it were a constant. This just means that we are thinking about how u varies along the bar at a particular moment in time, and not worrying about how u changes with t. This is called a partial derivative, because we are only measuring part of the variation of u. It is usual to change notation slightly to indicate this, and write u/x instead of du/dx to indicate this. Thus, our conclusion is that x exp(x2 /4t) u = . x 2t3/2 We can instead consider a xed point on the bar, and consider how the temperature at that point varies over time. This means that we treat x as a constant, and dierentiate with respect to t: > diff(u,t);

N. P. STRICKLAND

I claim that u satises the dierential equation u 2 u 2 = 0. t t This is easy to check with Maple: > simplify(diff(u,t) - diff(u,x,x)); Can you verify this by hand? Mobius functions A Mbius transformation is a function of the form o ax + b f (x) = cx + d where a, b, c and d are nonzero numbers with ad > bc. (Here we have slightly changed the proper denition, to make our lives easier.) A typical example is the function f (x) = (x 1)/(x + 1) (here a = 1, b = 1, c = 1 and d = 1). Let us look at the graph: > f := (x) -> (x-1)/(x+1); > plot(f(x),x=-10..10,y=-10..10); Something funny is happening at x = 1. That seems reasonable, because the formula gives f (1) = 2/0, which is undened. You might want to say that 2/0 = , but all statements about innity are fraught with diculty and subtle pitfalls. It is hard even to say things that are truly meaningful, and harder still to say things that are actually correct. The issues involved are in fact well understood (this is one of the great triumphs of 19th century mathematics) but we will leave them for future courses (starting with PMA113, Introduction to Analysis). For the moment, we will simply say that f (1) is undened. When x is very large, the curve appears to atten o, but it is not too easy to see the level that it is tending to. To work this out, we rearrange the formula a little: (x + 1) 2 2 x1 = =1 . x+1 x+1 x+1 When x is large (either positive or negative) then 2/(x + 1) will be very small, so f (x) will be close to 1. We can try this numerically, either directly: f (x) = > f(1000.0); or by plotting the graph for large values of x: > plot(f(x),x=1000..2000); Alternatively, we can ask Maple to work out the limit symbolically: > limit(f(x),x=infinity); It is also interesting to work out the derivative of f (x). We can write f (x) = u(x)/v(x), where u(x) = x 1 and v(x) = x + 1. The derivative is then given by the quotient rule: f (x) = u (x)v(x) v (x)u(x) 1.(x + 1) 1.(x 1) 2 = = v(x)2 (x + 1)2 (x + 1)2

If you ask Maple to do this, you will get a more complicated expression, and you will need to simplify it explicitly to get our answer: > diff(f(x),x); > simplify(%);

SAMPLE EXERCISES USING MAPLE

Here is a puzzle. As (x + 1)2 is never negative, the derivative as always positive, so the function slopes upwards and is always increasing. Out at the left hand end of the graph, x is a large, negative number, so f (x) is approximately 1. The function increases from there, so f (x) should always be at least 1. But f (0) = (0 1)/(0 + 1) = 1. What is going on? The answer is clear from the pictures that we drew earlier. The graph does indeed slope upwards everywhere, but it also has a massive jump at x = 1. The idea that a function with positive derivative is increasing relies on there being no jumps. You will study such questions in more detail in PMA113, but for the moment, you should just keep an eye open for the kinds of things that can happen when innities crop up. Now suppose we have a Mbius transformation g(x) = (ax + b)/(cx + d), but we do not know o what a, b, c and d are. Then g(0) = b/d and we do not know what b/d is, so we cannot even begin to plot the graph. Instead, we must work with formulae to understand the features of the function. (1) When is g(x) undened? A problem occurs if we try to divide by zero, which happens when cx + d = 0, or equivalently x = d/c. The formula g(x) = (ax + b)/(cx + d) is ne for all other values of x. (2) What happens when x is very large? You will see more about this in PMA113, but this example is quite easy to understand directly. The trick is to rewrite the formula slightly, by dividing the top and bottom of the fraction by x: g(x) = a + b/x ax + b = . cx + d c + d/x

When x is very large, b/x and d/x will be very small so we can neglect them and we see that g(x) a/c. (3) What is the derivative? We can again use the quotient rule: g (x) = acx + ad acx bd ad bc a(cx + d) c(ax + b) = = . 2 2 (cx + d) (cx + d) (cx + d)2

We assumed at the beginning of this section that ad > bc, so adbc is positive, and (cx+d)2 is nonnegative (because it is the square of something). We have the usual problem when x = d/c because then cx + d = 0 and our formula g (x) = (ad bc)/(cx + d)2 involves division by zero and so is not dened. For all other values of x, we see that g (x) > 0 and so g(x) is increasing. We conclude that the picture is the same as for the special case that we studied before. The function starts o at approximately a/c when x is very large and negative, then increases to as we approach x = d/c, jumps instantly down to and then increases back to a/c as x becomes large and positive. Star Trek In this section we will measure speed in terms of warp factors, as in Star Trek. Warp factor 1 is the speed of light; 70 mph is warp factor 0.0000001. Sadly for Star Trek, warp factors greater than one are forbidden by Einsteins theory of relativity, and do not occur in the real world. Remark 8. Warp factor 1 is about a foot per nanosecond. A nanosecond is a billionth of a second, which is the clock cycle time of a computer chip running at 1 GHz, which is fairly typical in 2001. Suppose we have an object moving at warp factor v. If v is close to 1 (ie the object is approaching the speed of light), then all sorts of strange relativistic eects occur: time slows down, lengths contract, light changes colour, and so on. The strength of these eects is given by the function g(v) = (1 v 2 )1/2 = 1 . 1 v2

In this section, we investigate this function as a purely mathematical exercise. First, we draw the graph:

N. P. STRICKLAND

> g := (v) -> (1-v^2)^(-1/2); > plot(g(v),v=-2..2,y=0..10); (1) First, note that the graph is symmetrical, ie g(v) = g(v). This is obvious from the formula, because (v)2 = v 2 , and is also physically reasonable: the strength of relativistic eects does not depend on the direction that you are moving in. (2) Next, we have g(0) = 1. This means that if you are not moving (ie v = 0) then time slows down by a factor of 1 (ie it does not slow down at all). This seems very reasonable. (3) As v approaches 1 we see that 1 v 2 approaches 0, so 1 v 2 also approaches 0, so g(v) tends to innity, and the relativistic eects become extremely strong. (4) What has happened to the graph for v > 1? If v > 1 then v 2 > 1 so 1 v 2 < 0. We cannot take the square root of a negative number, so g(v) is undened. This is reasonable, given Einsteins dictum that speeds greater than warp 1 are impossible. (If you know about complex numbers, you may want to say that g(v) is dened but imaginary. This is mathematically sensible but not physically meaningful.) (5) We next look at g(v) when v is small. Of course we then have g(v) g(0) = 1, but we want something a bit more accurate than that. We will use the approximation given by Taylors theorem: g(v) g(0) + g (0)v + g (0)v 2 /2 Weve seen that g(0) = 0. The function has a minimum at v = 0, so we must have g (0) = 0, but we will check this algebraically anyway. We rst work out the relevant derivatives: > diff(g(v),v); > diff(g(v),v,v); We then work out their values at v = 0. > eval(diff(g(v),v),v=0); > eval(diff(g(v),v,v),v=0); Putting this back into Taylors formula gives g(v) 1 + v 2 /2.

We could actually have asked Maple this more directly: > series(g(v),v,3); The motorway speed limit is approximately v = 107 , so g(v) 1 + 5 1015 . Your car clock is slowed by a factor of g(v), which works out to about 0.15 microseconds per year. This kind of eect is measurable with atomic clocks. Absolute zero There is a certain temperature called absolute zero (about 273 degrees centigrade) that is the coldest that anything can possibly get. (Warmth is caused by random movement of molecules; at absolute zero they are not moving at all, so it is not possible to cool down any further.) Well use a scale of temperature in which t = 0 corresponds to absolute zero. In this context, the following function comes up regularly: f (t) = exp(1/t) = 1/ exp(1/t). It has some interesting mathematical properties, which we will now investigate. We start by plotting the graph: > f := (t) -> exp(-1/t); > plot(f(t),t=0..1);

SAMPLE EXERCISES USING MAPLE

(We have only drawn the curve for t > 0, because negative temperatures do not make sense.) The main feature of interest is that the curve is very at near t = 0, which is a reection of the rather strange nature of absolute zero. To see this more clearly, plot the graph for 0 t 0.1. How can we understand this mathematically? (1) Firstly, what is f (0)? We cannot just put t = 0 in the formula, because then we would have a term 1/0 which is undened. Instead we need to look at what happens as t approaches 0. For this, we need to recall the formula for exp(x): exp(x) = 1 + x + so exp(1/t) = 1 + t1 + t2 /2 + t3 /6 + . As t is positive, all the terms are positive, so we have exp(1/t) = t1 + other positive terms t1 . It follows that 0 f (t) = 1/ exp(1/t) t. If you like, you can check this by plotting f (t) and t together: > plot({f(t),t},t=0..1); In any case, as 0 f (t) t, it is clear that f (t) goes to zero as t goes to zero. (2) Next, the curve seems to be at at x = 0, which means that f (0) should be 0. We can dierentiate f by the chain rule: f (t) = exp(1/t) d dt 1 t = t2 exp(1/t) = t2 f (t). x3 x2 + + , 2 6

What happens as t approaches zero? Certainly t2 tends rapidly to innity, but it is multiplied by f (t), which (as we saw above) tends to 0. Naively we have something like 0, which is horribly undened. We need to be more careful. The solution is to reuse the method of the previous section. From the formula for exp(1/t), we see that exp(1/t) = It follows that 0 f (t) = 1/ exp(1/t) 6t3 , and so 0 t2 f (t) 6t. If you like, you can check this by plotting t2 f (t) and 6t together: > plot({f(t)/t^2,6*t},t=0..1); In any case, as 0 f (t) = t2 f (t) 6t, it is clear that f (t) goes to zero as t goes to zero, as expected. (3) In the last section, we had a competition between f (t) and t2 . The factor f (t) was tending to zero, and t2 was tending to innity, and f (t) won because the product t2 f (t) turned out to tend to zero. What would happen if we replaced t2 by something that tended to innity much faster, like t5 or t10 ? We can investigate by plotting the graphs of tk f (t) for various values of k: > plot(f(t)/t^4,t=0..1); > plot(f(t)/t^6,t=0..1); > plot(f(t)/t^10,t=0..1); > plot(f(t)/t^20,t=0..1); t3 t3 1 + other positive terms = 3 6 6 6t

10

N. P. STRICKLAND

In each case, the graph becomes very large at around t = 1/k, but it then gets small again as t decreases, and tends to 0 at t = 0: the f (t) factor denitely wins. We can verify this mathematically by the same method as above: in the full series for exp(1/t), one of the terms is tk1 /(k + 1)!, so exp(1/t) It follows that 0 f (t) = 1/ exp(1/t) (k + 1)!tk+1 so 0 tk f (t) (k + 1)!t so tk f (t) tends to 0 as t tends to 0. This is the precise sense in which the graph is very at at t = 0; it tends to zero more rapidly than tk for any k. Remark 9. The statement that all movement stops at absolute zero ignores some important complications from quantum mechanics; but we will not go into details here. Matrices Exercise 10. (a) Consider the following matrix 1 1 A=a b a2 b2 (called a Vandermonde matrix): 1 c c2 1 t1k . = (k + 1)! (k + 1)!tk+1

Write down the transpose matrix AT , then calculate AAT by hand. Can you see the pattern? You may want to introduce some new notation, so that you do not have to write the same expression several times over. (b) Now consider the matrix 1 1 1 1 a b c d B= 2 2 2 a b c d2 a3 b3 c3 d3 Can you guess what BB T looks like? Check your guess with Maple. (c) Use Maple to calculate the determinant of B, and factorize it. (d) What is the determinant of B when a = b? Give an answer based on just looking at the matrix, and another answer based on part (c). Exercise 11. Consider the following matrix: a b c A = b c a c a b (a) Calculate det(A) by hand. (b) Ask Maple to factor det(A). (c) What is det(A) when a + b + c = 0? Give an answer based on part (b), and another answer based on row operations. Exercise 12. Consider the following matrix: 0 a b A = a 0 c b c 0

(a) What is A + AT ? (b) What is det(A)?

SAMPLE EXERCISES USING MAPLE

11

(c) The answer to (b) actually follows from the answer to (a) can you see why? (This is quite challenging) Exercise 13. Consider the following matrix (called a Jordan block): 1/a 1 0 0 0 1/a 1 0 A= 0 0 1/a 1 0 0 0 1/a Ask Maple to nd A1 , and then check by hand that AA1 = I and A1 A = I. Calculate det(A) by hand. Exercise 14. Consider the following matrix: 1 0 1 a 1 c A= b a d 0 b 0 Check that det(A) = (b d)2 + (a c)(ad bc). Exercise 15. Find the inverse of the following 1 1 A= 1 1 matrix: 1 1 1 0 1 1 0 0 1 0 0 0

0 1 c d

Exercise 16. Consider the following matrix: 1 2 1 1 3 3 A= 1 4 6 1 5 10

0 1 4 10

0 0 1 5

0 0 0 1

Do you recognize the entries? Row-reduce the matrix. Exercise 17. Suppose that x2 + y 2 + z 2 = 1 (make sure that you use this fact to simplify all your answers). Consider the matrix 2 x xy xz A = xy y 2 yz xz yz z 2 (a) (b) (c) (d) What What What What is is is is det(A)? trace(A)? AT ? A2 ? Quaternions Dene a b c b a d Q(a, b, c, d) = c d a d c b d c b a

(Matrices like this are called quaternion matrices.) We can teach the rule to Maple as follows: > Q := (a,b,c,d) -> matrix([[a,-b,-c,-d],[b,a,-d,c],[c,d,a,-b],[d,-c,b,a]]);

12

N. P. STRICKLAND

Get Maple to calculate the product P = Q(a, b, c, d)Q(w, x, y, z). Let p, q, r and s be the entries in the rst column of P . If you check carefully, you will see that the entry at the top of the second column is the same as q. Similarly, every other entry in P is either p, q, r or s, so the whole matrix can be written in terms of p, q, r and s. Do this, and thus check that P = Q(p, q, r, s). Euler angles Consider the matrix cos() cos() sin() cos() sin() A = sin() cos() cos() sin() sin() sin() 0 cos() (a) Calculate det(A), using the identities sin()2 + cos()2 = sin()2 + cos()2 = 1 to simplify your answer. (b) Calculate the inverse of A by the cofactor method. Check that A1 = AT . Quadratics Let a, b and c be real numbers, and put f (t) = t2 (a + c)t + (ac b2 ). This is a quadratic function, so it might have no real roots, or one real root, or two real roots, and if it has any real roots, they might be positive or negative or zero. For functions f (t) as above, it turns out that not all of these possibilities can occur. (a) Enter the following in Maple: > f := (t) -> t^2 - (a + c) * t + (a * c - b^2); Now, if you want to see the graph of f (t) when a = 1, b = 2 and c = 3 you can do > a := 1; b := 2; c := -3; > plot(f(t),t=-10..10); You will see that there are two real roots, one negative and one positive. Try some other values of a, b and c, and see what happens. (You may also want to plot a dierent range of values instead of 10 t 10.) Find the number t0 where f (t0 ) = 0. Where is t0 in the above pictures? Show that f (t0 ) = (a c)2 /4 b2 . When is this equal to zero? Show that if a = c and b = 0, then f (t) has precisely one real root, but for all other values of a, b and c, the function f (t) has two dierent real roots. Suppose that a > 0 and ac b2 > 0. Show that c > 0 and thus that t0 > 0. Observe that f (0) = ac b2 > 0. By looking through the pictures for cases that are compatible with these facts, show that all the roots of f (t) are real and positive. Conversely, suppose that all the roots of f (t) are real and positive. This means that there are numbers , > 0 such that f (t) = (t )(t ). Show that ac b2 > 0 and a + c > 0. Deduce that ac > 0, and say what this means about the signs of a and c. Given what you know about the sign of a + c, deduce that a > 0. This gives the converse of part (e).

(b) (c) (d) (e)

(f)

b Remark 18. This exercise is really about the matrix A = a c . Note that A is symmetric, ie b T A = A. Note also that a is the top left entry in A, and that ac b2 is the determinant, and that f (t) is the characteristic polynomial. Thus parts (e) and (f) say that the eigenvalues of A are positive i the determinant and the top left entry are positive. This can be generalized as follows. Let A be a real symmetric n n matrix; then all eigenvalues of A are real. Now let Ak be

SAMPLE EXERCISES USING MAPLE

13

the submatrix formed by the rst k rows of the rst k columns, so Ak is a real symmetric k k matrix. Then it turns out that all eigenvalues of A are positive i det(Ak ) > 0 for all k. It would be very painful to prove this directly, as we did in the case n = 2. Instead, we use more abstract methods, covered in the Level 2 Linear Mathematics courses.

Sequences Exercise 19. Fix two positive real numbers a and b. Put xn = (an+1 + bn+1 )/(an + bn ) yn = (an b + abn )/(an + bn ). We can teach this to Maple, and try it out when a = 10 and b = 20 as follows: > > > > > > x := (n) -> (a^(n+1) + b^(n+1))/(a^n + b^n); y := (n) -> (a^n * b + a * b^n)/(a^n + b^n); a := 10.0; b := 20.0; seq(x(n),n=1..50); seq(y(n),n=1..50);

What do you think are the limiting values? Experiment with other values of a and b if you wish, work out the general rule, and then prove it using the Algebra of Limits.

Wilsons theorem Investigate the numbers n! (mod p), where p is prime. (a) Start with a reasonably large prime p, say p = 19. Choose a random value of n, say n = 10, and then work out n! (mod p) like this: > modp(10!,19); Repeat this for some more values of n. (b) Now do the same thing more systematically, by listing the sequence of values of n! (mod 19) for n = 1, . . . , 100: > seq(modp(n!,19),n=1..100); There is a very obvious pattern when n is large. When exactly does this pattern start? What do you see shortly before the pattern starts? (c) Repeat part (b) for some other prime numbers in place of 19. You can get a list of the rst 20 prime numbers like this: > seq(ithprime(i),i=1..20); Write down your conclusions. Your answer should consist of complete sentences, and should be essentially self-contained, and comprehensible to someone who has not read the question. You should distinguish clearly between things that you can prove and things that you are merely guessing. (It does not matter if you cannot prove anything, provided that you say so.) (d) Explain the pattern for large n. You should have noticed something that happens just before the large-n pattern starts. This phenomenon is explained by Wilsons Theorem, which will be proved a little later in the course.

14

N. P. STRICKLAND

Binomial coefficients Let b(n) be the greatest common divisor of the numbers teach Maple this denition? (1) The binomial coecient n k can be expressed as binomial(n,k). n k with 0 < k < n can be expressed as n k for 0 < k < n. How can we

(2) The sequence of binomial coecients

seq(binomial(n,k),k=1..n-1). (3) The greatest common divisor of a sequence of integers is given by the function igcd. (4) Thus, the required denition is > b := (n) -> igcd(seq(binomial(n,k),k=1..n-1)); We will now investigate the behaviour of this function. (a) List the values of b(n) for n = 1, . . . , 100: > seq(b(n),n=1..100); You should observe than b(n) is usually equal to 1. Our next task is to understand the special values of n where b(n) = 1. (b) Roughly in the middle of the list in (a), you should see the number 53. Unfortunately, it is not immediately clear what was the value of n that gave the answer b(n) = 53. To cure this, try > seq([n,b(n)],n=1..100); In the middle of the resulting list, we see . . . , [52, 1], [53, 53], [54, 1], . . . . This means that b(52) = 1, b(53) = 53, and b(54) = 1. (c) Look through the list in (b) for some numbers n where b(n) = 1. Factorize these numbers n (you can use the ifactor function for this). What do you notice? (d) To be more systematic, we can get Maple to do step (c) for us: > B := (n) -> (b(n) <> 1); > L := select(B,[seq(n,n=1..100)]); This sets L to be the list of all integers n in the range n = 1, . . . , 100 where b(n) = 1. We can factorize all these integers n as follows: > map(ifactor,L); What do you notice? Your answer should consist of complete sentences, and should be essentially self-contained, and comprehensible to someone who has not read the question. You should distinguish clearly between things that you can prove and things that you are merely guessing. (It does not matter if you cannot prove anything, provided that you say so.) Fermats little theorem In this exercise we investigate the numbers an (mod p), where p is prime and 0 < a < p. (a) Take p = 7 and a = 3. Calculate 3n for n = 0, . . . , 100: > seq(modp(3^n,7),n=0..100); Examine the resulting sequence carefully; you should notice that it repeats itself. How long is the cycle? What do we get at the beginning of each cycle?

SAMPLE EXERCISES USING MAPLE

15

(b) Keeping p = 7, repeat part (a) for a = 2, then for all other values of a in the range a = 1, . . . , 6. (c) Now repeat everything for a couple of other primes. What do you notice? Power sums For any prime p and any integer k 0 we put f (p, k) = 1k + 2k + . . . + (p 1)k (mod p).

We can teach this denition to Maple as follows: > f := (p,k) -> modp(sum(a^k,a=1..p-1),p); Take p = 7 and calculate f (p, k) for k = 0, . . . , 100. What do you notice? Try the same thing for some other primes in place of p = 7. Write down what you observe, as a complete and self-contained sentence or set of sentences. Dirichlets theorem In this exercise we investigate the numbers p (mod 12) for all primes p. We will need to use some statistical functions, so we load in Maples statistical packages: > with(stats); > with(transform); (a) We can nd all prime numbers using the function ithprime(); for example, the 20th prime number is ithprime(20), which is 71. (b) Set P to be the sequence of the rst hundred prime numbers, like this: > P := seq(ithprime(i),i=1..100); Now set M to be the sequence of these primes mod 12: > M := seq(modp(ithprime(i),12),i=1..100); Of course, because we are working mod 12, this sequence contains only numbers in the range 0, . . . , 11. For each of the numbers in this range, say whether (i) it occurs repeatedly in M ; or (ii) it occurs only once in M ; or (iii) it does not occur in M at all. There is a simple rule in terms of prime factors that says which numbers occur repeatedly; can you see it? (c) Explain why 0 does not occur in the sequence M . (d) Explain why 4 not occur. (e) We next look at the numbers that occur repeatedly, and check how many times they occur: > tally(M); The entry Weight(1,22) in the resulting list means that 1 occurs 22 times in M . The entry 2 just means that 2 occurs precisely once. You should notice that the entries that occur repeatedly occur roughly the same number of times each. The precise and general version of this statement is called Dirichlets theorem on primes in arithmetic progression; the proof is hard, but it may be covered in a level four course. With a little work, we can display this graphically. (Unfortunately the peculiarities of Maples statistical plotting functions make this harder than it should be.) > ranges := [seq(i-0.5..i+0.5,i=0..11)]: > T := tallyinto(M,ranges): > histogram(T,area=1,numbars=12,xtickmarks=12);

16

N. P. STRICKLAND

To get a clearer picture, you can repeat the process with the rst thousand (or even ten thousand) primes. (f) Repeat parts (a), (b) and (e) working modulo 15 instead of modulo 12. Before you start, try to predict what will happen in as much detail as possible.