Chapter9.mws

• Algebra on matrices filled with expressions

---

Chapter 9: Symbolic algebra

Debugging

This is what Maple was actually designed to do: symbolic mathematics. We have already seen some of its ability to do this as we have taken derivatives, integrated, and solved differential equations. I have placed this material in the last chapter because it's almost an art form and you need to be pretty comfortable with Maple before you tackle it. You can use this chapter as a source for what to try when Maple gives you answers that are correct, but not in the form you want. It may take many months to get good at this, but it is worth the effort because Maple will not make minus sign mistakes or lose factors of 2. The trick, however, is to coax it to combine precision with a sense of style and good taste. This is done via a partnership: Maple supplies the precision; you supply the style and taste. I will give you just enough of an introduction to get you started. After that it's up to you and online help.

copy and paste:

Maple almost never gives you exactly what you want, so you will still need to do some algebra by hand. But the editing capability of Maple using the mouse, copy ( ctrl-c ) and paste ( ctrl-v ) often makes it possible for you to stay in the Maple worksheet and do algebra by hand.

expand:

Here are some examples of what expand does.

> expand((x-1)^2*(x+2));

> expand(cos(x+y));

> expand(exp(x-y));

> assume(x,real);expand(ln(x/(1+x^2)));

> expand(tan(x-y));

This command is half of Maple's all-purpose algebra simplification duo-- simplify(expand(%)); . Always try this combination first when you get an algebraic expression that needs help.

combine:

This is the flip side of expand . The two of them together know some trig identities, but not as many as they should.

> restart;

> expand(cos(x+y));

> combine(%);

> expand(exp(x-y));

> combine(%);

> assume(x,real);expand(ln(x/(1+x^2)));

> combine(%);

> restart;

> expand(tan(x-y));

> combine(%);

Hmm---and sometimes it doesn't work. Let's try expand and simplify and see what happens

> expand(%);simplify(%);

Well, sometimes you just can't get where you want to go.

>

simplify:

Maple tries to read your mind with this command, but doesn't always succeed. You usually have to be creative with the other commands discussed here to get exactly what you want. A good thing to try is a combination of expand , simplify , and combine . And don't forget that you can change the forms of expressions by using copy and paste.

> simplify(cos(x)^2+sin(x)^2);

> simplify((x^2-2*x+1)/(x-1));

> simplify(cos(x)^2-sin(x)^2);

> simplify(2*sin(x)*cos(x),trig);

Simplify doesn't know it, but these last two can be rewritten as and . Simplify doesn't give these to you, but for some reason combine does

> combine(cos(x)^2-sin(x)^2);combine(2*sin(x)*cos(x),trig);

So you will have to be creative and play around with these commands to master them.

simplify with side-relations:

You can give simplify some extra help by using side-relations (see ?simplify/siderels ). The idea is to tell Maple what identities you want to use to put your expression into the form you want. For example, suppose you wanted to put the expression into a new form involving instead of by using the identity . You could do this

> f:=cos(x)^4-3+cos(x)^2;

> simplify(f,{cos(x)^2 = (1+cos(2*x))/2});

(Notice that the side-relations must be inside curly brackets.) You may think that this is just a substitution (see the subs command, discussed below), but it isn't the same:

> subs(cos(x)^2=(1+cos(2*x))/2,f);

Using side-relations is like an intelligent substitution where Maple tries to change everything it can using the identity you gave it. The side-relations go in curly brackets, as shown in the example above, and several identities can be specified at once using the syntax simplify(expression,{side1,side2,side3});

> f:=cos(x)^4+1-cos(x)^2-sin(x)^2+3*sin(x)^6;

> s1:=cos(x)^2=(1+cos(2*x))/2;

> s2:=sin(x)^2=(1-cos(2*x))/2;

> simplify(f,{s1,s2});



sort:

Sort puts expressions, especially polynomials, in order from high powers down to low

> sort(2-3*x^2-x+x^4);

> sort(1-cos(x)^2+cos(x));

factor:

Factor tries to write polynomials as products of terms

> factor(x^2-6*x+8);

collect:

This is sort of a half-way factor command. Perhaps you see a common factor within a larger expression and you would like to collect together all terms that contain it, like the term below. Then do this

> F:=x*ln(y)-5*x + 3*ln(y)+5*exp(y);

> collect(F,ln(y));

And if you want to collect more than one factor at once, give collect a list of factors, like this

> F:=a*x*sin(y)+b*a*x/cos(y);

> collect(F,[x,a]);

solve:

We have seen this before; it solves an equation for the specified variable. If the variable to solve for is obvious because there is only one, you don't have to give it; and if you just give solve an expression without an = sign, it solves it by setting it equal to zero.

> solve(x^2-6*x+8);

> solve(ln(y/x)=3,y);solve(ln(y/x)=3,x);

> solve(1-tan(x)^2=-2,x);

Can you spot the difference between this next one and the previous one?

> solve(1-tan(x)^2=-2.,x);

If you give Maple an equation that has no solution, it just hands it back to you, or gives you nothing.

> plot(1-tan(x)^2,x=-Pi..Pi,-3..3);

> solve(1-tan(x)^2=2,x);

As you can see from the plot the function can never be 2, so Maple can't solve the equation.

>

normal:

The normal form for rational functions is to have a common denominator

> normal(1/x + 1/(x+1));

And what happens if you try to expand this result?

> expand(%);

It's expanded, but not back to where it came from. To get it back to where it started you have to be more specific and use the convert command, which is coming up.

numer and denom:

Maple will let you refer directly to the numerator or denominator of a fractional expression, like this

> denom((1+x)/(3-x^2));

> numer((1+x)/(3-x^2));

but with perhaps not exactly the form you expected. You actually have much more control over parts of expressions like this by using the Maple editor's copy and paste ability:

> (1+x)/(3-x^2);

Use the mouse to highlight the numerator of the expression above, copy it ( ctrl-c ), then paste ( ctrl-v ) in this execution group

>

op:

This command picks out pieces of expressions. Its syntax is op(n,expr) , which picks out the nth piece of the expression expr. Just experiment with it whenever you need to take something apart.

> kk:=(1+x^2)/(sin(x)+tan(x))*cosh(x)/(1+tanh(x));

> op(1,kk);op(2,kk);op(3,kk);op(4,kk);

>

subs:

This is the substitute-one-variable-for-another command

> Eq1:=cos(x)^2-5*cos(x)+1;

> subs(cos(x)=y,Eq1);

> Eq2:=exp(2*x)-t^2=x*t;

> subs(x=xi,t=tau,Eq2);

As you use subs you may become frustrated because it seems stupid, like this:

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

> subs(x^2=y,f);

Why couldn't it recognize that ? Well, subs is just very literal, not algebraic. There is, however, an algebraic form of subs called algsubs which is much smarter. You will probably have lots better luck with it because it recognizes that your substitution was really an substitution.

> algsubs(x^2=y,f);

>

convert:

Finally, there is the convert command, which is actually many commands depending on what kind of conversion you want to do. Use ?convert to see all of the possibilities. When you get stuck and the commands discussed above don't do the job, try looking through the list of convert options in online help. Here are a few of what I think are useful or interesting options.

Put rational functions in partial fraction form with the parfrac option

> convert((2*x+1)/(x*(x+1)),parfrac,x);

Convert rational functions to continued fraction form with confrac

> convert((2*x+1)/(x*(x+1)),confrac,x);

Express functions in exponential form with exp

> convert(sin(x),exp,x);

> convert(sinh(x),exp,x);

Convert expressions with sines and cosines into tangent form with tan

> restart;

> a:=convert((2+cos(x)^2)/(2-sin(x)^2),tan);

> a:=simplify(a);

Convert expressions with tangents into sine-cosine form with sincos

> convert(a,sincos);

> simplify(%);

Convert English units to metric ones with convert(,metric)

> convert(1.*acre,metric);convert(1.*MPH,metric);

To see what English quantities can be converted to metric try

> ?metric

Convert expressions containing abs , Heaviside , signum , etc., into piecewise functions

> F:=abs(x);

> F:=convert(F,piecewise,x);

For instance, convert will give you a nice definition of Maple's signum function

> G:=signum(x);

> G:=convert(G,piecewise,x);

Convert complex expressions to polar form

> convert(x+I*y,polar);

Now this is probably not exactly what you had in mind. Fortunately Maple's map command can put it in a more useful form. The command map tells Maple to apply the procedure in the first slot of map(a,b) to each of the operands of the object in the second slot, in this case the magnitude and argument expressions, The conversion procedure is evalc, so we would do this

> map(evalc,%);

>

Map is quite powerful; see ?map for details.

ratpoly:

The ratpoly option converts a series into a rational function. This is often a very convenient and compact way to get an approximate expression for a complicated function. The syntax is

convert(series,ratpoly,numdeg,dendeg)

where numdeg is the order of the numerator polynomial and where dendeg is the order of the denominator polynomial. These two arguments are optional. If you don't specify them convert will just choose the appropriate orders for the number of terms in your series. Here are a few examples with illustrative plots using taylor to generate the series.

Approximate the sine function. Compare the exact result (red) with the rational aproximation (blue).

> sinap:=convert(taylor(sin(x),x=0,10),ratpoly);

> plot([sin(x),sinap],x=0..6,color=[red,navy]);

Approximate the tangent function. In this example notice how the ability of rational functions to be singular beautifully matches the singularity of tan(x) at . The rational expression assigned below to tanap is an excellent approximation to . (Exact in red, approximation in blue.)

> tanap:=convert(taylor(tan(x),x=0,10),ratpoly);

> plot([tan(x),tanap],x=0..5,-10..10,color=[red,navy]);

I know this may seem a little strange--why would I want to build a rational function approximation to a function I already have? Well, suppose you had to have Bessel functions in a C or Java program. You could use Maple to get a simple approximation to it using ratpoly , then use this approximation in your other program. (Be careful that you don't use the approximation outside of the domain where it works.)

> J0ap:=convert(taylor(BesselJ(0,x),x=0,10),ratpoly);

> plot([BesselJ(0,x),J0ap],x=0..5,color=[red,navy]);

Or even better, plot the difference between the exact function and the approximation so you know over what range the approximation is sufficiently good.

> plot(BesselJ(0,x)-J0ap,x=0..5,-.001..0.001);

>

Strategy:

Because Maple will not do exactly what you want it to do you will get very frustrated. Perhaps the best way to avoid ulcers as you use it is not to expect too much. In complicated problems you will need to do a lot of pencil-and-paper work to guide the calculation, but you can use Maple to do the details. Maple does lots better on small chunks of algebra than it does on big ones, so feed it bite-size pieces, then find creative ways of combining the pieces on your paper. And when you find what looks like a promising approach, use Maple to implement it. Finally, when you have used Maple, your head, and your hand to get a final beautifully simplified answer, plot it against the original ugly form to make sure that they are really the same.

Sometimes you will need to do algebra on matrices, and you will discover that it is a little tricky. Here is a section on how to handle problems like this.

Algebra on matrices filled with expressions

Maple will let you do all of the things discussed in the section above to expressions contained inside matrices and arrays, but you have to work a little harder. Suppose we have the following matrix of expressions

> restart;with(LinearAlgebra):

> F:=Matrix([[x^2+x*y^2,1/x+y/x],[(x-y)*(x+y),x^2+2*x*y+y^2]]);

and suppose we want to change its form using the algebra commands discussed in the previous section. We might try to substitute q for x in each term of the matrix using subs

> subs(x=q,F);

which is exactly what we wanted. (Note however: if you had used with(linalg) instead of with(LinearAlgebra) this wouldn't have worked. In general, use with(LinearAlgebra) whenever you can.)

Now let's try expanding each term using expand

> expand(F);

Why does subs work and not expand ? I don't know, but for algebra commands like expand that fail like this, the map command will do the job. This would be a good time to look up map , but here is basically how it works. The command map(command,symbol) applies the first argument ( command ) to each part of the expression(s) in the second argument (each term in the matrix in our case).

> map(expand,F);

> map(factor,evalm(F));

etc.. With the help of map you should be able to do about anything you want to matrices and vectors containing expressions.

>

Go to top of section

Well, this is a very brief introduction to a very large subject. You will only get good at algebra in Maple by doing your homework with it in spite of the pain and frustration. Here is an example of how it's done, followed by some problems for you to practice on.

Worked example:

Here's an example from wave physics. Carefully step through it watching what happens, then proceed to the problems below. The dispersion relation for electromagnetic waves in the ionosphere is

.

I want to find both the phase velocity and the group velocity as functions of . So first we solve for :

> restart:Eq:=omega^2=omega[p]^2 + k^2*c^2;

> solve(Eq,omega);

Hmm, that didn't work well at all. I'm going to show you what went wrong now, and it's very important. I was seduced by Maple's beautiful typesetting capability to write the plasma frequency with a subscript, the way it would be done in a textbook. Do not ever use this for anything but comments in your worksheet. When Maple got to omega[p] in the solve it got very confused because it has been asked to find the pth element of the array omega, but omega isn't an array and p isn't an integer. When I code I always use non-subscripted, but suggestive, variables like wp:

> restart:Eq:=omega^2=wp^2 + k^2*c^2;

> s:=solve(Eq,omega);

I just want the first one and will assign to the name w

> w:=s[1];

Ok, getting the phase velocity is now not too hard:

> vp:=w/k;

But I want it as a function of , so now I need . No problem

> s:=solve(Eq,k);

Choose the positive one and assign it to k

> k:=s[1];

Now look at the expression for vp with this form for k

> vp;

Maple can't simplify the square root in the top unless it knows something about , in this case that it is positive

> assume(omega>0);vp:=simplify(vp);

Now it's time to find the group velocity as a function of .

> vg:=diff(w,k);

Error, wrong number (or type) of parameters in function diff

Now what went wrong? Well, we did something bad that you will be tempted to do often, but which you should not do. At the beginning we had a function , but now we have solved for and assigned a formula to it, so is now ruined:

> w;

What we should have done is assigned the solution for to another variable, say kk. So let's do everything all over again

> restart:Eq:=omega^2=wp^2 + k^2*c^2;

> s:=solve(Eq,omega);

I just want the first one and will assign to the name w

> w:=s[1];

Ok, getting the phase velocity is not too hard:

> vp:=w/k;

But I want it as a function of ; no problem, just solve for

> s:=solve(Eq,k);

Choose the positive one and assign it to kk

kk:=s[1];

Now look at the expression for vp with this form for k by substituting kk for k

> vp:=subs(k=kk,vp);

Maple can't simplify the square root in the top unless it knows something about , in this case that it is positive

>

> assume(omega>0);vp:=simplify(vp);

Now find the group velocity

> vg:=diff(w,k);

Substitute for k again

> vg:=subs(k=kk,vg);

> vg:=simplify(vg);

>

and we are done.

Problem 9.1

(a) Express the function as a polynomial in

(b) Express the function as a sum of terms containing and .

(c) Coax Maple into giving you the half angle identities for , , and . Hint: a single simple command (one of the ones discussed above) will do this problem for you if you give it the squares of sine and cosine to work with. But it isn't obvious (to me anyway) which one to try. So just experiment until you find the magic one.

----------------------------------------------------------------------------

Problem 9.2

The formula for the effective resistance of a certain combination of resistors is given by

Beat on this with Maple until you can get a semi-nice formula for with in the numerator and a sum of terms in the denominator, one of which is a fraction with denominator and with a factor of in the numerator, times another factor, and the other two added terms in the denominator being . This is an exercise in the creative use of copy and paste.

----------------------------------------------------------------------------

Problem 9.3

In your first physics course you studied collisions between two particles. Let's do a simple case in which the two particles are constrained to move along the x-axis. The two particles have masses and and initial velocities and with the fast particle behind the slow one so they can collide. Find the final velocities and .

(a) Elastic collisons:

The conditions that allow this problem to be solved are conservation of momentum and energy:

and , where the has been cancelled out of the kinetic energy equation. Coax Maple into solving these two equations for the final velocities. You will get two solutions; interpret each set physically. Then specialize to the cases (i) , (ii) much larger than and so that the light particle comes up from behind and collides with the massive particle, and (iii) much larger than and so that the massive particle comes up from behind and whacks the light particle. Does each case make physical sense? Talk about it with someone until you understand each case intuitively.

----------------------------------------------------------------------------

Problem 9.4

In Freshman physics you also studied projectile motion with gravity obtaining the following equations describing the horizontal and vertical motion:

and . Suppose that and are given as well as the initial speed and the angle that the initial velocity makes with the horizontal plane. Also suppose that the elevation of the region in which the projectile will fall is known. Find a formula for the horizontal range and for the "hang-time" of the projectile (the time it spends in the air). Do not be satisfied with some mess involving RootOf--get workable formulas. And if you get two solutions, figure out which one is the physically relevant one. When you get final expressions for the horizontal range and the hang-time, choose , , m/ , and m/s and make plots of hang time and range as functions of if the elevation of the target region is (shooting down into a valley). Use in the range (straight down to straight up). Carefully examine the graphs to be sure they make physical sense. In particular, why doesn't 45 degrees give the best range in this case?

----------------------------------------------------------------------------

Problem 9.5

The formula relating the focal length , the object distance , and the image distance for a focusing lens or mirror is . Given and , this equation can be solved for the distance to the image. Assume that m (about right for the curved mirror in the ESC lobby) and make a plot of for ranging from 15 m to 0.5m, about the range that is possible in the lobby of the ESC. When becomes negative it means that the image is virtual (meaning that it's as if there were an image behind the mirror), and because your eye has a lens in it you can see virtual images, and they appear to be at the place given by this formula. Look at your graph, then go down to the lobby and see if what you see is what you got. A good movable object to use to define is yourself. In particular, think about what you should see if you let your face be the object and place it so that the object and the image are at the same place.

What you see in the lobby will not correspond in an obvious way with your graph because your eye has a lens in it too. This means, among other things, that you can see images at infinity! If you have a partner to argue with, this lobby exercise could be both fun and instructive.

----------------------------------------------------------------------------

Problem 9.6

Here is a thermodynamics problem that is a pain to do by hand, but Maple takes it right in stride. Suppose you have a quantity of water of mass with specific heat J/kg/K and a quantity of ice of mass and specific heat J/kg/K. The specific heat of transformation is J/kg. If the water starts out at temperature and the ice starts out at temperature , and if the freezing point is at temperature (we in are degrees celsius) do the following.

(a) Assuming that all of the ice melts (i) find a formula for the final temperature of the combined mass of water and (ii) also find a formula for the mass of ice such that the final water temperature is To (the final mixture has just barely melted). In case you have forgotten, or never studied, heat, here's some help. The heat lost by the water is given by the formula

and the heat gained by the ice that melts is given by the formula

Just use conservation of energy and the rest is not too hard. To get definite numbers out, use kg. (Notice that part (ii) is quite practical: given water at 20 C and ice at -10 C, how much ice does it take to bring 1 kg of water to be freezing cold?)

(b) Assuming that all of the water freezes completely and then cools down below (i) find a formula for the final temperature of the combined mass of ice and (ii) also find a formula for the mass of ice such that the final ice temperature is To (the final mixture is just barely frozen). Similar formulas to those used in part (a) also apply here, except this time it is the ice formula that is simple (it just warms up as ice from to , while the water cools to , then freezes into ice, then cools down further as ice to . Another practical question: How much ice would you need at -10 C to just barely freeze 1 kg of water at 20 C? The answer is impressive and shows that a lot of heat has to be removed to freeze water.

----------------------------------------------------------------------------

Problem 9.7

Light is incident from transparent region a into transparent region b, which in turn transmits light into region c, as shown below. The indices of refraction in each region are denoted by , , and .

> with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

> l2:=listplot([[0,-2],[0,2]]):

> l3:=listplot([[1,-2],[1,2]]):

> t1:=textplot([-.5,.5,"a"]):

> t2:=textplot([.5,.5,"b"]):

> t3:=textplot([1.5,.5,"c"]):

> l4:=plot(.5*sin(20*(x+1))-.5,x=-1..-.2):

> a1:=arrow([-.7,-1.3],[-.4,-1.3],.05,.2,.2):

> plots[display]([t1,l2,t2,l3,t3,l4,a1],axes=none,thickness=3,view=[-2..2,-2..3],font=[HELVETICA,14]);

>

Maxwell's equations tell us that for normal incidence, which we assume we have here, the ratios between the incident wave in region 1 and the reflected wave in region and the transmitted wave in region 2 are given by and . For convenience we will call these ratios and , respectively. Your job is to follow all of the successive transmissions and reflections in the figure above and arrive at the final ratio of the light reflected back into region a as a function of the three indices of refraction, the thickness of region b, and the wavenumber of the light in region b. To do this follow the steps outlined below.

(a) Note that the first reflection and transmission ratios are simply and . (Don't do anything except make sure you understand this statement because there are a lot more just like it coming up.)

(b) Now think about what happens to the light transmitted into region b. It first reflects from region c, is phase shifted by the complex phase factor due to the round trip through b from a to c and back, then is transmitted back into a. Show that the factor for this whole process, including the first transmission, is . (This is not too difficult, but be patient--hard stuff is coming up.)

(c) Now show that each successive multiple reflection ratio is given by where the first multiple reflection in part (b) corresponds to .

(d) Now combine all of this to get the total reflection ratio formula

. (Here's where the Maple actually starts--copy this into your exercise worksheet and and go to work.)

(e) Use the formulas for the reflection and transmission coefficients and Maple's result for the sum to build a formula for . Once you have it, manipulate it into a form where you can easily obtain separate equations for its real and imaginary parts. The reflected wave will be exactly zero if both real and imaginary parts vanish. You will save yourself some trouble if you define a variable and work with it. It will also simplify the look of some of the equations if you let . The formulas will become horrible, and then at the end they become simpler and you will discover that so that , and once you know this it becomes not too hard to get , the same result you get if you ignore the problem of multiple reflections and transmissions. This is a hard problem--it took about all the Maple skill I have, and I learned some new things along the way.

Go to top of chapter

>

>