Linear algebra in physics is the art of using matrices and vectors to solve physical problems. Maple knows a lot about linear algebra and you can explore what it knows by typing ?LinearAlgebra. The sections below just give a brief introduction, but before you start doing any of this stuff you have to tell Maple you want it to use linear algebra, like this:

> with(LinearAlgebra):

Matrices, vectors, and basic operations

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

To refer to the elements of this matrix use the notation A[row,column] , i.e. A[1,2] like this:

> A[1,2];

And if you want to just define a 2x2 matrix B but not load values yet you would use

> B:=Matrix(1..2,1..2);

Once B is defined you can load its elements any way you want, including by hand, like this

> B[1,1]:=2;B[1,2]:=3;B[2,1]:=7;B[2,2]:=4;

Another one of Maple's ways of entering matrices is so handy that it should be mentioned here. Click View on the toolbar, then select Palettes , and then Matrix Palettes . You will now be looking at a template for selecting various sizes of matrices. Click on the size you want and the following Maple command will appear, just waiting for you to enter matrix elements in place of the question marks and then to assign it a name with := .

> Matrix([[%?, %?, %?], [%?, %?, %?], [%?, %?, %?]]);

You might want to look at B to check that it is right, and try this

> B;

For more details about the Matrix command use ?Matrix.

The default is to make a column vector

> b:=Vector([5,6]);

To make a row vector you have to specifically tell Maple, like this

> b:=Vector[row]([5,6]);

> c:=Column(A,2);

To load c with the first row of A use

> c:=Row(A,1);

There is also a handy short syntax for defining matrices and column vectors. Column vectors can be defined without a special command using the <...> syntax:

> s:=<1,2,3>;

And since a matrix can be thought of as a collection of column vectors, matrices can be defined similarly

> R:=<<1,2,3>|<4,5,6>|<7,8,10>>;

This short form is sometimes used in the examples in online help so I thought you should know about it.

Since this is Maple, you don't have to just work with numbers. Your matrices and vectors can also be formulas, like this

> d:=Matrix([[x11,x12],[x21,x22]]);

or

> e:=Matrix([[x^2*y,y*x],[sin(x),1/x]]);

Now suppose you want to give x and y values and then see what the matrices are

> x:=1;y:=2;

> e;

Well, we've been here before and we know what to do--just use evalf :

> evalf(e);

Hmmm--apparently not. The problem is that e is not an expression but a matrix of expressions, so Maple needs to use two eval commands: evalm on the inside to refer directly to the elements of e and evalf on the outside to force floating point evaluation.

> evalf(evalm(e));

Use ?evalm to see what this command is about, or just remember to use decimal points on your numbers when you want a numerical answer to come out.

>

> B:=Matrix([[1,2,3],[4,5,6],[7,8,9]]);

> C:=Matrix([[3,2,1],[5,6,4],[9,7,8]]);

> E:=Add(B,C);

Subtraction is just adding with a minus sign

> E:=Add(B,-C);

You do have to be careful that addition and subtraction make sense. For instance, you can't do this

> F:=Matrix([[5,7],[1,3]]);

> Add(E,F);

Error, (in LinearAlgebra:-MatrixAdd) matrix dimensions don't match: 3 x 3 vs 2 x 2

because the matrices are not the same size.

> G:=Multiply(B,C);

There is a wordier way to do matrix multiplication which you can find in the command list by going to ?LinearAlgebra.

Problem 5.1

Verify that Maple did this right by computing by hand.

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

You can also multiply matrices and vectors together, like this

> b:=<5,6,7>;

> Multiply(B,b);

As expected, a square matrix times a column vector is a column vector. A row vector times a square matrix would be a new row vector, and Maple knows how to do this too:

> c:=Vector[row]([1,2,3]);

(Or you could use the shorthand form)

> c:=<1|2|3>;

> Multiply(c,C);

You can also multiply matrices and vectors by scalars using the ordinary multiplication operator

> 3*b;

> 5*B;

>

Finally, you will notice that division is not mentioned here. That's because division with matrices and vectors is rather tricky and actually falls under the heading of "Solving systems of linear equations", which comes a little later.

Problem 5.2

> Rx:=Matrix([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]]);

Similarly, the rotation matrices about the y- and z-axes through angles and are given by

> Ry:=Matrix([[cos(phi),0,sin(phi)],[0,1,0],[-sin(phi),0,cos(phi)]]);

> Rz:=Matrix([[cos(psi),-sin(psi),0],[sin(psi),cos(psi),0],[0,0,1]]);

>

(Note: there is not a sign mistake on the sines in the Ry matrix. The sine signs are different in this matrix than they are in Rx and Rz because rotation is always defined with respect to the right-hand rule on angle: put your thumb in the direction of the rotation axis and your fingers curl in the direction of positive angle. Try it with your hand with and check that Ry is correct.) Use these rotation matrices to check where the vector [1,0,0] ends up after the following sequence of rotations: (a) by about the y-axis, then (b) by about the x-axis, and finally (c) by about the z-axis. After Maple gives you the answer verify that it is correct by making a coordinate system with your fingers and doing the rotations. (Make sure the coordinate system is right-handed). Then start with the same vector [1,0,0] but do the rotations in the reverse order. Do you end up in the same place? (You don't--we describe this situation by saying that rotation operations don't commute.) Rotation should not change the magnitude of the vector; does this sequence of matrix multiplications preserve the magnitude of the vector? You can get the magnitude of a vector v with the command Norm(v,2) .

> MatrixInverse(C);

To check it do this

> Multiply(C , MatrixInverse(C));

> Determinant(B);Determinant(C);

> Transpose(C);

> E:=Matrix([[1,2,3]]);

> Transpose(E);

And in quantum mechanics you will want to take the transpose of the complex conjugate, which is called the HermitianTranspose :

> F:=<<1,I>|<2*I,2>>;

> HermitianTranspose(F);

> Trace(C);

> Norm(E,2);

>

For matrices the norm is more complicated; look under ?LinearAlgebra[Norm] .

Problem 5.3

Find the inverse of the matrix B defined in this section. Explain why you get an error message instead of an answer.

[Alice, Bob, Chuck]. Here are A and b:

> restart;with(LinearAlgebra):

> A:=Matrix([[1,1,-3/2],[0,-1,1],[1,-2-1/4,1]]);

> b:=Vector([0,6,0]);

We would like to solve for the vector x of the friends' ages. Maple has a command LinearSolve that makes this kind of problem really easy:

> LinearSolve(A,b);

>

And that's all there is to it. Any linear system can be solved this way. We will learn later that Maple's equation solving command solve can also do problems like this.

Problem 5.4

Consider 5 springs with spring constants hooked together along the x-axis in this order. If each spring has known relaxed length , find the unknown lengths of each spring when the springs are stretched until the sum of their lengths is twice their total unstretched length. Have Maple give you the general algebraic solution to this problem (it is moderately horrible), then get a numerical answer using . Let all the relaxed lengths be the same and equal to 10 cm. Recall that the force exerted by a spring is given by the formula . You will also need to remember how to do equilibrium problems from freshman mechanics.

> restart:with(LinearAlgebra):

> v1:=<1|2|3>;v2:=<4|5|6>;

> DotProduct(v1,v2);CrossProduct(v1,v2);

> VectorAngle(v1,v2);

or, rather,

> evalf(VectorAngle(v1,v2));

>

Note that it comes out in radians.

>

The eigenvalue problem as applied to matrices is this: given a square matrix A , what are the vectors and numbers such that when A is multiplied into the result is the same vector , but multiplied by , or in other words, such that A ? Since matrix multiplication by a vector simultaneously rotates the vector and multiplies it by a factor, it is not clear that there is any solution at all to this problem. But it can be proven mathematically that in most cases, if the matrix is NxN, then there are N such vectors and values (called eigenvectors and eigenvalues). But they may be complex, which sometimes (but not always) means they are unphysical. Here's an example.

> restart:with(LinearAlgebra):

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

> Eigenvectors(A);

Ok, it's clear we've got to do something about this. We usually just want numbers, not gigantic algebraic expressions. An easy way to get them (for the 142nd time) is to enter the elements of the matrix as floating-point numbers, i.e., they have to have decimal points, like this

> A:=Matrix([[1.,2.,3.],[4.,3.,2.],[-1.,0.,2.]]);

> v:=Eigenvectors(A);

Good--much better. Now look at the output. According to online help for Eigenvectors the first set of numbers is a column vector containing the eigenvalues, and the second set of numbers is a matrix whose columns are the eigenvectors. So if we wanted to study the third eigenvector we would start by picking off the third eigenvalue like this

> lambda3:=v[1][3];

(I know this looks weird, but since v[1] is a column vector it makes sense to ask for its third element). Now we want the third column of v[2], a job tailor made for Maple's Column command

> v3:=Column(v[2],3);

Now we are ready to see if this really is an eigenvector by multiplying by A . If it works, then A will just give us back.

> y:=Multiply(A ,v3);

To see if it got it right let's divide y by and subtract from it to see if we get zero

> y/lambda3-v3;

which is exactly right, to the precision of Maple floating point numbers.

>

Problem 5.5

Here's a simple eigen-idea that shows up in numerical analysis all the time. Think about what happens if you pick a vector and just keep multiplying it by A in the example above all the time. Your original vector could be written as a linear combination of the eigenvectors (usually), so you can think about what multiplication by A does to each of the three eigenvector components. Now look, this is not too hard--eventually, the vector with the largest eigenvalue will completely take over because every time you multiply by A you get more of it back than you get of the other two. In the example above the biggest eigenvalue is about twice its nearest competitor, and you already know what repeatedly multiplying by 2 can do. So now do this: choose a random vector; then multiply it by A and divide it by its norm using Norm(v,2) so that the result is a unit vector; then do the multiplication and renormalization about 20 times. Now find the angle between your final result and the eigenvector with the largest eigenvalue; did your calculation get this eigenvector? Can you explain the angle you get?

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

Problem 5.6

Go get the rotation matrix Rx from Problem 5.2 and have Maple find its eigenvalues and eigenvectors. Then think about what rotation means and see if you can make physical sense out of Maple's answers. Maple will require you to choose a value for theta, so try (a) , (b) , and (c) . (Note: you will encounter (c) when you study spin in quantum mechanics.) Warning: Due to a bug in Maple 6 you will need to do each part of this problem in this order: (a) define the angle (b) build Rx (c) find the eigenvectors. Or you could define the matrix as a Maple function like this:

> M:=theta->Matrix([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]]);

and then use this function to define the three matrices in parts (a)-(c).

>

Cartesian:

> f1:=3*x^2+y*sqrt(z);

> grad(f1,[x,y,z]);

Cylindrical

> f2:=cos(theta)/sqrt(r^2+z^2);

> grad(f2,[r,theta,z],coords=cylindrical);

Spherical

> f3:=sin(phi)*r^2*(3*cos(theta)^2-1);

> grad(f3,[r,theta,phi],coords=spherical);

>

Cartesian

> g1:=[x^2,y^2,z^2];

> diverge(g1,[x,y,z]);

Cylindrical

> g2:=[r^2*cos(theta),r^2*sin(theta),z^2];

> diverge(g2,[r,theta,z],coords=cylindrical);

Spherical

> g3:=[r^2*cos(theta),r^2*sin(theta),cos(phi)];

> diverge(g3,[r,theta,phi],coords=spherical);

>

Cartesian

> curl(g1,[x,y,z]);

Cylindrical

> curl(g2,[r,theta,z],coords=cylindrical);

Spherical

> curl(g3,[r,theta,phi],coords=spherical);

>

Cartesian

> laplacian(f1,[x,y,z]);

Cylindrical

> laplacian(f2,[r,theta,z],coords=cylindrical);

Spherical

> laplacian(f3,[r,theta,phi],coords=spherical);

>

This capability might be nice, but usually in physics the problem is not calculating these vector derivatives, but in undoing them to find a field. This is much harder and Maple is not much help in this regard.

Problem 5.7

Two of Maxwell's equations for static electric and magnetic fields are and . For each of the following E or B fields use Maple's vector calculus routines to find formulas for the charge density and the current density . Note: and are vectors, so you need to define them properly to do this problem. For instance, if then you need to define this way: B:=array([0,k*r,0]);

(a) Spherical coordinates: (point charge) (b) Spherical coordinates: (solid sphere of charge)

(c) Cylindrical coordinates: (line charge) (d) Cylindrical coordinates: (solid cylinder of charge)

(e) Cylindrical coordinates: (long wire) (f) Cylindrical coordinates: (uniform current in a wire)

(g) Cylindrical coordinates:

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

Problem 5.8

The electrostatic potential of a point dipole is given in spherical coordinates by

. Find the electric field of this dipole in spherical coordinates by using the defining relation for the potential . What is the power law for the way E falls off with distance from the origin?