> diff(tan(x),x);

> diff(tan(x),x$2);

> d:=Diff(tan(x),x$3);

> d:=value(d);

> d:=simplify(d);

Now let's see how to do functions.

> f:=x->tan(x)/x;

> fp:=D(f);

Note that assigning the result of D(f) to fp produces the function fp(x). There is more than one way to take higher derivatives, but here is the most general one.

> fpp:=D[1$2](f);

>

In the square brackets the 1 means differentiate with respect to the first (in this case the only) variable in the argument list and $2 means do it twice just as in diff .

Well, that's about it. Here are some exercises to practice on.

Problem 3.2

Take the indicated derivatives of the functions listed below. Use expression form on most of them but use function form on (a) and (d). When you get a mess try using simplify to clean it up. You will discover that simplify doesn't work on functions; to make the result look nice use the mouse to copy the mess you want to simplify to the clipboard, assign it to a new variable, cut off the extraneous junk, execute the command, then rebuild the derivative function using cut and paste again. This combination of Maple and editing is a good way to do mistake-free algebra.

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

(f) (g) (this is the complete elliptic integral of the first kind, Maple's EllipticK .)

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

Problem 3.3

Here's a find-the-maxima-and minima problem like ones you used to do in high school. Consider the function . (Note: I am using the word function in the mathematical sense, not in the Maple sense. This problem is easier if you just use a Maple expression to define the function above.) (a) First plot this function on the interval [0,10]. (b) Look at the graph and find, approximately, the values of where the function has maxima and minima. Then take the derivative of the function and use the fsolve command to find these values of precisely. If your expression for the derivative is called , and if you wanted to find a zero of it near 1.1, you would do this: fsolve(f,x=1.1) ;

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

Problem 3.4

. This definition is too cumbersome for most computer languages, but Maple handles it with ease because it thinks symbolically instead of numerically. Here is a function that evaluates it

> with(orthopoly);

> Pnm:=(n,m,x)-> (-1)^m*(1-x^2)^(m/2)*diff(P(n,x),x$m);

Before we do anything fancy with it we ought to test it, so let's put in numbers for , , and .

> Pnm(3,1,.5);

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

Well, we're in trouble again. The problem is with what P(n,x) returns. As we saw in the section in Chapter 2 on this function, it doesn't return a number; it returns a polynomial. When we put in 0.5 for x it goes into our Pnm function, defined above, in place of x, and then diff tries to take a derivative with respect to 0.5, which makes no sense. Watch what happens when we put in a variable in place of x instead of a number

> Pnm(3,1,t);

So if you want a numerical result you have to do something like this

> a:=Pnm(3,1,t);t:=0.5;a;

This is annoying, but on the other hand, just think about it; why do you need a number in Maple anyway? You want to plot the function, differentiate it, integrate it, use it in differential equations, etc.. What could be better than having an explicit expression for the thing? So Maple doesn't think this is a problem; it's a feature! And this feature is shared by all of the orthogonal functions you get when you use with(orthopoly).

There is one other annoying thing about the function Pnm. Watch what happens when we try to evaluate it with m=0

> Pnm(5,0,x);

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

>

With m=0 it is supposed to give us Pn(x) back, but it doesn't. The reason it doesn't is that we just asked it to take the 0th derivative of a function, which Maple's diff command can't handle. Later when we study procedures we will come back to this problem and fix it so it works with m=0.

Ok, I have shown you how this can be done. Now I want you to make a plot of the 5 associated Legendre functions that go with P(5,x), i.e., n=5 and m=1,2,3,4,5. Plot them from x=-1..1. Put all five plots on one set of axes with different colors so you can see what happens as m ranges from 1 up to n=5. After you look at the picture you will probably want to rescale the functions so they all look about the same size. In the upcoming section on integration we will redo this plot and learn a natural way to make the functions nearer to the same size.

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

Problem 3.5

Here is an electricity problem borrowed from the integration section that follows this one. The potential function for a charged hemisphere of radius R with surface charge density as a function of z , where z goes up the symmetry axis of the hemisphere, is given by the expression

> V:=-1/2*sigma*R*(-sqrt(R^2+z^2)+sqrt((z-R)^2))/(z*e0);

where e0 is shorthand for the electrical constant .

>

The electric field component can be obtained from the potential V by differentiating: . Use Maple to take this derivative and obtain an (ugly) expression for . Simplify it. You will see an unfamiliar function called csgn ; look this function up with ?csgn and make sure you understand what it does. Then set , , and and plot both V and from to . It is a theorem in electromagnetism that when you cross surface charge density, the electric field jumps up by . (You will notice that in the definition of V above I used e0 in place of . This was intentional. Avoid subscripted variables whenever possible because subscripts in Maple reference array elements.) Verify on your plot that the correct jump is obtained. On your graph negative z is below the circular rim of the hemisphere; positive z from 0 to R is inside the hemisphere; and positive z from R to infinity is above the dome. Interpret your graph physically and convince yourself that it makes sense.

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

Problem 3.6

> restart;

> eq:=x^2 + y(x)^3 = 3;

> deq:=diff(eq,x);

> dydx:=solve(deq,diff(y(x),x));

>

If you don't want to type y(x) all of the time, you can use Maple's alias command to tell it to turn y into y(x) (only in Maple's internal processing) whenever it is encountered.

> restart;

> alias(y=y(x));

> eq:=x^2 + y^3 = 3;

> deq:=diff(eq,x);

> dydx:=solve(deq,diff(y,x));

>

Here's an example where this shows up in physics. The dispersion relation for electromagnetic waves in a plasma is , where is a frequency called the plasma frequency. The phase velocity of these waves is given by while the group velocity is given by . First use Maple to find formulas for both the phase and group velocities in terms of , , and by solving for and differentiating. Then use implicit differentiation to get the group velocity in terms of , , and .

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

> restart;f:=cos(x*y)/y;

> diff(f,x);diff(f,y);diff(f,x,y);

And here is the same thing using Maple function notation

> restart;f:=(x,y)->cos(x*y)/y;

> D[1](f);D[2](f);D[1,2](f);

>

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

Problem 3.7

Find the two first derivatives and all three second derivatives (double x, double y, and mixed xy) of the following function

, where K is the complete elliptic integral EllipticK . Use expression notation with diff. Try using expand and simplify to clean up the messes that result.

> int(sin(x),x);

or

> f:=sin(x)*x;int(f,x);

Note: do not use f(x) as the argument if f is an expression.

With functions the integration command works like this

> g:=(x,y)->sin(x*y)*x;

> int(g(x,y),x);

T

There is also an inert form of int, namely Int , which is used to display integrals. You would use this form for documenting a worksheet. Try this:

> s1:=Int(exp(x),x);

But be warned: Int only displays, it doesn't do math. "But if it doesn't do anything, why would I want to use it?", you ask? Because it will help you see if you have entered the integral properly, that's why. This makes Int a very valuable debugging tool. And after you have seen that it looks right in displayed form, use value(s1) to get the answer. So the right way to do the simple integral above and get an answer is this

> s1:=Int(exp(x),x);

> s1:=value(s1);

I suggest that you always do integrals this way, combining Int and value . It is one of those habits that will decrease the number of hours you will spend looking for stupid mistakes.

> s2:=Int(tan(x),x=0..1);

> s2:=value(s2);

And if you need a number there is always

> evalf(s2);

Oh, and if you just want to get a numerical answer without going through evalf , just give int floating point limits and you will get it right away.

> s2:=Int(tan(x),x=0..1.);

> value(s2);

>

You should also know that Maple can do integrals where the limits are infinite, but you may need to give it some guidance by using the assume command. Well, that's about all you have to know if Maple can do the integral. Try ?int to see some of the extra integration options Maple has available. Now let's do some exercises.

Problem 3.8

Do the following integrals with Maple using expression notation for (a)-(d) and function notation for (e)-(g). Obtain numerical values for (e) and (f). You will have trouble with (g), and even when you get it to work the answer will look ugly. Try simplifying it with simplify .

(a) (b) (c) (d)

(e) (f) (try using both 1/2 and 1./2. as the upper limit) (g) (Don't know how to handle ? Use online help.)

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

Problem 3.9

I will get you started. The vector distance to the observation point on the -axis is r =z k , where k is the unit vector in the -direction. The vector distance to a little bit of charge is r '=R r' where r ' is a unit vector pointing from the origin to the little bit of charge. Using the law of cosines we can then write

where the integral sign here means integrate over both and . The integral just gives and we are left with

Use Maple's Int command to build a function V(z) that we can use later. The syntax for function building is this (setting V(z) to the dying exponential function).

> V:=z->exp(-z);

Note: recall that you had better not use the symbol V(z); Maple interprets this literally, i.e., V(z) is the name of the function you are defining, requiring you to use the awkward syntax V(z)(z) when using it. As defined above you type V(z) when you use it. You can do the same thing with int , like this V:=z-> int(....) . In place of the symbol use e0, and define the constants so Maple can plot numbers. Here, I'll get you started; you fill in what goes in (...)

> restart:e0:=8.854e-12;sigma:=1e-10;R:=0.5;V:=z->int(sigma*R^2*sin(theta)/sqrt(R^2+z^2-2*R*z*cos(theta)),theta = 0 .. Pi/2)/(2*e0);

>

You will notice that Maple doesn't give you a formula for the answer this time; it just displays the integral. I don't know why. But now that it is a function you can plot it as a function of z. Do so, from z=-5 to z=5.

Now calculate the total charge on the hemisphere and overlay plots of the true potential and the corresponding point charge potential. The two plots should agree for |z|>>R. You can get a better fit by putting the point charge somewhere other than at the origin; can you graphically find about the best value of z at which to put the point charge? Recall that the formula for the potential of a point charge is .

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

Problem 3.10

where is the small-oscillation frequency of the pendulum. If we multiply this equation by and integrate once with respect to time, we get the conservation of energy theorem for this problem:

where is the initial angular position of the pendulum before we release it from rest. We would like to know how long it takes the pendulum to get from its initial angle down to because this time is , one-quarter of the period of a pendulum with initial amplitude . (a) Show that the energy equation above can be rewritten in the form

.

This means that the quarter period of the motion is given by the integral

.

(b) Have Maple do this integral in an attempt to get a formula for the period T of a pendulum, using s and s0 instead of and . Once you try it you will probably be unimpressed with the beauty of the result. Part of the problem is that Maple doesn't know what you are going to use for , so it gives you something incredibly general (as far as it knows, you might even want to be complex). To tie Maple down a little, use the assume command, like this

> assume(theta0,real,0<theta0,theta0<Pi);

>

(Note: when you see these variables again they will have little ~ signs attached to them. This is Maple's way of showing that they have been restricted in meaning by an assume command.) Now you can try integrating again; it will look a little better, but it is still ugly.

You can often help Maple out by using integration tricks, like u-substitutions. Try this integral again, but use the identity to change the functions inside the square root. You and Maple are a team; don't be afraid to use what you have learned in your math classes.

(c) Plot the period T as function of initial angle from 0 to . Does your plot make physical sense? Check the plot for both small and near .

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

Problem 3.11

You might first try this

> restart;

> s1:=Int(BesselJ(0,x),x=0..1);

> value(s1);

and be annoyed because you haven't ever heard of Struve functions. What you really want is just a number, which you can get via evalf , or by the following subtle change in the command (watch the decimal points).

> s1:=Int(BesselJ(0,x),x=0..1.);

> value(s1);

>

That little decimal point on the 1. makes all the difference. Ok, here are some exercises for you to practice on.

Problem 3.12

Try doing the following Bessel function integrals both at indefinite integrals and as definite integrals on the interval [1.0..2.0].

(a) (b) (c) (d)

(Note that you get these functions with BesselJ(0,x) , BesselK(0,x) , and BesselI(1,x) .

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

Problem 3.13

(a) (b) (Note: to get these functions you have to use the command with(orthopoly) first, then P(3,x) , etc..)

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

Problem 3.14

If I were to ask you to check whether two vectors were orthogonal, you would know just what to do. You would take their dot product, and if it turned out to be zero, then you would know that they are orthogonal. And how would you take the dot product? You would multiply the like components of the two vectors together and add them up, e.g.,

[1,2,3] dot [4,5,6] = 4 + 10 + 18 = 32

Well, an integral like this is really just the same thing. At each value of the two function values are multiplied together and the integral adds them all up. So it makes sense to say that two functions are orthogonal. It turns out that many functions have nice orthogonality properties, and in this exercise you will use Maple to see what they are like.

(a) Show, using Maple, that the functions and are orthogonal to each other when integrated over the interval by finding the values of the following integrals both with and unequal and with .

(Hint: you will get a mess unless you use assume to tell Maple that and are integers. Apparently Maple assumes more than that they are integers; it also assumes that they are different integers, so you will have to do the case by using the command m:=n .)

(b) Show, using Maple (this problem is an extension of Problem 3.13 above), that the Legendre functions are orthogonal to each other when two functions with different -values are integrated over the interval [-1.0..1.0]. Also find the value of . You will need to use the with(orthopoly) command to tell Maple you are going to use these functions. And when you try to find the value of the integral of , you will probably need to try several different values of and see if you can find the general formula empirically. This will be easier if you don't force floating point answers by using decimal points on the limits, i.e., use 1 and -1 instead of -1. and 1..

(c) Show, using Maple, that the Bessel functions are orthogonal to each other with respect to the integral , where the numbers and are the zeros of . These zeros are available from the Maple function BesselJZeros(n,m) , where is the order of the Bessel function, i.e., for , for , etc., and where counts the zeros ( and are -values), starting with the first one. (For example BesselJZeros(0.,1) gives 2.404825558; notice the decimal point--it is important.) Then see if you can find a formula for this integral when . (Hint: when you set you will just get meaningless numbers because BesselJZeros gives you floating point numbers back. You can do better than this. Replace BesselJZeros(0.,n) with BesselJZeros(0,n) (no decimal point) and let Maple try to do the integral symbolically. You will get a simple answer.

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

Problem 3.15

Let's go back to the associated Legendre functions and plot them so that they all have about the same size. We could scale them by their maximum values, but then we would have to find the zeros of large polynomials, which is hard. It is actually easier to "calculate their magnitudes" in the same function-dot-product sense used in Problem 3.14. To find the magnitude of a vector we square the components, add the squares up, then take the square root. The corresponding operation for functions is: square the function, integrate it over the relevant interval, then take the square root. Here is the associated Legendre function definition again from Problem 3.4

> Pnm:=(n,m,x)-> (-1)^m*sqrt(1-x^2)^m*diff(P(n,x),x$m);

>

The magnitude of this function is obtained from . Redo the plots in Problem 3.4, but divide each of the plotted functions by this magnitude and check to see that they all have about the same size on the new plot.

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

Problem 3.16

Finally, let's see what Maple does with the integrals of a few miscellaneous functions. Try to get Maple to do (a-c) both as indefinite integrals and as definite integrals between .5 and 1. Do (d) only as a definite integral.

(a)

(b)

(c) (Maple's EllipticK )

(d) (Note: try this integral with both integer 1 and floating point 1. as the upper limit.)



This should teach you that sometimes Maple can't do indefinite integrals, that it can give you a floating point number for a definite integral but you might have to insist on it by using evalf , and that sometimes it will give you the wrong answer! (Plot between 0 and 1 to figure out which of the two answers Maple gave you for (d) is wrong.) Maple doesn't make many mistakes, but occasionally it louses up. This means that you have to be very careful and not trust Maple, or any other computer program, too much. Indeed, a good rule to follow upon receiving an answer from any computer program is to quietly mutter under your breath, "Well, that's probably wrong."

> restart;

> s1:=Int(Int(erf(x*sin(y)),x=0..1),y=0..1);

> s1:=value(s1);

As usual, if Maple just gives it back to you as a displayed integral it means it can't do it, but if you want a number you can use evalf

> evalf(s1);

Or, if you just wanted a number in the first place, you can tell Maple to do the integral numerically using this command

> evalf(Int(Int(erf(x*sin(y)),x=0..1.),y=0..1.));

A variation on this command tells Maple to only compute the integral numerically accurate to 6 (or any other number) of significant figures.

> evalf(Int(Int(erf(x*sin(y)),x=0..1.,6),y=0..1.,6));

>

You might want to use this if Maple is taking forever to give you the answer--it will go faster if you ask for fewer significant figures.

> with(student);

To see what these routines do and how to use them use

> ?student

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

Problem 3.17

Get the numerical answer to this double integral .

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

Advanced topic: procedure for triple integrals

> restart;

> s1:=int(int(int(erf(x*sin(y*z)),x=0..1),y=0..1),z=0..1);

And when Maple gives it back unevaluated, you could use evalf again. But when you do you will notice something rather annoying--Maple takes forever to do it. The problem is that for triple integrals there are so many calculations to do that Maple itself is just too slow. You need a real number cruncher program to get the job done. And for integrals of dimension 4 or higher, even good number crunchers can take a long time. Go ahead and try evalf, but be prepared to stop the computation with the stop sign on the toolbar when you can't take it anymore.

> evalf(s1);

Even using Maple's numerical integration command with a limited number of significant figures takes forever.

> evalf(Int(Int(Int(erf(x*sin(y*z)),x=0..1.,4),y=0..1.,4),z=0..1.,4));

Fortunately, there is an alternative within Maple: write your own procedure to turn Maple into a "number cruncher" (insofar as this is possible). This requires that we use Maple's procedure capability which will be discussed in a later section. Don't worry yet about how it works; just know that it's here and that you can come copy it if you have to use Maple to do a triple integral.

Here is the numerical integration procedure together with the definition of the function f(x,y,z) corresponding to the example we just did so you can see how much faster it is than three nested int's. In the int3d command after the procedure try different numbers of intervals nx, ny, and nz to see if the answer converges to what Maple got when it did it to 4 significant figures. For example, you could first have nx=ny=nz=4, then nx=ny=nz=8, etc..

> f:=(x,y,z)-> erf(x*sin(y*z));

> int3d:=proc(x1,x2,y1,y2,z1,z2,nx,ny,nz,f) # 3-d integration procedure

>

> # x limits: x1,x2 y limits: y1,y2 z limits: z1,z2

> # number of subintervals (3 points per subinterval): nx, ny, nz

> # function f(x,y,z) passed in to be integrated: f

>

> local x,y,z,dx,dy,dz,wx,wy,wz,i,j,k,s; # define variables local to this procedure

> wx:=array(1..3*nx);wy:=array(1..3*ny);wz:=array(1..3*nz); # declare the weight arrays

>

> dx:=(x2-x1)/nx/3;dy:=(y2-y1)/ny/3;dz:=(z2-z1)/nz/3; # set the integration step sizes

>

> x:=array(1..3*nx);y:=array(1..3*ny);z:=array(1..3*nz); # declare the integration point arrays

>

> for i from 1 to nx do # load x integration points and weights

> wx[3*i-2]:=9/8.*dx;wx[3*i-1]:=3/4.*dx;wx[3*i]:=9/8.*dx;

> x[3*i-2]:=x1+(3*i-2.5)*dx;x[3*i-1]:=x1+(3*i-1.5)*dx;x[3*i]:=x1+(3*i-.5)*dx;

> end do;

>

> for j from 1 to ny do # load y integration points and weights

> wy[3*j-2]:=9/8.*dy;wy[3*j-1]:=3/4.*dy;wy[3*j]:=9/8.*dy;

> y[3*j-2]:=y1+(3*j-2.5)*dy;y[3*j-1]:=y1+(3*j-1.5)*dy;y[3*j]:=y1+(3*j-.5)*dy;

> end do;

>

> for k from 1 to nz do # load z integration points and weights

> wz[3*k-2]:=9/8.*dz;wz[3*k-1]:=3/4.*dz;wz[3*k]:=9/8.*dz;

> z[3*k-2]:=z1+(3*k-2.5)*dz;z[3*k-1]:=z1+(3*k-1.5)*dz;z[3*k]:=z1+(3*k-.5)*dz;

> end do;

>

>

> s:=0.; # set the final answer to zero

> for i from 1 to 3*nx do

> for j from 1 to 3*ny do

> for k from 1 to 3*nz do

> # loop over all the points, evaluate the function, multiply by the weights, and

> # add each contribution into the final answer using Maple's fast hardware floating

> # point eval command

>

> s:=s + evalhf(wx[i]*wy[j]*wz[k]*f(x[i],y[j],z[k])):

>

> end do:end do:end do:

>

>

>

> end; # end of int3d procedure

Test the int3d procedure

> int3d(0,1,0,1,0,1,4,4,4,f);

>

I recommend that you use this procedure whenever you need to do triple integrals in Maple because it is much faster. And to tell whether you have the accuracy you desire, just try increasing (nx,ny,nz) and watch how many digits stay the same. But if you are really serious about doing a time-consuming problem like this and Maple is just taking way too long, there is still hope. There are lots faster ways of doing numerical calculations like this and you ought to learn about them sometime. For instance I coded this triple integral in two other computing languages, MATLAB and FORTRAN. All three give the same answer for the numerical approximation to the integral as long as the same number of points is used. For example, with a 12,12,12 in the int3d command above Maple takes 48 seconds to do this integral on my workstation. MATLAB takes 3.7 seconds and FORTRAN takes 0.24 seconds, 200 times faster than Maple! Number crunching like this using Maple turns your Pentium into a clunky old 286.

Problem 3.18

Use the procedure int3d defined above to get the numerical value of the following triple integral in spherical coordinates

.

Make sure your answer is accurate to 4 significant figures. Even this level of accuracy will take some time to achieve. It would be worth your while to try to find a way to do one of these integrals analytically so that you only have to do a double integral instead of a triple one.

Taylor's theorem is: that if a function is defined and has well-defined derivatives at the point , then for near the function can be approximated by the series

+ . . .

Maple knows this theorem and can use it to generate all kinds of series with the command taylor . For instance, here is the one for the function discussed above

> taylor(sech(x),x=0,20);

>

The x=0 option tells Maple to take the expansion about 0 and the 20 tells it to include terms until the error term is of order .

Problem 3.19

Find the Taylor expansions of the following functions about and keeping terms out to .

(a) (b) (c) (d) (e) (f)

Your answers will look pretty impressive, but you will find that in most physics problems where expansions are important, we only keep the first two terms, so you can often get what you want without all of the clutter by using taylor(f(x),x=0,3) , which just goes through .

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

Problem 3.20

Maple can also do Taylor expansions of special functions. Find the expansions of (a) , (b) , (c) (complete elliptic integral EllipticK ), and (d) . Do (a)-(c) about and (d) about .Do each one through . The answer to (d) is pretty awful--see if you can get something more reasonable by using evalf , or by telling Maple to expand about x=1.0 (decimal point again).

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

Problem 3.21

(a) (b) (c) .

Try different numbers in place of the 3 in the series command above and see if you can tell precisely what it means.

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

Problem 3.22

(a) (b) (c) (d) .

I have suggested that you only put a 2 in the last slot because usually all we care about in an asymptotic expansion is the leading term because it answers the most important questions: (1) does the function blow up or die? and (2) in what way does it do it?

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

Problem 3.23

Here's an example of how we use this idea to solve physics problems. The equation of