Enns R.H. Computer Algebra Recipes for Mathematical Physics

Подождите немного. Документ загружается.

3.2. VECTORS: CURVILINEAR COORDINATES 107

curl:=Del &x A; div:=simplify(Del . A);

curl := 2 cos(φ) cos(θ)

− 2sin(θ) cos(φ) e

+3sin(θ) e

div := 5 cos(θ) − sin(φ)



A has a non-zero curl, and therefore is not conservative. Now, let’s check the

divergence theorem. The left-hand side



2π



π/2



(∇·



A) dv of the divergence

theorem is entered for the volume contained in the hemispherical bowl

LHS:=Int(Int(Int(div*rˆ2*sin(theta),r=0..R),

theta=0..Pi/2),phi=0..2*Pi);

LHS :=



2 π



sin(θ) (5 cos(θ) − sin(φ)) dr dθ dφ

andthenevaluatedwiththevalue command.

LHS:=value(LHS);

LHS :=

5 R

On the hemispherical upper surface of the bowl, the unit vector is radial every-

where, i.e., is given by ˆe

. This unit vector is entered in n1 as a vector ﬁeld. On

the planar bottom surface of the bowl, the unit vector points in the direction

ˆe

. This unit vector is entered in n2, also as a vector ﬁeld.

n1:=VectorField(<1,0,0>); n2:=VectorField(<0,1,0>);

n1 :=

n2 := e

The vector ﬁeld



A is evaluated on the top (hemisphere) surface in At and on

the bottom (planar) surface in Ab.

At:=eval(A,r=R); Ab:=eval(A,theta=Pi/2);

At := R cos(θ)

+ R sin(θ) e

+ R sin(θ) cos(φ) e

Ab := r e

+ r cos(φ) e

Now, we can evaluate the surface integral on the right-hand side of the diver-

gence theorem. The contribution from the hemispherical top is of the form



2π



π/2

(



At · ˆn1) R

sin θdθdφ, which is now entered and then evaluated.

RHS1:=Int(Int((At . n1)*Rˆ2*sin(theta),theta=0..Pi/2),

phi=0..2*Pi);

RHS1 :=



2 π



cos(θ)sin(θ) dθ dφ

RHS1:=value(RHS1);

RHS1 := R

The contribution to the surface integral from the planar bottom is given by



2π



(



Ab · ˆn2) r sin(π/2) dr dφ which is entered and evaluated.

RHS2:=Int(Int((Ab . n2)*r*sin(Pi/2),r=0..R),phi=0..2*Pi);

RHS2 :=



2 π



dr dφ

108 CHAPTER 3. VECTORS AND MATRICES

RHS2:=value(RHS2);

RHS2 :=

2 R

Adding the two surface integral contributions, RHS1 and RHS2 , completes the

evaluation of the right-hand side of the divergence theorem.

RHS:=RHS1+RHS2;

RHS :=

5 R

The results RHS and LHS are identical, conﬁrming the divergence theorem.

3.3 Matrices

A general matrix A of order m ×n is of the form

A =

⎛

⎜

⎝

... a

... ... ... ... ...

... a

⎞

⎟

⎠

the index m labeling the row and n the column. Each number a

in A is called

an element. If the number of columns equals the number of rows, the matrix

is said to be square. A matrix having only one row is called a row matrix or a

row vector, while a matrix having only one column is called a column matrix or

column vector. Some basic properties of matrices are as follows:

(a) Addition and Subtraction: If two matrices A=(a

)andB =(b

)areof

the same order, then A ± B =(a

± b

(b) Multiplication: If A=(a

)andλ is any scalar, then λA= Aλ=(λa

If A=(a

)isanm ×n matrix and B =(b

)isann ×p matrix, then the

product AB (or A · B) is a matrix C =(c

), where c



i=1

The new matrix C is of order m ×p.

of a matrix A =(a

)isA

=(a

), i.e.,

interchange rows and columns.

(d) Hermitian matrix: If A =(a

), then the complex conjugate matrix is

∗

=(a

∗

). The Hermitian conjugate (or adjoint) matrix A

†

=(A

)

∗

square matrix is Hermitian if A

†

= A.

(e) Inverse of a matrix: If A is a non-singular matrix (i.e., its determi-

nant det(A) = 0), then the inverse matrix A

−1

is of the form A

−1

)

/det(A), where (A

) is the matrix of cofactors A

.(Thecofactor

is equal to (−1)

j+k

times the resulting determinant of A obtained by

removing all the elements of the jth row and kth column.) It follows that

−1

= A

−1

A = I,whereI is the unit matrix with each element along

its principal diagonal equal to 1 and all oﬀ-diagonal elements 0.

Other matrix properties will be introduced in the recipes of this section.

3.3. MATRICES 109

3.3.1 Some Matrix Basics

The basic tool for the manipulation of reality is the manipulation of

words. If you can control the meaning of words, you can control the

people who must use the words.

Philip K. Dick, American science ﬁction writer, (1928–82)

In this recipe, the basic Maple syntax for dealing with matrices is introduced.

If A=

⎛

⎝

21−1

1 −23

−212

⎞

⎠

and B =

⎛

⎝

1 −12

−213

2 −11

⎞

⎠

show that

(a) A and B are non-singular matrices

(b) (A + B)

= A

+ AB+ BA+ B

−1

= B

−1

(d) (AB)

= B

The basic Maple library package for dealing with matrices and determinants is

the LinearAlgebra package. Replace the colon with a semicolon in with(Linear

Algebra) if you wish to see the command structures available in this package.

restart: with(LinearAlgebra):

A matrix can be entered in diﬀerent ways. For example, A is now entered using

the Matrix command with the argument consisting of the elements arranged in

order in a list of lists.

A:=Matrix([[2,1,-1],[1,-2,3],[-2,1,2]]);

A :=

⎡

⎣

21−1

1 −23

−212

⎤

⎦

The matrix is clearly a 3 × 3 (square) matrix. The dimensionality of A can be

conﬁrmed using the Dimension command.

Dimension(A);

3, 3

A matrix such as B can also be entered by using the following syntax.

B:=<<1|-1|2>,<-2|1|3>,<2|-1|1>>;

B :=

⎡

⎣

1 −12

−213

2 −11

⎤

⎦

(a) To prove that A and B are non-singular matrices, their determinants are

calculated and seen to be non-zero.

Determinant(A); Determinant(B);

−19 − 4

(b) On the left-hand side, (A + B)

is calculated.

110 CHAPTER 3. VECTORS AND MATRICES

LHS:=(A+B)ˆ2;

LHS :=

⎡

⎣

90 6

−2111

00 9

⎤

⎦

To multiply two matrices A and B, one can use either the “long-hand” syntax

Multiply(A,B) or the shorter dot notation A.B, again leaving spaces either

side of the dot as we did for vector multiplication. The right-hand side A

AB + BA+ B

is now entered using the “short-hand” syntax. For clarity, I

have place round brackets around each matrix multiplication.

RHS:=Aˆ2+(A . B)+(B . A)+Bˆ2;

RHS :=

⎡

⎣

90 6

−2111

00 9

⎤

⎦

The matrix result in RHS isthesameasinLHS , thus proving the relation.

−1

and B

−1

are determined.

Ainverse:=MatrixInverse(A); Binverse:=MatrixInverse(B);

Ainverse :=

⎡

⎢

⎣

−1

−2

⎤

⎥

⎦

Binverse :=

⎡

⎢

⎣

−1

−2

⎤

⎥

⎦

Then (AB)

−1

− B

−1

is calculated in eq1 .

eq1:=MatrixInverse(A . B)-(Binverse . Ainverse);

eq1 :=

⎡

⎣

000

⎤

⎦

The result is a zero or null matrix, thus conﬁrming the relation.

(d) The transposes A

and B

of A and B are determined.

Atranspose:=Transpose(A); Btranspose:=Transpose(B);

Atranspose :=

⎡

⎣

21−2

1 −21

−132

⎤

⎦

Btranspose :=

⎡

⎣

1 −22

−11−1

231

⎤

⎦

Calculating (AB)

−B

in eq2 yields a null matrix, conﬁrming the relation.

eq2:=Transpose(A . B)-(Btranspose . Atranspose);

eq2 :=

⎡

⎣

000

⎤

⎦

3.3. MATRICES 111

3.3.2 Eigenvalues and Eigenvectors

A man can become so accustomed to the thought of his own

faults that he will begin to cherish them as charming little

“personal characteristics”.

Helen Rowland, American journalist, A Guide to Men,“Brides” (1922)

Let A be an n × n matrix and X a column vector with n elements. The

equation AX=λX,or(A −λI) X =0,whereλ is a number and I the unit or

identity matrix, has non-trivial solutions if and only if det(A −λI)=0. When

expanded, the determinant yields a polynomial equation of degree n in λ, called

the characteristic polynomial equation. The roots of this polynomial equation

are called the characteristic values or eigenvalues of the matrix A. Correspond-

ing to each eigenvalue, there will be a non-trivial (X = 0) solution X,whichis

called a characteristic vector or eigenvector.

Find the eigenvalues and corresponding eigenvectors of the matrix

A :=

⎛

⎝

57−5

04−1

28−3

⎞

⎠

Determine a set of unit eigenvectors. A unit eigenvector has unit length, i.e.,

the sum of the squares of its components is 1.

To solve this problem, let’s load the LinearAlgebra package and enter

the given matrix A.

restart: with(LinearAlgebra):

A:=<<5|7|-5>,<0|4|-1>,<2|8|-3>>;

A :=

⎡

⎣

57−5

04−1

28−3

⎤

⎦

This problem can be tackled in two equivalent ways, either by mimicking the

hand calculation or using speciﬁc Maple commands which accomplish the same

task. Both ways will be demonstrated. The 3 × 3 identity matrix is entered.

IM:=IdentityMatrix(3);

IM :=

⎡

⎣

100

010

001

⎤

⎦

The identity matrix is multiplied by λ in S.Thesimplify command must be

applied to actually accomplish the multiplication. The same result could be

alternatively achieved by entering ScalarMatrix(lambda,3,3).

S:=simplify(lambda*IM);

S :=

⎡

⎣

λ 00

0 λ 0

00λ

⎤

⎦

112 CHAPTER 3. VECTORS AND MATRICES

If we were proceeding by hand, the characteristic matrix A − λI would be

formed. This is illustrated in CM .

CM:=A-S;

CM :=

⎡

⎣

5 − λ 7 −5

04− λ −1

28−3 − λ

⎤

⎦

The characteristic matrix may be more simply obtained directly from A by

using the CharacteristicMatrix command as in CM2 .

CM2:=CharacteristicMatrix(A,lambda);

CM2 :=

⎡

⎣

5 − λ 7 −5

04− λ −1

28−3 − λ

⎤

⎦

Proceeding by hand, we then would determine the characteristic polynomial

equation by taking the determinant of the characteristic matrix and setting it

equal to zero. This is done in CP, a cubic polynomial in λ resulting.

CP:=Determinant(CM)=0;

CP := 6 − 11 λ +6λ

− λ

The same equation follows on applying the CharacteristicPolynomial com-

mand directly to A and setting the result equal to zero.

CP2:=CharacteristicPolynomial(A,lambda)=0;

CP2 := λ

− 6 λ

+11λ − 6=0

In the hand calculation, the eigenvalues are obtained by solving the character-

istic polynomial equation CP for λ, 3 eigenvalues resulting in L.

L:=solve(CP,lambda);

L := 1, 2, 3

The same three eigenvalues again are more simply obtained by applying the

Eigenvalues command to A. Unless otherwise speciﬁed, the eigenvalues appear

by default in a column format.

Eigenvalues(A);

⎡

⎣

⎤

⎦

Instead of a column format, including the option output=list will produce the

eigenvalues in a list format as shown in the following command line.

Eigenvalues(A,output=list);

[1, 2, 3]

Now that the eigenvalues are determined, the corresponding eigenvectors will

be calculated. If proceeding by hand, the eigenvalues are extracted separately

from L and labeled L1 , L2 ,andL3 .

L1:=L[1]; L2:=L[2]; L3:=L[3];

3.3. MATRICES 113

L1 := 1 L2 := 2 L3 := 3

A general eigenvector X, with three elements, is then entered as a column

vector. Note that the elements are separated by commas when inputting a

column vector. The goal is to determine the value of x1, x2, and x3 for each of

the eigenvalues.

X:=<<x1,x2,x3>>;

X :=

⎡

⎣

⎤

⎦

Then the column vector X is multiplied by the characteristic matrix CM ,the

latter being evaluated for a speciﬁc eigenvalue, e.g., λ =L1 .

eq1:=eval(CM,lambda=L1) . X;

eq1 :=

⎡

⎣

4 x1 +7x2 − 5 x3

3 x2 − x3

2 x1 +8x2 − 4 x3

⎤

⎦

Extracting each row entry from eq1 and setting it equal to 0 yields a set of

equations in system1 for x1, x2, and x3.

system1:={eq1[1,1]=0,eq1[2,1]=0,eq1[3,1]=0};

system1 := {4 x1 +7x2 − 5 x3 =0, 3 x2 − x3 =0, 2 x1 +8x2 − 4 x3 =0}

The system of equations is solved in sol1 for the unknowns, the answer being

expressed in terms of the arbitrary quantity x2.

sol1:=solve(system1,{x1,x2,x3});

sol1 := {x3 =3x2 , x1 =2x2 , x2 = x2 }

The eigenvector X1 corresponding to the eigenvalue L1 is now obtained by

evaluating X with sol1 .

X1:=eval(X,sol1);

X1 :=

⎡

⎣

2 x2

3 x2

⎤

⎦

There is an inﬁnite family of eigenvectors X1 depending on the value chosen for

x2 . Proceeding by hand, the other two eigenvectors could be similarly obtained.

Thechoiceofaspeciﬁcvalueforx2 will be postponed for the moment.

The multitude of steps in deriving the three eigenvalues and eigenvectors

can be completely bypassed by applying the Eigenvectors command to A.

(eiv,V):=Eigenvectors(A);

eiv,V :=

⎡

⎣

⎤

⎦

⎡

⎢

⎣

−1

11 1

⎤

⎥

⎦

114 CHAPTER 3. VECTORS AND MATRICES

In the output, the eigenvalues are presented in column format (the order of the

entries may vary from one execution of the worksheet to the next) before the

comma. Corresponding to the ﬁrst (top) entry, the corresponding eigenvector

is given by the ﬁrst column in the matrix format after the comma. The second

column corresponds to the second eigenvalue entry, and so on. So that the eigen-

values and eigenvectors can be extracted separately, the assignment (eiv,V) is

made, eiv corresponding to the eigenvalues and V to the eigenvectors.

Choosing x2=1/3, the “hand calculated” eigenvector X1 , corresponding to

the eigenvalue 1, is identical to the ﬁrst column in V .

X1:=eval(X1,x2=1/3);

X1 :=

⎡

⎢

⎣

⎤

⎥

⎦

The three eigenvectors are now extracted from V by using the Column command.

The second argument speciﬁes the column to be extracted.

V1:=Column(V,1); V2:=Column(V,2); V3:=Column(V,3);

V1 :=

⎡

⎢

⎣

⎤

⎥

⎦

V2 :=

⎡

⎢

⎣

⎤

⎥

⎦

V3 :=

⎡

⎣

−1

⎤

⎦

Finally, it was requested that a set of unit eigenvectors be obtained. The sum

of the squares of the elements of each eigenvector can be obtained by applying

the Norm(V,2) command. Dividing each of the above eigenvectors by its norm,

produces a unit eigenvector.

U1:=V1/Norm(V1,2); U2:=V2/Norm(V2,2); U3:=V3/Norm(V3,2);

U1 :=

⎡

⎢

⎣

√

⎤

⎥

⎦

U2 :=

⎡

⎢

⎣

√

⎤

⎥

⎦

U3 :=

⎡

⎢

⎣

−

√

⎤

⎥

⎦

When ﬁrst learning how to ﬁnd eigenvalues and eigenvectors, its important to

go through the steps outlined in “our hand calculation”. But, after doing this

a few times, you will be grateful for Maple’s specialized commands.

3.3. MATRICES 115

3.3.3 Diagonalizing a Matrix

An idea is a point of departure and no more. As soon as you

elaborate it, it becomes transformed by thought.

Pablo Picasso, Spanish artist, (1881–1973)

Find a matrix C which reduces the matrix

A =

⎛

⎝

24−6

42−6

−6 −6 −15

⎞

⎠

to diagonal form (all elements not on the main diagonal are zero) by the trans-

formation C

−1

AC.

After loading the LinearAlgebra package, the matrix A is entered.

restart: with(LinearAlgebra):

A:=<<2|4|-6>,<4|2|-6>,<-6|-6|-15>>;

A :=

⎡

⎣

24−6

42−6

−6 −6 −15

⎤

⎦

A general matrix C with doubly subscripted elements is entered.

C:=<<c[1,1]|c[1,2]|c[1,3]>,<c[2,1]|c[2,2]|c[2,3]>,

<c[3,1]|c[3,2]|c[3,3]>>;

C :=

⎡

⎣

1, 1

1, 2

1, 3

2, 1

2, 2

2, 3

3, 1

3, 2

3, 3

⎤

⎦

The eigenvalues of A are obtained.

L:=Eigenvalues(A);

L :=

⎡

⎣

−2

−18

⎤

⎦

L gives the three distinct eigenvalues. It is desired to ﬁnd a matrix C which

reduces A to a diagonal matrix S with the main diagonal elements given by the

eigenvalues and the oﬀ-diagonal elements equal to zero. The DiagonalMatrix

command is used to form S.

S:=DiagonalMatrix([L[1],L[2],L[3]],3,3);

S :=

⎡

⎣

−20 0

09 0

00−18

⎤

⎦

Now, we want C

−1

AC=S. Multiplying this matrix equation from the left by

C and noting that CC

−1

= I,wehaveIAC= AC= CS,orAC − CS=0.

The left-hand side of this last matrix equation is now entered in E.

E:=(A . C)-(C . S);

116 CHAPTER 3. VECTORS AND MATRICES

E :=



4 c

1, 1

+4c

2, 1

− 6 c

3, 1

, −7 c

1, 2

+4c

2, 2

− 6 c

3, 2

, 20 c

1, 3

+4c

2, 3

− 6 c

3, 3





4 c

1, 1

+4c

2, 1

− 6 c

3, 1

, 4 c

1, 2

− 7 c

2, 2

− 6 c

3, 2

, 4 c

1, 3

+20c

2, 3

− 6 c

3, 3





− 6 c

1, 1

− 6 c

2, 1

− 13 c

3, 1

, −6 c

1, 2

− 6 c

2, 2

− 24 c

3, 2

, −6 c

1, 3

− 6 c

2, 3

+3c

3, 3



We have to solve the system of equations resulting from setting the 9 entries

of E equal to zero. This system of 9 equations is now obtained by using the

sequence command twice to sum over the rows and columns of E.

System:={seq(seq(E[i,j]=0,j=1..3),i=1..3)};

System := {4 c

1, 1

+4c

2, 1

− 6 c

3, 1

=0, −7 c

1, 2

+4c

2, 2

− 6 c

3, 2

=0,

20 c

1, 3

+4c

2, 3

− 6 c

3, 3

=0, 4 c

1, 2

− 7 c

2, 2

− 6 c

3, 2

=0,

4 c

1, 3

+20c

2, 3

− 6 c

3, 3

=0, −6 c

1, 1

− 6 c

2, 1

− 13 c

3, 1

=0,

−6 c

1, 2

− 6 c

2, 2

− 24 c

3, 2

=0, −6 c

1, 3

− 6 c

2, 3

+3c

3, 3

=0}

The double sequence command is used again to enter the 9 unknown matrix

elements to be solved for.

Elements:={seq(seq(c[i,j],j=1..3),i=1..3)};

Elements := {c

1, 1

1, 2

1, 3

2, 1

2, 2

2, 3

3, 1

3, 2

3, 3

}

The system of 9 equations is solved for the 9 unknown matrix elements.

Solution:=solve(System,Elements);

Solution := {c

1, 3

= c

2, 3

3, 3

=4c

2, 3

1, 2

= −2 c

3, 2

1, 1

= −c

2, 1

= c

2, 1

2, 3

= c

2, 3

3, 2

= c

3, 2

3, 1

=0,c

2, 2

= −2 c

3, 2

}

Then Solution is used to evaluate C.

C:=eval(C,Solution);

C :=

⎡

⎣

−c

2, 1

−2 c

3, 2

2, 3

2, 1

−2 c

3, 2

2, 3

0 c

3, 2

4 c

2, 3

⎤

⎦

The matrix elements c

2,1

, c

3,2

,andc

2,3

are undetermined. We are free to

choose values for these elements. There are an inﬁnity of choices and therefore

of C matrices which will diagonalize A with the transformation C

−1

AC.Let’s

check that this is so, keeping the three elements unspeciﬁed for the moment.

The inverse matrix C

−1

is calculated.

Cinverse:=MatrixInverse(C);

Cinverse :=

⎡

⎢

⎣

−

2, 1

−

3, 2

−

3, 2

2, 3

⎤

⎥

⎦