Unit 4: Vector Calculus

Chapter 22: Non-Cartesian Coordinates

Sections22.2: vector operators in polar coordinates

Copyright

Copyright * 2001 by Addison Wesley Longman, Inc.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any form or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher. Printed in the United States of America.

Initializations

> restart;

> with(linalg):
with(plots):
with(plottools):
read(`C:/Program Files/Maple 6.01/pvac.txt`):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>

Gradient in Polar Coordinates

In cartesian coordinates, the gradient of the scalar is the vector

> grad(u(x,y),[x,y]);

>

which is interpreted as the vector

i + j

where i and j are unit basis vectors for the cartesian system. What does it mean to compute the gradient of the scalar in polar coordinates? What does Maple give?

> grad(u(r,theta),[r,theta],coords=polar);

>

Next, we explore how this vector was obtained, and what it actually means.

>

Gradient in Polar Coordinates - Derivation

Recall the equations defining polar coordinates

> X := r*cos(theta);
Y := r*sin(theta);

>

and the equations defining the inverse transformation

> R := sqrt(x^2+y^2);
T := arctan(y/x);

>

Notice the care with which name clashes were avoided.

Write the radius vector R = x i + y j in Maple as Rc , where and , obtaining

> Rc := vector([X,Y]);

>

If is held constant and r allowed to vary, R = R ( r ), so an r -coordinate curve is traced in the xy -plane. For example, set = 1,

> R1 := subs(theta=1,op(Rc));

>

and plot the resulting r -coordinate curve which turns out to be a radial line emanating from the origin.

> plot([R1[1], R1[2], r = 0..1], color=black, labels=[x,y]);

>

Alternatively, set r = 1 and let vary, so the radius vector defines the -coordinate curve

> R2 := subs(r = 1, op(Rc));

>

which, when plotted, is, a circle centered at the origin.

> plot([R2[1], R2[2], theta = 0..2*Pi], color=black, labels=[x,y], scaling=constrained);

>

A fundamental idea of vector calculus, now long familiar, is that differentiation of a curve with respect to its parameter yields a vector tangent to that curve. Apply that principle to the radius vector R( r , ), first differentiating with respect to r , and then with respect to . In the first case, obtain a vector tangent to an r -coordinate curve (a radial line) and in the second, a vector tangent to a -coordinate curve, a circle about the origin.

Figure 22.7, created below, shows a sketch of the vectors R and R drawn at .

> p1 := circle([0,0],2):
p2 := plot([[0,0],[4*cos(Pi/4),4*sin(Pi/4)]],color=black):
p3 := arrow([sqrt(2),sqrt(2)],vector([1,1]),.1,.3,.2,color=black):
p4 := arrow([sqrt(2),sqrt(2)],vector([-2,2]),.1,.3,.2,color=black):
p5 := textplot([.5,.2,`q`],font=[SYMBOL,10]):
p6 := textplot({[2.25,1.78,`R `],[-1,2.9,`R `]},font=[TIMES,BOLD,12]):
p7 := textplot([2.45,1.72,`r `],font=[TIMES,ITALIC,10]):
p8 := textplot([-.8,2.8,`q `],font=[SYMBOL,10]):
p9 := textplot([.2,3.6,`y`], font=[TIMES,ITALIC,12]):
display([p||(1..9)],scaling=constrained, xtickmarks=5, ytickmarks=[-2,-1,0,1,2,4], view=[-2..3,-2..4], labels=[x,``], labelfont=[TIMES,ITALIC,12]);

>

The vector tangent to the r -coordinate curve is automatically a unit vector. The vector tangent to the -coordinate curve needs to be normalized since it does not automatically have unit length. Typical notation for these unit vectors in polar coordinates is

e = R = i + j

and

e R i + j

However, to avoid violating Lopez's Large Law (Never use as a name on the left, a variable in use on the right), name these unit tangent vectors and in Maple.

First, differentiate with respect to r .

> E[r] := map(diff,Rc,r);

>

Second, differentiate with respect to

> Q := map(diff,Rc,theta);

>

and normalize, obtaining

> E[theta] := map(simplify,normalize(Q),symbolic);

>

Solving for i and j in terms of e and e , we have

i = e e

and

j = e e

which we obtain in Maple as follows.

Introduce the i and j basis vectors of the cartesian frame by the following device.

> q1 := e[r] = dotprod(E[r],[i,j], orthogonal);
q2 := e[theta] = dotprod(E[theta],[i,j], orthogonal);

>

These are now two equations in the two unknowns i and j which can be solved for i and j in terms of the polar coordinate basis vectors and .

> q3 := solve({q1,q2}, {i,j});

>

To emphasize this result, induce Maple to write i and j in terms of and .

> i = subs(q3, i);
j = subs(q3, j);

>

The point here is to discover the proper expressions for the gradient and laplacian of a scalar-valued function, and the divergence of a vector-valued function, when each is given in polar coordinates. Since the expression for the gradient in cartesian coordinates is already known, use that knowledge to deduce the result for polar coordinates.

Start with the expression for the gradient in cartesian coordinates.

> q4 := grad(u(x,y),[x,y]);

>

This gradient vector, given explicitly in terms of the fixed cartesian basis vectors i and j , is

> q5 := q4[1]*i + q4[2]*j;

>

The conversion process begins by replacing the vectors i and j with their equivalents in terms of and .

> q6 := subs(q3, q5);

>

This looks a bit better if the coefficients of the basis vectors and are collected via

> q7 := collect(q6, [e[r], e[theta]]);

>

Then, we have to transform the partial derivatives by the use of the chain rule.

The second step expresses as prior to using the chain rule for determining the equivalents of the partial derivatives and in polar coordinates. For insight, consider the function

> h := (x,y) -> x*y^2;

>

and change to polar coordinates, obtaining , which is then

> H := unapply(h(X,Y),r,theta);

>

The introduction of is essential, since it is a different function of and than was a function of and .

In the first variable, , simply multiplies the square of the second.

In , the first variable, , is cubed.

Thus, two different letters, and are more than appropriate, they are necessary.

Now, go back the other way by computing , obtained in Maple via

> simplify(H(R,T),symbolic);

>

This second identity, namely,

is the starting point for the partial differentiations about to take place.

Thus, start with

entered into Maple as

> u := U(R, T);

>

Differentiate with respect to x and y to obtain the partials

and

Since

the gradient in cartesian coordinates becomes

[ e e ] + [ e e ]

= e e

To implement these calculations in Maple, begin with the derivatives and in the form

> ux := diff(u,x);
uy := diff(u,y);

>

These results are nothing more than the chain rule, which in subscript notation would be written as

Use these expressions for and in q7 , namely, in

to obtain

> q8 := subs(diff(u(x,y),x)=ux, diff(u(x,y),y)=uy, q7);

>

This represents as a vector in the polar coordinate basis vectors and . However, all traces of the variables and are to be removed. So, to begin, replace and in the arguments of U with r and , respectively.

> q9 := subs(R=r,T=theta,q8);

>

There are still other 's and 's left. Replace them with and respectively,obtaining

> q10 := simplify(subs(x=X,y=Y,q9), assume=positive);

>

A conversion from D-notation to the partials notation yields

> q11 := convert(q10,diff);

>

Since this is actually a vector, split the fraction into two parts to isolate each component of the gradient vector now expressed completely in polar coordinates.

> q12 := expand(q11);

>

Expression q12 represents the polar coodinate equivalent to the cartesian gradient vector grad( u ) = i + j . Writing the result as the column vector

> vector([coeff(q12,e[r]),coeff(q12,e[theta])]);

>

allows a direct comparison with Maple's

> grad(U(r,theta),[r,theta], coords=polar);

>

Divergence in Polar Coordinates

In cartesian coordinates, the vector field

> F := vector([f(x,y),g(x,y)]);

>

has as its divergence the scalar field

> diverge(F,[x,y]);

>

The divergence of the vector

B = e + e

that is, of the vector

> U:='U':
B := vector([U(r,theta),V(r,theta)]);

>

is the scalar

that is,

> expand(diverge(B,[r,theta],coords=polar));

>

Notice the term without any derivative! There is something to discover here!

>

Divergence in Polar Coordinates - Derivation

To derive the expression for the divergence of the polar vector

B = e + e

transform B completely to cartesian coordinates, compute the divergence in cartesian coordinates, then transform that result back to polar coordinates. Therefore, in terms of i and j , we have

B = i + j

where

and

Since and are functions of and , define

and

so that the cartesian version of B is really

B = i + j

Then, the partial derivatives and , computed by the chain rule, are precisely what are needed for determining the divergence of B .

By the chain rule,

and

so

and

The sum

= div( B )

is the divergence in polar coordinates.

To implement these calculations in Maple, we begin by writing an arbitrary vector in polar coordinates, using the unit basis vectors e and e in polar coordinates.

> B := U*e[r] + V*e[theta];

>

The basis vectors and have already been given in terms of the cartesian basis vectors, i and j by

> q1;
q2;

>

Transform the given polar vector completely to cartesian coordinates, compute the divergence in cartesian coordinates, then transform that result back to polar coordinates. That will clarify what replaces div( B ) in polar coordinates.

This is the vector B in terms of i and j .

> Bij := subs(q1, q2, B);

>

Collecting the i and j components is helpful.

> Bij_collected := collect(Bij, [i,j]);

>

Next, extract the i and j components of this vector which is now in cartesian form.

> F := coeff(Bij_collected,i);
G := coeff(Bij_collected,j);

>

The polar vector B = + has become the cartesian vector i + j . Unfortunately, and are functions of and , so define

so that the cartesian version of B is really i + j . Then, the partial derivatives and are precisely what are needed for determining the divergence of B .

To express and , as well as sin( ) and cos( ), in terms of and , define

and

Efficiency is achieved by setting up these substitutions in a template.

> q13 := {U=U(R,T),V=V(R,T),theta=T};

>

Making these substitutions in the polar components and yields the cartesian components and .

> f := simplify(subs(q13, F), symbolic);
g := simplify(subs(q13, G), symbolic);

>

The partial derivatives and can now be taken, yielding

> fx := diff(f,x);
gy := diff(g,y);

>

The sum is div( B ) in cartesian coordinates, the system in which the divergence operator was originally defined. In that sum, begin the process of returning to polar coordinates.

> q14 := subs(R=r,T=theta, fx + gy);

>

There are yet more 's and 's to convert to polar coordinates.

> q15 := simplify(subs(x=X,y=Y,q14),symbolic);

>

Next, change notation for the derivatives, from the D-notation to the partials notation, obtaining

> q16 := convert(q15,diff);

>

Making separate fractions yields the final result.

> q17 := expand(q16);

>

Once again turning to Maple for a comparison, write the general vector as

> B2 := vector([U(r,theta), V(r,theta)]);

>

and apply the built-in diverge command, obtaining

> `div(B2)` := diverge(B2,[r,theta],coords=polar);

>

This is the result just derived from first principles. Again creating separate fractions,

> expand(`div(B2)`);

>

the results are seen to agree.

> q17;

>

Laplacian in Polar Coordinates

In rectangular coordinates, the laplacian of is

> unassign('f','F');
laplacian(f(x,y),[x,y]);

>

In polar coordinates, the laplacian of is

that is,

> expand(laplacian(F(r,theta),[r,theta],coords=polar));

>

Example 22.2

For example, the laplacian of the function

> f := x + y^3;

>

is

> qc := laplacian(f, [x,y]);

>

Converting the function to polar coordinates gives , that is,

> F := subs(x = r*cos(theta), y = r*sin(theta), f);

>

In polar coordinates, the laplacian of is

> qp := laplacian(F, [r,theta], coords = polar);

>

an expression which actually simplifies to

> combine(qp,trig);

>

If , the laplacian in rectangular coordinates, were converted to polar coordinates, would result, in perfect agreement with the laplacian computed directly in polar coordinates. In fact, this comparison is made in Maple via

> simplify(6*r*sin(theta) - qp);

>

Maple Derivation

Instead of the very general , let x and y be specifically given in polar coordinates. Then,

=

> unassign('F');
f := F(sqrt(x^2+y^2), arctan(y/x));

>

The cartesian derivatives can computed as

> fx := diff(f,x);

> fxx := diff(f,x,x);

> fy := diff(f,y);

> fyy := diff(f,y,y);

>

Hence, the laplacian in cartesian coordinates has the following polar equivalent.

> q := fxx + fyy;

>

This unwieldy expression can be simplified to the much smaller

> q1 := simplify(q);

>

Further simplification comes notationally, by expressing the arguments of F as and , rather than as and , respectively. This step is done by making the following replacements.

> q2 := subs(R=r,T=theta,q1);

>

Still more simplification is achieved by changing the remaining apearances of and to their counterparts in polar coordinates, obtaining

> q3 := simplify(subs(x=X,y=Y,q2),symbolic);

>

Finally, convert from the D-notation for derivatives to the partials notation, getting

> expand(convert(q3,diff));

>

Comparing this result to Maple's built-in version

> expand(laplacian(F(r,theta),[r,theta],coords=polar));

>

shows the results to be identical.

>

Traditional Chain Rule Derivation

The Maple derivation above might seem intimidating and complex, but it is actually easier than the traditional derivation based on a direct use of the chain rule. The following is a Maple-implemented version of the traditional computation.

The starting point for the differentiations is the statement

An application of the chain rule gives

and

Second derivatives pose a greater challenge unless the dependence of the derivatives and on and is properly appreciated. Again using the chain rule, the second derivatives become

=

and

=

We can represent these derivatives in Maple as

> fxx := (F[rr]*r[x]+F[r*theta]*theta[x])*r[x] + F[r]*r[xx] + (F[r*theta]*r[x] + F[theta,theta]*theta[x])*theta[x] + F[theta]*theta[xx];
fyy := (F[rr]*r[y]+F[r*theta]*theta[y])*r[y] + F[r]*r[yy] + (F[r*theta]*r[y] + F[theta,theta]*theta[y])*theta[y] + F[theta]*theta[yy];

>

where we have utilized the equality of the derivatives and .

Then, the derivatives must all be computed and expressed in terms of the polar coordinate variables and . Thus, the derivatives are computed and converted via

> rx := simplify(subs(x=X,y=Y,diff(R,x)), symbolic);
rxx := simplify(subs(x=X,y=Y,diff(R,x,x)),symbolic);
ry := simplify(subs(x=X,y=Y,diff(R,y)), symbolic);
ryy := simplify(subs(x=X,y=Y,diff(R,y,y)),symbolic);
theta_x := simplify(subs(x=X,y=Y,diff(T,x)), symbolic);
theta_xx := simplify(subs(x=X,y=Y,diff(T,x,x)),symbolic);
theta_y := simplify(subs(x=X,y=Y,diff(T,y)), symbolic);
theta_yy := simplify(subs(x=X,y=Y,diff(T,y,y)),symbolic);

>

Making these substitutions into the sum , we get

> expand(simplify(subs(r[x]=rx, r[xx]=rxx, r[y]=ry, r[yy]=ryy, theta[x]=theta_x, theta[xx]=theta_xx, theta[y]=theta_y, theta[yy]=theta_yy, fxx+fyy)));

>

which we recognize as the equivalent of

> unassign('f');
expand(laplacian(f(r,theta),[r,theta], coords=polar));

>