Examples in Physics - Maple Help

はじめる前に
新機能一覧
Maple ワークシートの作成
Mapleワークシートを共有
Maple ウィンドウのカスタマイズ
Maple システムのカスタマイズ
コネクティビティ
Mathematics
数学
物理
- Physics パッケージ
プログラミング
グラフィックス
Student パッケージ
Science and Engineering
科学・エンジニアリング
アプリケーションと例題
- アプリケーション
- 概要
- 例題
リファレンス
システム
MapleSim
MapleSim ツールボックス
MapleSim ヘルプ
Tasks
Toolboxes
エラーメッセージガイド
マニュアル
数学アプリ

Examples in Physics

Vectors and Analytical Geometry

Vectorial equation of a line

Problem

1.	Derive the vectorial equation of a straight line passing through two generic points A and B.

2.	Choose two concrete points A and B, and use that vectorial equation to obtain the parametric equation of the line.

3.	Generate a 3D plot of this line.

Solution

The vectorial equation of a line is the equation satisfied by the position vector of any point of this line; in this case we want to express that equation in terms of the position vectors of the points A and B. The parametric equation is obtained from the vectorial equation with the coefficients of the unit vectors.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(1.1.1)

1. To construct the vectorial equation of the line, let and be vectors pointing to A and B. The line we want to represent with a vectorial equation is then parallel to the vector difference

A_ - B_;

(1.1.2)

Now let be the position vector of any point of this line. Then the vector is also parallel to the line and so to ; therefore

>	Eq := (r_ - A_) &x (1.1.2) = 0;

(1.1.3)

This is the vectorial equation of the line passing through A and B. The position vector of any point belonging to this line satisfies this equation.

2. To construct the parametric representation for this line, turn the generic points A and B into concrete points; that is, give values to the components of and . For example:

>	A_ := 2_i - 3_j + 4* _k;

(1.1.4)

>	B_ := -4* _i + 2*_j - _k;

(1.1.5)

Regardless of the values of and , for the position vector of a point of the line, we always have

>	r_ := x_i + y_j +z*_k;

(1.1.6)

Thus the vectorial equation for these particular points A and B becomes

Eq;

(1.1.7)

Since the unit vectors [] are independent, this vectorial equation represents three equations, constituted by the coefficients of these unit vectors

>	map(`=`, [coeffs(lhs(Eq), [_i, _j, _k])], 0);

(1.1.8)

This is a system of three equations in terms of three unknowns x, y, and z, but by construction is the equation of a line, so two of the variables can be expressed in terms of the third one, acting as the parameter of the parametric equations; you get these parametric equations by solving the system.

>	solve((1.1.8), [x, y, z]);

(1.1.9)

The interpretation is: given a value for z, the values of x and y are those you see above.

3. To plot this parametric representation (see plots[spacecurve]), first construct the appropriate input; select the curve's parameter:

>	select(evalb, (1.1.9)[1]);

(1.1.10)

So the input for the spacecurve command, say for z from -4 to 4, is

>	input := [op(map(rhs,(1.1.9)[1]))], lhs((1.1.10)[1]) = -4..4;

(1.1.11)

Some options improving the visualization are

>	opts := axes = boxed, scaling = constrained, orientation = [-130, 70], labels = [x,y,z];

(1.1.12)

>	plots[spacecurve](input, opts);

Vectorial equation of a plane

Problem

1.	Derive the vectorial equation of a plane passing through three generic points A, B, and C.

2.	Choose three concrete points A, B, and C and plot this plane.

Solution

The vectorial equation of a plane is the equation satisfied by the position vector of any point of this plane.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(1.2.1)

1. To construct this equation, let , , and be the vectors pointing to A, B, and C, respectively, and be the position vector of any point of this plane. The differences and are vectors parallel to the plane we want to represent, and also the difference between and any of , or is a vector parallel to the plane. So the equation of the plane can be obtained by taking the cross product of differences involving , and to construct a vector perpendicular to the plane,

>	G_ := (A_ - B_) &x (A_ - C_);

(1.2.2)

then equating to zero the scalar product of with any of the differences parallel to the plane involving . So for instance, one way of writing this vectorial equation of the plane is

>	Eq := (r_ - A_) . G_ = 0;

(1.2.3)

2. To plot this plane, turn the generic points A, B, and C into concrete points; that is, give values to the components of , and . For example:

>	A_ := 2_i - 3_j + _k*4;

(1.2.4)

>	B_ := 5* _i + 4_j - 7_k;

(1.2.5)

>	C_ := 30/4* A_ + 90/7*B_;

(1.2.6)

For we always have

>	r_ := x_i + y_j +z*_k;

(1.2.7)

The vectorial equation for these particular points A, B, and C is thus:

Eq;

(1.2.8)

As a verification that the surface represented by this equation contains the points A, B, and C, you can substitute in the values of the coordinates of these points and see that the equation is satisfied. These are the coordinates of the three points:

>	A, B, C := seq([seq(coeff(v, ui), ui=[_i, _j, _k])], v = [A_, B_, C_]);

(1.2.9)

This shows that at these points, the equation is satisfied.

>	for P in [A, B, C] do eval(Eq, [x = P[1], y = P[2], z = P[3]]); od;

(1.2.10)

That this surface is a plane is clear from the fact that Eq is linear in all of x, y, and z. One way of plotting this plane is to use the command implicitplot3d.

>	opts := axes = boxed, scaling = constrained, orientation = [125, 65], style = surface;

(1.2.11)

>	plots[implicitplot3d](Eq, x=-2..2, y=-2..2, z=-2..2, opts);

Vectorial equation of a plane tangent to a sphere of radius a

Problem

Derive the vectorial equation of a plane tangent to a sphere of radius a.

Solution

The vectorial equation of this plane is the equation satisfied by the position vector of any point of it.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(1.3.1)

Let represent the position vector of any point of the plane, be a vector pointing to the center of the sphere, and be a vector pointing to the point B where the plane is tangent to the sphere. So the difference is a vector on the plane, and the difference is a vector from the center of the sphere to the point of contact between the sphere and the tangent plane (that is, a vector perpendicular to the plane). Hence, these two vectors are perpendicular, and so their scalar product is equal to zero.

>	Eq := (r_ - B_) . (B_ - A_) = 0;

(1.3.2)

This is already the vectorial equation of the tangent plane but not yet expressed in terms of the radius of the sphere. Now since is a vector from the center of the sphere to the point of contact with the plane, the norm of this vector is the radius a of the sphere.

>	key := Norm(B_ - A_) = a;

(1.3.3)

Expand both expressions to use them together.

>	expand(Eq);

(1.3.4)

>	map(u -> u^2, key);

(1.3.5)

>	expand((1.3.5));

(1.3.6)

So simplify one with respect to the other one, eliminating (see simplify/siderels)

>	simplify((1.3.4), {expand((1.3.5))}, {Norm(B_)});

(1.3.7)

Collecting terms in the above, this vectorial equation requested can be rewritten more compactly as . As an exercise, consider choosing three concrete values for the positions of A, B, and the radius a of the sphere, then insert these values into the equation derived, and plot the sphere and this tangent plane together (see plots[display] to merge the plots).

Volume element of a sphere

Problem

Determine the infinitesimal volume element of a sphere expressed in spherical coordinates.

Solution

Let be the vectorial equation, parameterized by u, v, and w, of a generic 3-D geometric object; in this case, we are dealing with a sphere of generic radius r. The volume element is derived from equation as .

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(1.4.1)

We want this volume element expressed in spherical coordinates (); we can always choose these coordinates themselves as parameters u, v, w. We are thus interested in the explicit form of (note the use of %diff, the inert form diff):

>	answer := %diff(r_,r) . (%diff(r_,theta) &x %diff(r_,phi));

(1.4.2)

The first step is to write the vectorial equation for a sphere of radius r; that is, the equation satisfied by the position vector of any point of the sphere. In spherical coordinates, and choosing the origin of the reference system at the center of the sphere, this vectorial equation has its simplest form:

>	r_ = R_(r, phi, theta);

(1.4.3)

>	r_ := r * _r;

(1.4.4)

where points to any point of the sphere, r is the radial coordinate (constant over the sphere) and is the radial unit vector. Now compute the partial derivatives entering , and in doing so, take into account that depends on and To make this dependence of explicit, we change the basis in the vectorial equation to the Cartesian basis (), where all the unit vectors are constant and so the partial derivatives can be performed directly.

>	r_ := ChangeBasis(r_, 1);

(1.4.5)

So, the answer introduced lines above becomes:

answer;

Physics:-Vectors:-`.`(%diff(r*sin(theta)*cos(phi)*_i+r*sin(theta)*sin(phi)*_j+r*cos(theta)*_k, r), Physics:-Vectors:-`&x`(%diff(r*sin(theta)*cos(phi)*_i+r*sin(theta)*sin(phi)*_j+r*cos(theta)*_k, theta), %diff(r*sin(theta)*cos(phi)*_i+r*sin(theta)*sin(phi)*_j+r*cos(theta)*_k, phi)))

(1.4.6)

>	value((1.4.6));

(1.4.7)

Hence the volume of element requested is

Mechanics

Static: reactions of planes and tensions on cables

Problem

A bar AB of weight w and length L has one extreme on a horizontal plane and the other on a vertical place, and is kept in that position by two cables AD and BC. The bar forms an angle with the horizontal plane and its projection BC over this plane forms an angle with the vertical plane. The cable BC is in the vertical plane where the bar is. Determine the reactions of the planes at A and B and the tensions on the cables.

Solution

There are two equations that contain information about the state of equilibrium of a system. The first one, saying that the sum of the forces acting on the body are equal to zero, tells that the center of mass of the body is not accelerated. The second one, saying that the sum of the moments of the forces acting on the body (that is, the total torque) is zero tells that the rotation of the body around its center of masses is not changing (if it is not rotating, it stays that way). These two equations involve the reactions of the planes and the tensions on the cables, so from them we can obtain the solution to the problem. There is no friction so it is also clear that the reactions and are perpendicular to the planes, as shown in the figure, and the tensions and on the cables have direction AD and BC, respectively.

The steps to solve this problem are:

1.	Determine each force acting on the bar and its application point

2.	Equate the sum of the forces to zero.

3.	Equate the sum of the moments to zero.

4.	Solve these two vectorial equations for , , and , representing the reactions of the planes at the points of contact A and B, and the tensions of the cables attached to the bar at A and B, respectively.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(2.1.1)

The forces acting on the bar are its weight , and the reactions and tensions , , and . So the two equilibrium equations are:

>	eq[1] := w_ + R_[A] + R_[B] + T_[A] + T_[B] = 0;

(2.1.2)

>	eq[2] := r_[w] &x w_ + r_[A] &x R_[A] + r_[B] &x R_[B] + r_[A] &x T_[A] + r_[B] &x T_[B] = 0;

(2.1.3)

Set the origin and orientation of the reference system to project these vectors; any choice will do, but a good one will simplify the algebraic manipulations. We choose the origin at the point B, the vertical z axis in the direction of the reaction , so that the y axis in the direction of , and the x axis in the remaining direction, anti-parallel to . With this choice, the vectors entering and are projected as follows:

>	R_[B] := abs(R[B])*_k;

(2.1.4)

where represents the norm of to be determined,

>	r_[B] := 0;

(2.1.5)

>	R_[A] := abs(R[A])*_j;

(2.1.6)

where is to be determined. This reaction is applied to the bar at A, represented by ; its component along the x axis is obtained by projecting the segment BA onto the horizontal plane (), resulting in BC, then into the intersection of the two planes. So:

>	r_[A] := Lcos(alpha)(cos(beta)_i + sin(2Pi - beta)_j) + Lsin(alpha)*_k;

(2.1.7)

For the other vectors we have

>	T_[A] := - abs(T[A])*_i;

(2.1.8)

>	T_[B] := abs(T[B])cos(beta)_i + abs(T[B])sin(2Pi - beta)*_j;

(2.1.9)

where and are to be determined,

>	w_ := -abs(w)*_k;

(2.1.10)

>	r_[w] := r_[A]/2;

(2.1.11)

The two equilibrium equations now appear as

eq[1];

(2.1.12)

eq[2];

(1/2)*L*cos(alpha)*sin(beta)*abs(w)*_i+(1/2)*L*cos(alpha)*cos(beta)*abs(w)*_j-L*sin(alpha)*abs(R[A])*_i+L*cos(alpha)*cos(beta)*abs(R[A])*_k-L*sin(alpha)*abs(T[A])*_j-L*cos(alpha)*sin(beta)*abs(T[A])*_k = 0

(2.1.13)

These two vectorial equations represent a system of six equations, obtained by equating each of the coefficients of , , and in each of the equations to zero; that is, taking the components of the vectorial equations along each axis:

>	Eq[1,2,3] := seq(Component(lhs(eq[1]), n) = 0, n = 1..3);

(2.1.14)

>	Eq[4,5,6] := seq(Component(lhs(eq[2]), n) = 0, n = 1..3);

(1/2)*L*cos(alpha)*sin(beta)*abs(w)-L*sin(alpha)*abs(R[A]) = 0, (1/2)*L*cos(alpha)*cos(beta)*abs(w)-L*sin(alpha)*abs(T[A]) = 0, L*cos(alpha)*cos(beta)*abs(R[A])-L*cos(alpha)*sin(beta)*abs(T[A]) = 0

(2.1.15)

So the system of equations to be solved is

>	sys := {Eq[1,2,3], Eq[4,5,6]}:

The unknowns are

>	var := {abs(R[A]), abs(R[B]), abs(T[A]), abs(T[B])};

(2.1.16)

and the solution is

>	solve(sys, var);

(2.1.17)

Work

Problem

A particle is submitted to a force whose Cartesian components are given by. Calculate the work done by this force when moving the particle along a straight line from the origin to a point ().

Solution

The work to be calculated is equal to the line integral of the force along the path (line) indicated:

where represents the integration domain.The steps to solve this problem are:

1.	Determine the vectorial (parameter t) equation of the line over which we are going to integrate.

2.	Express x, y, and z entering the components of in terms of the components of , then the scalar product becomes dt and the work can be expressed as .

3.	Compute the integral.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(2.2.1)

The work to be calculated is given by

>	W := Int(F_.v_, t=0..t);

Int(Physics:-Vectors:-`.`(F_, v_), t = 0 .. t)

(2.2.2)

The force acting on the particle is

>	F_ := (ax^3 + bxy^2 + cz)_i + (ax^3 + bxy^2)_j + cz*_k;

(2.2.3)

The vectorial equation of a line where t is a parameter, is represented by the vector position of any of its points. This line is passing through the origin and a generic point

>	r_[0] := x0_i + y0_j + z0*_k;

(2.2.4)

In such a case, can be written directly as the product of the parameter t and the vector representing .

>	r_ := t * r_[0];

(2.2.5)

So at the origin, t = 0, and when , t = 1. To construct the scalar product first shall be expressed in terms of the vectorial equation of this line; that is, x, y, and z entering its components shall be taken from the components of the equation .

>	x = Component(r_, 1), y = Component(r_, 2), z = Component(r_, 3);

(2.2.6)

>	F_ := eval(F_, [(2.2.6)]);

(2.2.7)

Compute now .

>	v_ := diff(r_, t);

(2.2.8)

So the line integral giving the work W is now

Int(a*t^3*x0^4+b*t^3*x0^2*y0^2+c*t*z0*x0+a*t^3*x0^3*y0+b*t^3*x0*y0^3+c*t*z0^2, t = 0 .. t)

(2.2.9)

When the particle is at , we have that t = 1, so the value of W to move the particle to is

>	eval(W, t=1);

Int(a*t^3*x0^4+b*t^3*x0^2*y0^2+c*t*z0*x0+a*t^3*x0^3*y0+b*t^3*x0*y0^3+c*t*z0^2, t = 0 .. 1)

(2.2.10)

>	value((2.2.10));

(2.2.11)

Lagrangian for a pendulum

Problem

Determine the Lagrangian of a plane pendulum having a mass m in its extremity and whose suspension point:

a) moves uniformly over a vertical circumference with a constant frequency

b) oscillates horizontally on the plane of the pendulum according to .

Solution

a) A figure for this part of the problem is

The Lagrangian is defined as

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(2.3.1)

>	L := T - U;

(2.3.2)

where T and U are the kinetic and potential energy of the system, respectively, in this case constituted by a single point of mass m. The potential energy U is the gravitational energy

>	U := -mgy;

(2.3.3)

where g is the gravitational constant and we choose the y axis along the vertical, pointing downwards, so that the gravitational force . The kinetic energy is:

>	T := 1/2m(v_ . v_);

(2.3.4)

To compute this velocity, the position vector of the suspension point of the pendulum,

>	r_ := x_i + y_j;

(2.3.5)

must be determined. Choosing the x axis horizontally and the origin of the reference system at the center of the circle (see figure above), the x and y coordinates are given by:

>	parametric_equations := [x = acos(omegat) + l* sin(phi(t)), y = -asin(omegat)+l*cos(phi(t))];

(2.3.6)

>	r_ := eval(r_, parametric_equations);

(2.3.7)

>	v_ := diff(r_, t);

(2.3.8)

(1/2)*m*((-a*sin(omega*t)*omega+l*cos(phi(t))*(diff(phi(t), t)))^2+(-a*cos(omega*t)*omega-l*sin(phi(t))*(diff(phi(t), t)))^2)

(2.3.9)

This expression contains products of trigonometric functions, so one simplification consists of combining these products.

>	combine(T, trig);

(1/2)*m*a^2*omega^2-m*a*omega*l*(diff(phi(t), t))*sin(omega*t-phi(t))+(1/2)*m*l^2*(diff(phi(t), t))^2

(2.3.10)

For the gravitational energy, expressed in terms of the parametric equations of the point of mass m, we have

>	U := eval(U, parametric_equations);

(2.3.11)

So the requested Lagrangian is

>	L := combine(L, trig);;

(1/2)*m*a^2*omega^2-m*a*omega*l*(diff(phi(t), t))*sin(omega*t-phi(t))+(1/2)*m*l^2*(diff(phi(t), t))^2-m*g*a*sin(omega*t)+m*g*l*cos(phi(t))

(2.3.12)

Taking into account that the Lagrangian of a system is defined up to a total derivative with respect to t, we can eliminate the two terms that can be rewritten as total derivatives; these are (1/2)*m*a^2*omega^2 and so

>	select(has, L, [omega^2, sin(omega*t)]);

(1/2)*m*a^2*omega^2-m*g*a*sin(omega*t)

(2.3.13)

>	L := L - (2.3.13);

-m*a*omega*l*(diff(phi(t), t))*sin(omega*t-phi(t))+(1/2)*m*l^2*(diff(phi(t), t))^2+m*g*l*cos(phi(t))

(2.3.14)

__________________________________________________________

b) The steps are the same as in part a:

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(2.3.15)

>	L := T - U;

(2.3.16)

>	U := -mgy;

(2.3.17)

>	T := 1/2m(v_ . v_);

(2.3.18)

>	r_ := x_i + y_j;

(2.3.19)

Now, regarding part a), the only change is in the expression of the y coordinate, which for this part b) is:

>	y = 1 + l*cos(phi(t));

(2.3.20)

So the parametric equations in this case are

>	parametric_equations := [x = acos(omegat) + l* sin(phi(t)), (2.3.20)];

(2.3.21)

>	r_ := eval(r_, parametric_equations);

(2.3.22)

>	v_ := diff(r_, t);

(2.3.23)

For the gravitational energy, expressed in terms of the parametric equations of the point of mass m, we have

>	U := eval(U, parametric_equations);

(2.3.24)

So the requested Lagrangian is

(1/2)*m*((-a*sin(omega*t)*omega+l*cos(phi(t))*(diff(phi(t), t)))^2+l^2*sin(phi(t))^2*(diff(phi(t), t))^2)+m*g*(1+l*cos(phi(t)))

(2.3.25)

>	L := combine(L, trig);

(1/4)*m*a^2*omega^2-(1/4)*m*a^2*omega^2*cos(2*omega*t)-(1/2)*m*a*omega*l*(diff(phi(t), t))*sin(omega*t+phi(t))-(1/2)*m*a*omega*l*(diff(phi(t), t))*sin(omega*t-phi(t))+(1/2)*m*l^2*(diff(phi(t), t))^2+m*g+m*g*l*cos(phi(t))

(2.3.26)

The terms in L that can be expressed as total derivatives can be discarded, so

>	select(has, L, [omega^2, cos(2omegat)]);

(1/4)*m*a^2*omega^2-(1/4)*m*a^2*omega^2*cos(2*omega*t)

(2.3.27)

So the Lagrangian is

>	L := L - (2.3.27);

-(1/2)*m*a*omega*l*(diff(phi(t), t))*sin(omega*t+phi(t))-(1/2)*m*a*omega*l*(diff(phi(t), t))*sin(omega*t-phi(t))+(1/2)*m*l^2*(diff(phi(t), t))^2+m*g+m*g*l*cos(phi(t))

(2.3.28)

Inertia tensor for a triatomic molecule

Problem

Determine the Inertia tensor corresponding to a triatomic molecule that has the form of an isosceles triangle with two masses in the extremes of the base and mass in the third vertex. The distance between the two masses is equal to a, and the height of the triangle is equal to h.

Solution

>	restart; with(Physics, KroneckerDelta): with(Physics[Vectors]): Setup(mathematicalnotation = true);

(2.4.1)

The general formula for the Inertia tensor is given by

>	InertiaTensor := Sum(m[k] * (Norm(r_[k])^2 * kd_[i,j] - Component(r_[k], i)*Component(r_[k], j)), k=1..N);

Sum(m[k]*(Physics:-Vectors:-Norm(r_[k])^2*Physics:-KroneckerDelta[i, j]-Physics:-Vectors:-Component(r_[k], i)*Physics:-Vectors:-Component(r_[k], j)), k = 1 .. N)

(2.4.2)

where N is the number of particles, is the mass of each one, is its position in a reference system with the origin at the "center of mass", is the component in the i direction of the position vector associated to the k particle in that reference system, and is the Kronecker delta, part of Physics. Make this definition of the InertiaTensor be an indexing function for the InertiaTensor matrix.

>	IT := unapply(InertiaTensor, i,j);

proc (i, j) options operator, arrow; Sum(m[k]*(Physics:-Vectors:-Norm(r_[k])^2*Physics:-KroneckerDelta[i, j]-Physics:-Vectors:-Component(r_[k], i)*Physics:-Vectors:-Component(r_[k], j)), k = 1 .. N) end proc

(2.4.3)

For example, for a component in the diagonal, we have

IT(1,1);

Sum(m[k]*(Physics:-Vectors:-Norm(r_[k])^2-Physics:-Vectors:-Component(r_[k], 1)^2), k = 1 .. N)

(2.4.4)

and outside of the diagonal we have

IT(1,2);

Sum(-m[k]*Physics:-Vectors:-Component(r_[k], 1)*Physics:-Vectors:-Component(r_[k], 2), k = 1 .. N)

(2.4.5)

At this point, we can proceed to setting the particularities of this problem. The number of particles is 3.

N := 3;

(2.4.6)

Hence the matrix is

>	IT_Matrix := Matrix(3,3, IT);

IT_Matrix:=[[[(∑)m[k] (||(r)[k]||^2-(((r)[k])[1])^2),(∑)-m[k] ((r)[k])[1] ((r)[k])[2],(∑)-m[k] ((r)[k])[1] ((r)[k])[3]],[(∑)-m[k] ((r)[k])[1] ((r)[k])[2],(∑)m[k] (||(r)[k]||^2-(((r)[k])[2])^2),(∑)-m[k] ((r)[k])[2] ((r)[k])[3]],[(∑)-m[k] ((r)[k])[1] ((r)[k])[3],(∑)-m[k] ((r)[k])[2] ((r)[k])[3],(∑)m[k] (||(r)[k]||^2-(((r)[k])[3])^2)]]]

(2.4.7)

Two of the masses are equal

>	m[3] := m[1];

(2.4.8)

We now choose any system of reference (not at the center of mass) where we are going to project the position vectors of each atom and the center of mass . The vectors entering the definition of the Inertia tensor (2.4.2) are related to and by

>	position := r_[k] = R_[k] - R_[CM];

(2.4.9)

For the , we choose a system of reference with the origin at the middle of the segment connecting the two atoms of mass equal to . Using Cartesian coordinates, we take the x axis along this segment and the z axis passing through the third atom, of mass . So in this referential, the positions of the three atoms are

>	R_[1] := -a/2*_i; #to the left of the origin

(2.4.10)

>	R_[2] := h*_k; #along the z direction

(2.4.11)

>	R_[3] := a/2*_i; #to the right of the origin

(2.4.12)

Compute the position of the "center of mass." By definition, this is

>	R_[CM] := Sum(m[k]*R_[k], k=1..N)/Sum(m[k], k=1..N);

(Sum(m[k]*R_[k], k = 1 .. 3))/(Sum(m[k], k = 1 .. 3))

(2.4.13)

Evaluating these sums, we have the value of :

>	R_[CM] := value(R_[CM]);

m[2]*h*_k/(2*m[1]+m[2])

(2.4.14)

So these are the positions of the three particles viewed from the center of mass and expressed in terms of the known quantities , a and h:

>	seq(eval(position, k=j), j = 1..N);

r_[1] = -(1/2)*a*_i-m[2]*h*_k/(2*m[1]+m[2]), r_[2] = h*_k-m[2]*h*_k/(2*m[1]+m[2]), r_[3] = (1/2)*a*_i-m[2]*h*_k/(2*m[1]+m[2])

(2.4.15)

The answer to the problem posed, that is, the inertia tensor for this triatomic molecule, is obtained now by evaluating the abstract expression for the IT_Matrix at these values of the position vectors

>	IT_answer := simplify(eval(value(IT_Matrix), [(2.4.15)]));

$IT_answer := Matrix(3, 3, {(1, 1) = 2*m[2]*h^2*m[1]/(2*m[1]+m[2]), (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = m[1]*(2*m[1]*a^2+4*m[2]*h^2+a^2*m[2])/(4*m[1]+2*m[2]), (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = (1/2)*m[1]*a^2})$

(2.4.16)

Electrodynamics

Potential and electric field of a charged disk

Problem

Calculate the potential and the electric field of a disc of radius a, loaded with a surface density of charge (constant), in points of an axis perpendicular to the disk and passing through its center.

Solution

Given the potential , the general expression of the electric field is

In turn, the general expression for the potential of a distribution of charges over a surface is given by

Phi(r)=(∫)[Sigma]^(sigma((R)))/|(R) -( r)| &DifferentialD;s

where is the position vector of any point in space, is the position vector of any point of the disk, and is the surface element; the above expression is a surface integral with representing the integration domain.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(3.1.1)

The expression for the electric field as the gradient of the potential can be entered as

>	E_ := -%Gradient(Phi);

(3.1.2)

where in the above we are using the inert form of the Gradient command. The expression for in turn can be entered as a double integral in cylindrical coordinates (); the element of surface of a disk in these coordinates is , with varying from 0 to a and from 0 to .

>	Phi := Int(Int(sigmarho/Norm(r_ - R_), rho = 0..a), phi = 0..2Pi);

Int(Int(sigma*rho/Physics:-Vectors:-Norm(r_-R_), rho = 0 .. a), phi = 0 .. 2*Pi)

(3.1.3)

We choose the origin of the system of references to be the center of the disk and the z axis oriented perpendicular to the disk. In this system of references, the position vector of a point over the z axis is

>	r_ := z * _k;

(3.1.4)

and the position vector of a point of the disk is

>	R_ := rho*_rho;

(3.1.5)

The value of the potential for z > 0 can now be computed.

>	Phi := value(Phi) assuming z > 0;

(3.1.6)

From this we get

E_;

(3.1.7)

>	E_ := value((3.1.7));

2*Pi*sigma*((a^2+z^2)^(1/2)-z)*_k/(a^2+z^2)^(1/2)

(3.1.8)

Magnetic field of a rotating charged disk

Problem

A disk of radius a, uniformly charged with a surface density of charge rotates around its axis with a constant angular velocity , where is the cylindrical coordinate (the polar angle). Calculate the magnetic field on the axis of the disk.

Solution

The expression of the magnetic field due to a current of charges is

H(r)=(∫)[Sigma]^(J(R)*(r-R))/(c |r-R|^3) &DifferentialD;s

where is the position vector of any point in space, is the position vector of any point where the current exists, in this case a disk of radius a, and is the surface element. represents the integration domain and the above expression is a surface integral.

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(3.2.1)

The expression for can be entered as a double integral in cylindrical coordinates (); the element of the surface of a disk in these coordinates is , with varying from 0 to a and from 0 to .

>	H_ := Int(Int(J_&x (r_ - R_)/c/Norm(r_ - R_)^3rho, rho = 0..a), phi = 0..2Pi);

Int(Int(Physics:-Vectors:-`&x`(J_, r_-R_)*rho/(c*Physics:-Vectors:-Norm(r_-R_)^3), rho = 0 .. a), phi = 0 .. 2*Pi)

(3.2.2)

We choose the system of references as in the previous problem, with the origin in the center of the disk, and the z axis oriented perpendicular to the disk. So again the position vector of a point over the z axis is

>	r_ := z * _k;

(3.2.3)

and the position vector of a point of the disk is

>	R_ := rho*_rho;

(3.2.4)

By definition, the current at a point is equal to the value of the density of charge times the velocity of this charge; that is,

>	J_ := sigma * V_;

(3.2.5)

Finally, the velocity of a point of the disk can be computed as the derivative of with respect to t (the time), and in doing so we need to take into account that the unit vector and the disk is rotating.

This derivative of can be computed in two different ways. One way is to make explicit this dependence of on by changing the basis onto which is projected from the cylindrical to the Cartesian basis.

>	ChangeBasis(R_, 1);

(3.2.6)

Now make depend on t and differentiate.

>	subs(phi = phi(t), (3.2.6));

(3.2.7)

>	diff((3.2.7) ,t);

(3.2.8)

Introduce , and remove the explicit dependence of with respect to t to arrive at an expression for .

>	factor(subs([diff(phi(t),t) = omega, phi(t) = phi], (3.2.8)));

(3.2.9)

Alternatively (simpler), knowing that and so depends on the time only through , you can compute = For this purpose, use VectorDiff, to automatically take into account this dependence of on .

>	V_ := omega * VectorDiff(R_, phi);

(3.2.10)

At this point, we have all the quantities defined in the system of coordinates chosen and in terms of the constant angular velocity, and the radius of the disk, a. The expression of the magnetic field looks like

H_;

Int(Int((sigma*omega*rho*z*_rho+sigma*omega*rho^2*_k)*rho/(c*(rho^2+z^2)^(3/2)), rho = 0 .. a), phi = 0 .. 2*Pi)

(3.2.11)

However, to perform the integrals, we still need to express as a function of , one of the integration variables. For that purpose, it suffices to change the vectors involved ( and ) to the Cartesian basis.

>	R_ := ChangeBasis(R_, 1);

(3.2.12)

>	V_ := ChangeBasis(V_, 1);

(3.2.13)

With this change, looks like

H_;

Int(Int((sigma*rho*cos(phi)*omega*z*_i+sigma*rho*sin(phi)*omega*z*_j+sigma*omega*rho^2*_k)*rho/(c*(rho^2*cos(phi)^2+rho^2*sin(phi)^2+z^2)^(3/2)), rho = 0 .. a), phi = 0 .. 2*Pi)

(3.2.14)

and so the integrals can be performed, leading to the desired value of the magnetic field on the axis of the rotating disk.

>	H_ := value(H_) assuming a > 0, z > 0;

-2*sigma*omega*_k*Pi*(2*z*(a^2+z^2)^(1/2)-2*z^2-a^2)/((a^2+z^2)^(1/2)*c)

(3.2.15)

Magnetic potential and field of an oscillating distribution of currents

Problem

A distribution of currents varies according to

for some constant vectors and , where represents a point in space. Calculate the corresponding magnetic potential and magnetic field .

Solution

Maxwell equations are gauge invariant, from where the potential is defined up to the 4-divergence of a 4-vector; we can always explore that freedom to "gauge" according to what may be convenient. For static problems, this means we can freely choose the value of ; we choose here the Coulomb gauge, that is . In this gauge, the expression of the magnetic potential and magnetic vector due to a current of charges are as follows:

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(3.3.1)

>	eq[1] := %Laplacian(A_) = -4PiJ_/c;

(3.3.2)

>	eq[2] := H_ = %Curl(A_);

(3.3.3)

where in the above we are using the inert versions of the Laplacian and Curl commands; the computations will be performed by using the value command afterwards, when convenient.

Given , to compute we need to invert the Laplacian of . Knowing that the Laplacian is a scalar operator, and are either parallel or antiparallel. In this example the solving approach can also be simplified by the fact that is a pure oscillating function of the position , the second derivative of the cosine function is proportional to itself, and the Laplacian involves second order derivatives with respect to the position. So we can directly search for proportional to ; that is,

>	A_ := alpha*J_;

(3.3.4)

where

>	J_:= C_*cos(k_ . r_);

(3.3.5)

Orienting the z axis of the reference system in the direction of the constant vector , we have

>	C_ := abs(C)*_k;

(3.3.6)

where represents the module of . The vectors and are generic.

>	r_ := x_i + y_j + z*_k;

(3.3.7)

>	k_ := k[1]_i + k[2]_j + k[3]*_k;;

(3.3.8)

At this point, all of the vectors of the problem are stated and we need to compute the value of the constant of proportionality entering

eq[1];

%Laplacian(alpha*abs(C)*_k*cos(k[1]*x+k[2]*y+k[3]*z)) = -4*Pi*abs(C)*_k*cos(k[1]*x+k[2]*y+k[3]*z)/c

(3.3.9)

>	value((3.3.9));

-alpha*abs(C)*_k*cos(k[1]*x+k[2]*y+k[3]*z)*(k[1]^2+k[2]^2+k[3]^2) = -4*Pi*abs(C)*_k*cos(k[1]*x+k[2]*y+k[3]*z)/c

(3.3.10)

Note that is a common factor in the left- and right-hand sides of this equation, so instead of taking coefficients of , we can proceed directly with solving for .

>	isolate((3.3.10), alpha);

alpha = 4*Pi/(c*(k[1]^2+k[2]^2+k[3]^2))

(3.3.11)

>	simplify((3.3.11));

alpha = 4*Pi/(c*(k[1]^2+k[2]^2+k[3]^2))

(3.3.12)

Let's assign this value of so that all expressions start taking it into account automatically.

>	assign((3.3.12));

So is given by

A_;

4*Pi*abs(C)*_k*cos(k[1]*x+k[2]*y+k[3]*z)/(c*(k[1]^2+k[2]^2+k[3]^2))

(3.3.13)

and is given by

eq[2];

H_ = %Curl(4*Pi*abs(C)*_k*cos(k[1]*x+k[2]*y+k[3]*z)/(c*(k[1]^2+k[2]^2+k[3]^2)))

(3.3.14)

>	value((3.3.14));

H_ = -4*Pi*abs(C)*sin(k[1]*x+k[2]*y+k[3]*z)*k[2]*_i/(c*(k[1]^2+k[2]^2+k[3]^2))+4*Pi*abs(C)*sin(k[1]*x+k[2]*y+k[3]*z)*k[1]*_j/(c*(k[1]^2+k[2]^2+k[3]^2))

(3.3.15)

Electric field around a conductive neutral sphere in an exterior homogeneous field

Problem

Determine the electric field around a conductive neutral sphere of radius a submersed in an exterior homogeneous electric field .

Solution

>	restart; with(Physics[Vectors]): Setup(mathematicalnotation = true);

(3.4.1)

The solution to the problem, generally speaking, is that the (total) electric field can be computed as the negative gradient of the (total) potential,

>	E_[T] := - %Gradient(Phi[T]);

(3.4.2)

where the potential is the sum of the potentials associated with the electric field induced in the neutral conductive sphere and the external homogeneous electric field, respectively.

>	Phi[T] := Phi[s] + Phi[0];

(3.4.3)

The expression for the potential for which is a constant vector can be obtained by reversing this formula; that is, by calculating a line integral with a constant integrand, arriving at:

>	Phi[0] := E_[0] . r_;

(3.4.4)

where is the position vector of any point outside the sphere. An expression for the potential of the sphere, , can be derived from the observation that, in the presence of a homogeneous external field, in the neutral conductive sphere, positive and negative charges are polarized (total charge remains 0) so that a dipole is induced. The potential of a dipole in an external field, in turn, is proportional (by a factor to be determined) to the scalar product of the external field with the potential of an electric "monopole." A multipole expansion of the potential shows that this result is independent of the geometry of the problem.

For the particular case of a sphere, choosing the origin of the system of references at the center of the sphere and aligning the z axis with the external field , we thus have, for :

>	Phi[s] := alpha * E_[0] . %Nabla(1/r);

(3.4.5)

That is,

>	Phi[s] := value(%);

-alpha*Physics:-Vectors:-`.`(E_[0], _r)/r^2

(3.4.6)

Note that the output above is directly in spherical coordinates and using the spherical basis of unit vectors; this is because the argument of is 1/r, where r is the spherical coordinate. In this orthonormal basis, the position vector is written as

>	r_ := r*_r;

(3.4.7)

The external homogeneous electric field is a constant vector along the z axis.

>	E_[0] := E[0] * _k;

(3.4.8)

In order to proceed, we need to project all of the quantities in the same orthonormal basis; that is, change the projection of from the Cartesian basis to the spherical basis.

>	E_[0] := ChangeBasis(E_[0], 3);

(3.4.9)

At this point all the relevant quantities have their values derived, so the total potential evaluates to

Phi[T];

-alpha*E[0]*cos(theta)/r^2+r*E[0]*cos(theta)

(3.4.10)

In the above we see that the problem has been mapped into determining the value of a single unknown, the constant of proportionality . This constant can be determined recalling that the potential on the surface of a neutral conducting sphere is always constant, and that the actual value of that constant is arbitrary (what we measure is , not ). Thus, choosing the value of over the surface of the sphere to be zero, and recalling that the radius of the sphere is equal to a, we have

>	Eval('Phi[T]', r = a) = 0;

Eval(Phi[T], r = a) = 0

(3.4.11)

Solving now for

>	alpha := solve(value((3.4.11)), alpha);

(3.4.12)

So is equal to

Phi[T];

-a^3*E[0]*cos(theta)/r^2+r*E[0]*cos(theta)

(3.4.13)

Hence the total electric field is

>	E_[T] := value(E_[T]);

-E[0]*cos(theta)*(2*a^3+r^3)*_r/r^3-E[0]*sin(theta)*(a^3-r^3)*_theta/r^3

(3.4.14)

Note in the above the use of value to evaluate the inert %Gradient.

Tensors in Special and General Relativity

>	restart; with(Physics);

(4.1)

>	Setup(mathematicalnotation = true);

(4.2)

To define a system of references (coordinates system) and the related spacetime vector you can use Setup or Coordinates;

>	Coordinates(X);

$`Default differentiation variables for d_, D_ and dAlembertian are: `*{X = (x1, x2, x3, x4)}$

(4.3)

You can now use X to represent function dependency, as in , equivalent to writing F(), and also as the spacetime vector (see also SpaceTimeVector). To indicate that an index is contravariant, prefix it with ~ (tilde)

>	X[~mu], X[mu];

(4.4)

The label X can also be used to select each of the components of the spacetime vector

>	X[1], X[3];

(4.5)

and is always mapped into

x0;

(4.6)

You can also set the coordinates to be any sequence of four names. Three predefined sets are 'cartesian', 'cylindrical' and 'spherical'

>	Coordinates(X = cartesian);

$`Default differentiation variables for d_, D_ and dAlembertian are: `*{X = (x, y, z, t)}$

(4.7)

In all cases you will still be abel to refer to each coordinate using the generic symbols x1, x2, ...

x1, x2;

(4.8)

x0 = x4;

(4.9)

You can set many systems of coordinates simultaneously - although only one is considered the 'default differentiation variables' for the d_ dAlembertian and D_ operators (and that is the system of references were are all the general relativity tensors are defined).

>	Coordinates(Y);

(4.10)

y0 = y4;

(4.11)

>	Setup(diff);

(4.12)

To change the differentiation variables enter for instance Setup(diff = Y);

When you load Physics, some package's commands that are automatically set as spacetime tensors, such as the metric g_[mu, nu] or the differential operator d_[mu]. Every other symbol that you want to be considered a spacetime tensor must be defined as such using the Define command, and nothing else will be considered a spacetime tensor during computations. To see the tensors predefined you can call Define without arguments

Define();

${Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(4.13)

This set includes the special and general relativity tensors as well as the Pauli and Dirac matrices.

The default dimension of spacetime is

>	Setup(dimension);

(4.14)

To change the value of the spacetime dimension use Setup(dimension = N) where N is any positive integer greater than 1.

To see the values of the components of the metric you can enter the spactime metric g_ without indices; when you load the package it is of Minkowski type

g_[];

g[mu, nu] = (Matrix(4, 4, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}))

(4.15)

The spacetime indexed differentiation operator

>	d_[mu](F(X));

(4.16)

This operator is also a representation for the differential of a function; for that purpose enter it without indics

d_(F(X));

(4.17)

Using d_ you can express the differential of any coordinate defined using Coordinates

d_(x0);

(4.18)

d_(y0);

(4.19)

The differential of everything else is zero

d_(a);

(4.20)

The diff and d_ operators are interconnected

>	diff(F(X), X[~mu]);

(4.21)

>	diff(F(X), X[mu]);

(4.22)

By applying it two times you get the dAlembertian

>	d_[mu](d_[mu](F(X)));

(4.23)

You can enter any pair of contracted indices one covariant and one contravariant, not being relevant which one is of which kind, so you can also enter both covariant or both contravariant and the system will automatically rewrite them as one covariant and one contravariant.

>	g_[mu, nu]^2;

(4.24)

Define two tensors for exploration purposes

>	Define(A, B);

${A, B, Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], Physics:-SpaceTimeVector[mu](X)}$

(4.25)

Note the automatic rewriting of repeated indices as one covariant and one contravariant

>	g_[alpha, mu] * A[mu] * g_[alpha, nu] * B[nu, sigma, sigma];

(4.26)

When a pair of contracted indices is already one covariant and one contravariant, it is left unchanged. Although it is not recommented, you can switch off the Physics evaluator that performs this automatic rewriting of contracted indices by entering Setup(usephysicsevaluator = false).

To determine the repeated and free indices in an expression use Check

>	Check((4.26), all);

(4.27)

The main simplifier in the Physics package is Simplify. The simplification of tensor indices takes into account the sum rule for repeated indices

>	Simplify((4.26));

(4.28)

Besides its use in quantum mechanics, the Physics `.` operator is also a handy shortcut for simplifying contracted indices - you use it replacing the product operator *. This is useful when you want the simplification to happen only in some places. Recalling that multiplication is left associative, you may need to put parenthesis to get what you want. For example, replace in only the third `*` by `.`

>	g_[alpha, mu] * A[mu] * (g_[alpha, nu] . B[nu, sigma, sigma]);

(4.29)

Set the spacetime metric to be the Schwarzschild metric in spherical coordinates: you can now do that in various manners, the simplest of which is to pass part of the identifying keyword directly to the metric tensor g_ and in this way you also set the corresponding coordinates all in one step (setting thus automatically redefines the coordinates label )

g_[sc];

$g[mu, nu] = (Matrix(4, 4, {(1, 1) = r/(-r+2*m), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -r^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -r^2*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (r-2*m)/r}))$

(4.30)

The metric can also be set passing any line element to Setup, as in Setup(metric = line_element). The line element is expected to be quadratic in the differentials of the coordinates that should have been defined using Coordinates or Setup, or can be defined simultaneously for instance as in Setup(coordinates = (Z = [...]), metric = line_element).

As a rule for all the relativity tensors, you can see their components for different convariant or contravariant indices by entering the tensor with the indices and the keyword matrix at the end, as in

>	g_[mu, ~nu, matrix];

g[mu]^(nu)=[[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]]

(4.31)

>	g_[~mu, ~nu, matrix];

g[]^(mu,nu)=[[[(-r+2 m)/(r),0,0,0],[0,-1/(r^2),0,0],[0,0,-1/(r^2 (sin(theta))^2),0],[0,0,0,-(r)/(-r+2 m)]]]

(4.32)

You can also check the value of any tensor component by giving numerical values to the indices; prefix them with ~ to indicate that they are contravariant

g_[1,1];

(4.33)

>	g_[~1, ~1];

(4.34)

A coordinate transformation to a conformal euclidean form, i.e. with line element proportional to

>	TR := r = (1+2m/(4rho))^2*rho;

r = (1+(1/2)*m/rho)^2*rho

(4.35)

To transform coordinates in a tensor expression - this is the standard tensor transformation related to its free indices - use TransformCoordinates

>	TransformCoordinates(TR, g_[mu, nu], [rho, theta, phi, t]);

$Matrix(4, 4, {(1, 1) = -(1/16)*(2*rho+m)^4/rho^4, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = -(1/16)*(2*rho+m)^4/rho^2, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -(1/16)*(2*rho+m)^4*sin(theta)^2/rho^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = (m-2*rho)^2/(2*rho+m)^2})$

(4.36)

Check the line element using the related option; the differentials of the coordinates are expressed using d_

>	TransformCoordinates(TR, g_[mu, nu], [rho, theta, phi, t], output = line_element);

-(1/16)*(2*rho+m)^4*(%d_(rho)^2+(%d_(theta)^2+sin(theta)^2*%d_(phi)^2)*rho^2)/rho^4+(m-2*rho)^2*%d_(t)^2/(2*rho+m)^2

(4.37)

TransformCoordinates does not set the metric to the result of the transformation; you can do that by passing its output to Setup, in this example that would be Setup(metric = .

The general relativity tensors automatically update their value when you set the metric. This is for instance contraction of all the indices of the Riemann tensor for the (current) Schwarzschild metric

>	Riemann[alpha,beta,mu,nu] . Riemann[alpha,beta,mu,nu];

48*m^2/r^6

(4.38)

When computing algebraically with the general relativity tensors of the Physics package, their indices are automatically sorted according to their symmetry properties to facilitate zero recognition

>	Riemann[mu, mu, beta, alpha];

(4.39)

>	Riemann[nu, mu, beta, alpha] - Riemann[alpha, beta, nu, mu];

(4.40)

Besides algebraic computations, it is possible to check the value of all or a subset of the tensor's components in a visual way. These are the values of and

>	Christoffel[1, 2, 2];

(4.41)

>	Christoffel[~1, 2, 2];

(4.42)

The 2 x 2 matrix corresponding to the Christoffel symbols with the first index contravariant is displayed when you pass the keyword matrix as an extra index

>	Christoffel[~1, alpha, mu, matrix];

GAMMA[alpha,mu]^(1)=[[[m/(r (-r+2 m)),0,0,0],[0,-r+2 m,0,0],[0,0,(-r+2 m) (sin(theta))^2,0],[0,0,0,(-2 m^2+m r)/(r^3)]]]

(4.43)

Compare with the matrix for (so the first index is covariant and then not prefixed with ~)

>	Christoffel[1, alpha, mu, matrix];

$GAMMA[1, alpha, mu] = (Matrix(4, 4, {(1, 1) = m/(-r+2*m)^2, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = r, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = r*sin(theta)^2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = -m/r^2}))$

(4.44)

The nonzero components of

>	Christoffel[~mu, alpha, beta, nonzero];

$Physics:-Christoffel[`~mu`, alpha, beta] = {(1, 1, 1) = m/(r*(-r+2*m)), (1, 2, 2) = -r+2*m, (1, 3, 3) = (-r+2*m)*sin(theta)^2, (1, 4, 4) = (-2*m^2+m*r)/r^3, (2, 1, 2) = 1/r, (2, 2, 1) = 1/r, (2, 3, 3) = -sin(theta)*cos(theta), (3, 1, 3) = 1/r, (3, 2, 3) = cos(theta)/sin(theta), (3, 3, 1) = 1/r, (3, 3, 2) = cos(theta)/sin(theta), (4, 1, 4) = -m/(r*(-r+2*m)), (4, 4, 1) = -m/(r*(-r+2*m))}$

(4.45)

For the all covariant you can pass all the indices covariant, or simpler: pass the keyword nonzero alone

>	Christoffel[nonzero];

$Physics:-Christoffel[mu, nu, alpha] = {(1, 1, 1) = m/(-r+2*m)^2, (1, 2, 2) = r, (1, 3, 3) = r*sin(theta)^2, (1, 4, 4) = -m/r^2, (2, 1, 2) = -r, (2, 2, 1) = -r, (2, 3, 3) = r^2*sin(theta)*cos(theta), (3, 1, 3) = -r*sin(theta)^2, (3, 2, 3) = -r^2*sin(theta)*cos(theta), (3, 3, 1) = -r*sin(theta)^2, (3, 3, 2) = -r^2*sin(theta)*cos(theta), (4, 1, 4) = m/r^2, (4, 4, 1) = m/r^2}$

(4.46)

With the Schwarzschild metric, the Ricci tensor has all its components equal to zero. To see the matrix form, either pass the keyword matrix, or simpler: pass no indices

Ricci[];

R[mu, nu] = (Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}))

(4.47)

As is the case in the literature, both the Ricci and the Riemann tensors are displayed with a capital R, but the correct value for each case is retrieved through copy and paste.

The Ricci tensor is defined in terms of the Riemann tensor as You can see how the four Riemann matrices add to form the Ricci matrix; for this purpose give values to the 1st and 3rd indices and pass the extra keyword matrix. This is the first matrix

>	Riemann[~1, beta, 1, nu, matrix];

R[beta,1,nu]^(1)=[[[0,0,0,0],[0,-m/r,0,0],[0,0,-(m (sin(theta))^2)/r,0],[0,0,0,(4 m^2-2 m r)/(r^4)]]]

(4.48)

To add the four of them use +

>	Riemann[~1, beta, 1, nu, matrix] + Riemann[~2, beta, 2, nu, matrix] + Riemann[~3, beta, 3, nu, matrix] + Riemann[~4, beta, 4, nu, matrix];

R[beta,1,nu]^(1)+R[beta,2,nu]^(2)+R[beta,3,nu]^(3)+R[beta,4,nu]^(4)=[[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,(4 m^2-2 m r)/(r^4)+(2 (-2 m^2+m r))/(r^4)]]]

(4.49)

>	simplify((4.49));

R[beta,1,nu]^(1)+R[beta,2,nu]^(2)+R[beta,3,nu]^(3)+R[beta,4,nu]^(4)=[[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]]

(4.50)

The Riemann scalars and (see the corresponding help page)

>	Riemann[scalars];

S[1] = m^2/r^6, S[2] = -m^3/r^9

(4.51)

The covariant derivative of

>	D_[mu](A[nu](X));

(4.52)

>	expand((4.52));

(4.53)

The divergence of

>	D_[mu](A[mu](X));

(4.54)

>	expand((4.54));

((1/2)*A[`~mu`](X)*(4*r^3*sin(theta)^2*Physics:-KroneckerDelta[mu, `~1`]+2*r^4*sin(theta)*cos(theta)*Physics:-KroneckerDelta[mu, `~2`])/(r^4*sin(theta)^2)^(1/2)+(r^4*sin(theta)^2)^(1/2)*Physics:-d_[mu](A[`~mu`](X), [X]))/(r^4*sin(theta)^2)^(1/2)

(4.55)

>	simplify((4.55)) assuming positive;

(2*A[`~mu`](X)*sin(theta)*Physics:-KroneckerDelta[mu, `~1`]+A[`~mu`](X)*r*cos(theta)*Physics:-KroneckerDelta[mu, `~2`]+sin(theta)*Physics:-d_[mu](A[`~mu`](X), [X])*r)/(r*sin(theta))

(4.56)

You can simplify the remaining contracted indices using Simplify.

The covariant derivative for a tensor with mixed indices; in some cases expand will expand not just the covariant derivative; you can be more selective using convert; to avoid redudant display of the functionality of use PDEtools:-declare

>	PDEtools:-declare(B(X));

(4.57)

>	D_[mu](B[~alpha, beta, ~nu](X));

(4.58)

>	convert((4.58), d_);

(4.59)

This expression can be rewritten in different ways. For example, rewrite Christoffel in terms of g_ and its derivatives, expand, then revert back the rewriting

>	convert((4.59), g_);

(4.60)

>	expand((4.60));

(4.61)

>	convert((4.61), Christoffel);

(4.62)

>	Simplify((4.62));

(4.63)

which is the same as (4.59).

For more details see Spacetime, tensor simplification and default settings, and the help pages for Coordinates, d_, D_, Einstein, g_, KroneckerDelta (abbreviation: kd_), LeviCivita (abbreviation: ep_), Ricci, Riemann, SpaceTimeVector and Weyl.

Classical field theory

Maxwell equations departing from the Action for electrodynamics

Maxwell equations result from taking the functional derivative of the Action.

Let X and Y represent two spacetime points

>	restart; with(Physics):

>	Coordinates(X, Y);

$`Default differentiation variables for d_ and dAlembertian are: `*{X = (x1, x2, x3, x4)}$

(5.1.1)

Define now the 4-D electromagnetic pontential

>	Setup(tensors = A[mu](X));

$[tensors = {Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, beta, mu, nu], A[mu](X), Physics:-SpaceTimeVector[mu](X)}]$

(5.1.2)

From herein avoid displaying the functionality whenever it is X

>	PDEtools:-declare(A(X));

(5.1.3)

Also to avoid having to repeat saying who are the differentiation variables of the operators and just set them

>	Setup(var = X);

$`Default differentiation variables for d_ and dAlembertian are: `*{X = (x1, x2, x3, x4)}$

(5.1.4)

Query about the dimension and signature ...

>	Setup(dimension, signature);

(5.1.5)

Everything is set. Let's define the electromagnetic field tensor in terms of the derivatives of

>	F[mu, nu] := d_[mu](A[nu](X)) - d_[nu](A[mu](X));

(5.1.6)

Define now the Action

>	S := Intc(F[mu,nu]^2, X);

Int(Int(Int(Int((Physics:-d_[mu](A[nu](X), [X])-Physics:-d_[nu](A[mu](X), [X]))^2, x1 = -infinity .. infinity), x2 = -infinity .. infinity), x3 = -infinity .. infinity), x4 = -infinity .. infinity)

(5.1.7)

Take now the functional derivative of S (use ' to delay the computation, to see what is that will be computed)

>	'Fundiff'(S , A[rho](Y));

Physics:-Fundiff(Int(Int(Int(Int((Physics:-d_[mu](A[nu](X), [X])-Physics:-d_[nu](A[mu](X), [X]))^2, x1 = -infinity .. infinity), x2 = -infinity .. infinity), x3 = -infinity .. infinity), x4 = -infinity .. infinity), A[rho](Y))

(5.1.8)

Substitute also Y = X to take advantage of the automatic compact display when A[mu], and the d'Alembertian depend on X instead of Y

>	subs(Y = X, (5.1.8));

(5.1.9)

To simplify the contracted spacetime indices, use the Simplify command, resulting in the Maxwell equations in their familiar 4-D tensorial form

>	Simplify((5.1.9));

(5.1.10)

The field equations for the model

The Lagrangean of the model, the corresponding Action, and the field equations:

>	PDEtools:-declare(Phi(X));

(5.2.1)

>	L := 1/2 * d_[mu](Phi(X)) * d_[mu](Phi(X)) - m^2/2 * Phi(X)^2 + lambda/4*Phi(X)^4;

(1/2)*Physics:-d_[mu](Phi(X), [X])*Physics:-d_[`~mu`](Phi(X), [X])-(1/2)*m^2*Phi(X)^2+(1/4)*lambda*Phi(X)^4

(5.2.2)

>	S := Intc(L,X);

Int(Int(Int(Int((1/2)*Physics:-d_[mu](Phi(X), [X])*Physics:-d_[`~mu`](Phi(X), [X])-(1/2)*m^2*Phi(X)^2+(1/4)*lambda*Phi(X)^4, x1 = -infinity .. infinity), x2 = -infinity .. infinity), x3 = -infinity .. infinity), x4 = -infinity .. infinity)

(5.2.3)

>	Fundiff(S, Phi(Y));

(5.2.4)

For a more compact display, take advantage that X has been defined as the default differentiation variables for the d'Alembertian:

>	subs(Y = X, (5.2.4));

(5.2.5)

Quantum Mechanics

Mathematical tools and their representation in Dirac notation

This section aims to summarize the mathematical objects entering formulations in quantum mechanics and to show how these objects are entered and represented in the Maple worksheet.

>	restart; with(Physics): Setup(mathematicalnotation = true);

(6.1.1)

Kets and Bras

The quantum state of a system, belonging to a space of quantum states, is represented in Dirac notation by a Ket state-vector.

Ket(u);

(6.1.1.1)

The above is the quantum analog of a non-projected vector of a 3-D Euclidean space. The norm of a generic Ket is not predefined and can be indicated by setting a bracket rule, as shown below. Every Ket can be projected onto a basis of state-vectors of the space to which the Ket belongs. The Kets conforming a basis (analogous to 3-D unit vectors) are distinguished from a generic Ket by the fact that they have one or many quantum numbers, as in

>	Ket(v, n);

(6.1.1.2)

Kets having quantum numbers are always assumed to belong to a basis of quantum states, are orthogonal to each other, have norm equal to 1, and are distinguished from each other by the values of these quantum numbers. For example, an orthonormal basis of a two dimensional space of quantum states is

>	Ket(v,0), Ket(v,1);

(6.1.1.3)

There are no restrictions on the number of quantum numbers that a Ket can have. This is a Ket belonging to a basis of a space that depends on four quantum numbers.

>	Ket(v, j,k,m,n);

(6.1.1.4)

You can associate a space of states with each quantum number, so a Ket with many quantum numbers represents a state in a space constructed as a tensor product of spaces.

There is a Bra associated with each Ket, obtained from the Ket by performing the Hermitian conjugate, or Dagger, operation.

>	Dagger(%);

(6.1.1.5)

>	Dagger(%);

(6.1.1.6)

You can enter Bras directly by using the Bra function.

>	Bra(v, j,k,m,n);

(6.1.1.7)

The space of Bras of a system is the dual of the space of Kets of that system.

Discrete and continuous basis of states

Kets belong to either a discrete or a continuous spaces of quantum states. A discrete space of states is one where the quantum numbers of its Kets vary discretely. These Kets thus belong to discrete bases of quantum states. A continuous space is one where the quantum numbers vary continuously and its Kets belong to continuous bases of quantum states. Unless explicitly stated otherwise, Kets are assumed to belong to discrete space of states.

You can indicate that a label R identifies a continuous space of states by using the Setup command.

>	Setup(continuous = R);

(6.1.2.1)

Kets of a continuous space of states can also have any number of quantum numbers.

>	Ket(R, x);

(6.1.2.2)

>	Ket(R, x,y,z);

(6.1.2.3)

>	Dagger(%);

(6.1.2.4)

Scalar product and orthonormalization relation

The scalar product is defined between a Bra and a Ket, in that order, and can be performed by using the dot operator `.` of the Physics package, or by using the Bracket function; both represent the same object.

>	Bra(u) . Ket(v);

(6.1.3.1)

>	Bracket(Bra(u), Ket(v));

(6.1.3.2)

Note that when the scalar product is just represented, not actually computed, as in the above, the result is always expressed in terms of the Bracket function. A shortcut notation for entering the scalar product using the Bracket function is

>	Bracket(u, v);

(6.1.3.3)

Under the Dagger operation, so that the Bracket becomes

>	Dagger(%);

(6.1.3.4)

For the Bracket, the same happens under conjugation:

>	conjugate(%);

(6.1.3.5)

Two generic Kets such as the above may or not belong to the same space. To make practical use of Kets, depending on the problem you may want to set a bracket rule, stating the value of the Bracket between them, by using the Setup command.

>	Setup(% = f(u,v));

(6.1.3.6)

After that, both the dot operator of the Physics package and the Bracket function know how to perform their scalar product.

>	Bracket(u, v);

(6.1.3.7)

Kets belonging to a discrete basis of states satisfy an orthonormalization relation involving the KroneckerDelta symbol, and when the quantum numbers are present, you do not need to specify a bracket rule.

>	%Bracket(Bra(u, n) , Ket(u, m)) = Bra(u, n) . Ket(u, m);

(6.1.3.8)

Note in the above the inert form %Bracket. It is sometimes useful to represent mathematical operations without having them actually performed. For that purpose, use any command with the name prefixed by the symbol %. To have the inert operation performed, use the value command.

value(%);

(6.1.3.9)

Evaluate the the output two operations above at .

>	eval(%%, m = n);

(6.1.3.10)

The Bracket of state-vectors depending on many quantum numbers results in products of KroneckerDelta symbols. The shortcut notation of the Bracket function also works in the presence of quantum numbers (tip: to avoid typographical mistakes, it is practical to group the objects visually, by leaving spaces or not after the commas).

>	Bracket(u, i,j,k, u, n,m,l);

(6.1.3.11)

Kets of a continuous basis of states satisfy an orthonormalization relation involving the Dirac function (recall that R has been set as a label of a continuous space, by using the Setup command, above).

>	Bracket(R, x, R, y);

(6.1.3.12)

>	%Bracket(R, x,y,z, R, a,b,c): % = value(%);

(6.1.3.13)

In the above, the 3-dimensional Dirac function can be expanded using expand

>	expand((6.1.3.13));

(6.1.3.14)

Closure relation, Projectors

Every Ket of a space of states can be expanded into a basis of that space. The operator that performs the expansion is called a Projector. To construct these projectors, information about the basis dimension is necessary. You can indicate this dimension directly to the Projector command, or set it by using Setup. The information available at this point is

>	Setup(basisdim, cont);

(6.1.4.1)

By default, continuous bases are assumed to range from to , so this information on R is enough to compute its Projector.

>	P[R] := Projector(Ket(R, x,y,z));

Int(Int(Int(Physics:-`*`(Physics:-Ket(R, x, y, z), Physics:-Bra(R, x, y, z)), x = -infinity .. infinity), y = -infinity .. infinity), z = -infinity .. infinity)

(6.1.4.2)

This expression for the projector is also called the closure relation; together with the orthonormalization relation , it tells that the set of Kets containing all the possible values of x, y, and z forms a basis, and so any has a unique expansion onto .

Note that the scalar product of with itself is equal to itself.

>	P[R] . P[R];

Int(Int(Int(Physics:-`*`(Physics:-Ket(R, x, y, z), Physics:-Bra(R, x, y, z)), x = -infinity .. infinity), y = -infinity .. infinity), z = -infinity .. infinity)

(6.1.4.3)

The following is the projector for a basis generically labeled u that has not been set to represent a continuous basis. By default, if nothing is known about the label of a basis, it is assumed to be related to a discrete space of states. The dimension of the basis can be indicated directly to the Projector command.

>	P[u] := Projector(Ket(u, n), dimension = N);

Sum(Physics:-`*`(Physics:-Ket(u, n), Physics:-Bra(u, n)), n = 0 .. N-1)

(6.1.4.4)

The information passed to Projector is automatically tracked by the system, so you do not need to give it again.

>	Setup(quantumbasisdimension, quantumcontinuousbasis);

(6.1.4.5)

To change this information, see Setup and its redo option. In order to compute scalar products of Kets belonging to a basis with other Kets of the same space, you can define a bracketrule .

>	%Bracket(Bra(R, x,y,z), Ket(psi)) = psi(x, y, z);

(6.1.4.6)

Setup(%);

(6.1.4.7)

Now Bracket and the `.` operator of the Physics package know how to compute a number of related operations.

>	Bracket(Bra(R, a,b,c), Ket(psi));

(6.1.4.8)

>	P[R] . Ket(psi);

Int(Int(Int(Physics:-Ket(R, x, y, z)*psi(x, y, z), x = -infinity .. infinity), y = -infinity .. infinity), z = -infinity .. infinity)

(6.1.4.9)

>	Bra(psi) . P[R];

Int(Int(Int(conjugate(psi(x, y, z))*Physics:-Bra(R, x, y, z), x = -infinity .. infinity), y = -infinity .. infinity), z = -infinity .. infinity)

(6.1.4.10)

>	Bra(psi) . P[R] . Ket(psi);

Int(Int(Int(abs(psi(x, y, z))^2, x = -infinity .. infinity), y = -infinity .. infinity), z = -infinity .. infinity)

(6.1.4.11)

This is a bracket rule for the scalar product of a state-vector of the discrete basis u and .

>	%Bracket(Bra(u, n), Ket(psi)) = psi(n);

(6.1.4.12)

Setup(%);

(6.1.4.13)

This rule for permits projecting onto the u basis, which is equivalent to inserting a projector between and . Note the use of delay evaluation quotes surrounding the Bracket; the operation is performed in the next line by using %.

>	' Bracket(psi, P[u], psi) ';

(6.1.4.14)

Sum(abs(psi(n))^2, n = 0 .. N-1)

(6.1.4.15)

A Ket can have differenty types of spaces associated with its quantum numbers. In the following example, Kets from a basis B have four quantum numbers, two of which, and , are associated with continuous spaces, and the dimension of the space associated with each quantum number is different.

>	Setup(quantumcontinuousbasis = {B[3], B[4]}, quantumbasisdimension = {B[1] = -1/2..1/2, B[2] = 0..N, B[3] = -a..a, B[4] = -infinity..infinity});

[quantumbasisdimension = {R = infinity, u = N, B[1] = -1/2 .. 1/2, B[2] = 0 .. N, B[3] = -a .. a, B[4] = -infinity .. infinity}, quantumcontinuousbasis = {R, B[3], B[4]}]

(6.1.4.16)

This is the projector onto the basis B:

>	Projector(Ket(B, n,m, x,y));

Int(Int(Sum(Sum(Physics:-`*`(Physics:-Ket(B, n-1/2, m, x, y), Physics:-Bra(B, n-1/2, m, x, y)), n = 0 .. 1), m = 0 .. N), x = -a .. a), y = -infinity .. infinity)

(6.1.4.17)

Quantum operators, eigenvectors, eigenvalues and commutators

To indicate to the system that a letter represents a quantum operator, use the Setup command; this sets B as a quantum operator.

>	Setup(op = B);

(6.1.5.1)

Note that after having entered > Setup(mathematicalnotation = true); the extended typesetting displays noncommutative objects in different colors; to change this color, see ?Setup. Because B is now a quantum operator, is an eigenvector of the four operators , and , with eigenvalues m, n, x, and y, respectively. For example,

>	B[2] . Ket(B, m,n,x,y);

(6.1.5.2)

Quantum operators can also be 3-D Euclidean vectors; for that purpose, you must load the Physics[Vectors] subpackage.

>	with(Vectors);

(6.1.5.3)

Set the vectors , , and as quantum operators (note the use of the option redo to erase previous definitions of quantum operators).

>	Setup(op = {L_, r_, p_, x,y,z, px, py, pz, Lx, Ly, Lz}, redo);

$[quantumoperators = {L_, Lx, Ly, Lz, p_, px, py, pz, r_, x, y, z}]$

(6.1.5.4)

Define as the angular momentum operator , and set commutation rules for the components of and .

>	L_ := r_ &x p_;

(6.1.5.5)

>	r_ := x_i + y_j + z*_k;

(6.1.5.6)

>	p_ := px_i + py_j + pz*_k;

(6.1.5.7)

Vector calculus with noncommutative components preserves and . This is achieved by symmetryzing and antisymmetrizing, respectively, each of these products. For example, below is the Component of along the x-axis.

_i . L_;

(6.1.5.8)

To enter the commutation rules between each component of and with each other, you can write these commutators and pass the whole set to Setup. When there are many, as in this case, it is more convenient to use a Matrix and an indexing function. Enter the core information as a procedure: C represents the Commutator of the Components of the vectors and .

>	C := (a_, i, b_, j) -> %Commutator(Component(a_, i), Component(b_, j));

(6.1.5.9)

So, given i and j from 1 to 3 identifying the components of and , an algebra can be set as set as follows.

>	algebra := (i, j) -> (

>	C(r_, i, p_, j) = I*KroneckerDelta[i,j],

>	C(r_, i, r_, j) = 0,

>	C(p_, i, p_, j) = 0);

(6.1.5.10)

Now all of the commutators between each component of and can be constructed with one call to Matrix.

>	Matrix(3, 3, algebra);

[[[[x,px][-]=&ImaginaryI;,[x,x][-]=0,[px,px][-]=0,[x,py][-]=0,[x,y][-]=0,[px,py][-]=0,[x,pz][-]=0,[x,z][-]=0,[px,pz][-]=0],[[y,px][-]=0,[y,x][-]=0,[py,px][-]=0,[y,py][-]=&ImaginaryI;,[y,y][-]=0,[py,py][-]=0,[y,pz][-]=0,[y,z][-]=0,[py,pz][-]=0],[[z,px][-]=0,[z,x][-]=0,[pz,px][-]=0,[z,py][-]=0,[z,y][-]=0,[pz,py][-]=0,[z,pz][-]=&ImaginaryI;,[z,z][-]=0,[pz,pz][-]=0]]]

(6.1.5.11)

Pass this Matrix to Setup to set the algebra rules.

Setup(%);

[algebrarules = {%Commutator(py, px) = 0, %Commutator(pz, px) = 0, %Commutator(pz, py) = 0, %Commutator(x, px) = I, %Commutator(x, py) = 0, %Commutator(x, pz) = 0, %Commutator(y, px) = 0, %Commutator(y, py) = I, %Commutator(y, pz) = 0, %Commutator(y, x) = 0, %Commutator(z, px) = 0, %Commutator(z, py) = 0, %Commutator(z, pz) = I, %Commutator(z, x) = 0, %Commutator(z, y) = 0}]

(6.1.5.12)

Set, for instance, the values of , and , the components of .

>	Lx := _i . L_;

(6.1.5.13)

>	Ly := _j . L_;

(6.1.5.14)

>	Lz := _k . L_;

(6.1.5.15)

Verify the commutator algebra for these components of .

>	Commutator(Lx, Ly) = I * Lz;

(6.1.5.16)

>	Commutator(Lz, Lx) = I * Ly;

(6.1.5.17)

>	Commutator(Ly, Lz) = I * Lx;

(6.1.5.18)

The three equations above are identically true.

Other operators frequently used in different contexts are the Annihilation and Creation operators: they augment or diminish the value of a quantum number by one. These operators are suitable, for instance, for working with multi-particle vector states; in that context the quantum numbers are called occupation numbers. This constructs a pair of annihilation/creation operators acting on the basis A involving only one quantum number.

>	am := Annihilation(A);

(6.1.5.19)

>	ap := Creation(A);

(6.1.5.20)

Annihilation and Creation operators act on Kets belonging to discrete bases and assume that the "lower" state happens when the quantum number is equal to zero (frequently called "vacuum": a ket with occupation number equal to zero represents a state with "no particles").

>	am . Ket(A, n);

(6.1.5.21)

am . %;

(6.1.5.22)

>	ap . Ket(A, n);

(6.1.5.23)

ap . %;

(6.1.5.24)

The Commutator of the operators and are automatically set when these operators are constructed, and satisfy (note the use of the inert form %Commutator):

>	%Commutator(am, ap);

(6.1.5.25)

value(%);

(6.1.5.26)

To indicate that the Kets of a basis are fermionic, use an anticommutative variable to label the basis. To set the prefix identifier of anticommutative variables use the Setup command.

>	Setup(anticommutativeprefix = Theta);

$[anticommutativeprefix = {Theta, _lambda}]$

(6.1.5.27)

>	type(Theta, anticommutative);

(6.1.5.28)

>	Ket(Theta);

(6.1.5.29)

Construct Annihilation and Creation operators acting on this basis; use the option notation = explicit so that the basis and the quantum numbers onto which these operators act are explicit.

>	Am := Annihilation(Theta, notation = explicit);

(6.1.5.30)

>	Ap := Creation(Theta, notation = explicit);

(6.1.5.31)

The AntiCommutator of these operators satisfy

>	%AntiCommutator(Am, Ap);

(6.1.5.32)

value(%);

(6.1.5.33)

According to Pauli's exclusion principle, only one fermionic particle can be in a given state, so starting from the vacuum,

>	Ket(Theta, 0);

(6.1.5.34)

Ap . %;

(6.1.5.35)

Ap . %;

(6.1.5.36)

Am . %%;

(6.1.5.37)

And as is always the case, the annihilation operator acting on the vacuum returns zero

Am . %;

(6.1.5.38)

implies

>	restart; with(Physics); Setup(math = true);

(6.2.1)

Consider two conjugate observables Q, P, and the corresponding Hermitian operators satisfying .

>	macro(h = `&hbar;`):

>	Setup(Hermitian = {Q,P}, %Commutator(Q,P) = I*h);

(6.2.2)

Suppose now that the system where and act is in some state normalized to 1, and set as the default state for computing Brackets.

>	Setup(%Bracket(psi, psi) = 1, bracketbasis = psi);

(6.2.3)

The mean values of the operators and in the state are then given by:

>	Qm := Bracket(Q); #Shortcut for Bracket(psi, Q, psi) after having set the bracketbasis to psi

(6.2.4)

>	Pm := Bracket(P);

(6.2.5)

Let's introduce another Hermitian operator, , and denote and the operators representing the observable deviations from these mean values by and .

>	Setup(her = Delta);

(6.2.6)

>	DQ := Delta(Q) = Q - Bracket(Q);

(6.2.7)

>	DP := Delta(P) = P - Bracket(P);

(6.2.8)

The value of the Commutator between and is a consequence of the value of the Commutator between and , and so it can be computed by rewriting the deviations in terms of and .

>	%Commutator(Delta(Q), Delta(P));

(6.2.9)

>	eval(%, {DQ, DP});

(6.2.10)

value(%);

(6.2.11)

Track this result as an algebra rule, so that in what follows we compute directly with and .

>	Setup(%%% = %);

(6.2.12)

To show now that implies , consider the action of these deviation operators and on the state of the system , and construct with them a new Ket involving a real parameter .

>	Ket(Psi, lambda) := (Delta(Q) + IlambdaDelta(P)) . Ket(psi);

(6.2.13)

The square of the norm of is

>	Dagger(%) . %;

(6.2.14)

Simplify this norm, taking into account the commutators and the fact that is real.

>	Setup(alg, quiet);

(6.2.15)

>	Simplify(%%) assuming lambda::real;

(6.2.16)

This is a polynomial in of second degree; its discriminant is negative or zero.

>	discrim(%, lambda) <= 0;

(6.2.17)

isolating , we obtain the lower bound for .

>	isolate(%, h^2)/4;

(6.2.18)

Note that this result is a consequence of , which in turn is a consequence of , so that and too satisfy , and in fact the product of any two conjugate Hermitian operators, as well as of the root-mean square deviations of them, satisfy this inequality.

Angular momentum: and

1. Consider the angular momentum operators , , , and in quantum mechanics. We want to verify that the Commutator of with any of the components of is zero (see, for instance, Chapter VI of Cohen-Tannoudji). For that purpose, a 3-D vector-quantum-operator representation of is constructed with the Vectors package (vectorpostfix identifier is '_'), and , and as well as and and their components are set as quantum operators.

>	restart; with(Physics); with(Physics[Vectors]); Setup(math = true);

(6.3.1)

To set the components and as quantum operators, it suffices to set and .

>	Setup(quantumoperators={L_, L, r_, x,y,z, p_, p});

$[quantumoperators = {L, L_, p, p_, r_, x, y, z}]$

(6.3.2)

So for and for itself in terms of the vector operators and , you have

>	LL := L_ . L_;

(6.3.3)

>	L_ := r_ &x p_;

(6.3.4)

where

>	r_ := x_i + y_j + z*_k;

(6.3.5)

>	p_ := p[x]_i + p[y]_j + p[z]*_k;

(6.3.6)

The Commutator rules for the components of are a consequence of the Commutator rules for the components of and . These rules can be set by using the Setup command, entering all of the commutators between any two of , as in > Setup(algebrarules = %Commutator(x, p[x]) = I, %Commutator(x, p[y]) = 0, %Commutator(x, y) = 0, ... ). A convenient alternative for situations such as this, where there are many commutators to be entered, is to work with indexed (tensor) notation (see problem 2 below) or create an indexing procedure for a Matrix. For example,

>	algebra := (i,j) -> ( %Commutator(Component(r_, i), Component(p_, j)) = I*KroneckerDelta[i,j], %Commutator(Component(r_, i), Component(r_, j)) = 0, %Commutator(Component(p_, i), Component(p_, j)) = 0);

(6.3.7)

The commutators are then generated by the Matrix constructor and the whole Matrix can be passed to Setup.

>	Matrix(3, 3, algebra);

(6.3.8)

Setup(%);

(6.3.9)

The components of are:

>	L[x] := Component(L_, 1);

(6.3.10)

>	L[y] := Component(L_, 2);

(6.3.11)

>	L[z] := Component(L_, 3);

(6.3.12)

To verify that commutes with each of , and , an expansion of the Commutator is not sufficient.

>	%Commutator(LL, L[x]); # inert representation of Commutator(LL, Lx)

(6.3.13)

>	value(%); # activates the inert computation

(6.3.14)

>	zero := expand(%);

(6.3.15)

To verify that the above is actually equal to 0, these commutator rules:

>	Setup(alg);

(6.3.16)

must be taken into account. For that purpose, use Simplify.

>	Simplify(zero);

(6.3.17)

>	%Commutator(LL, Ly);

(6.3.18)

>	Simplify( % );

(6.3.19)

>	Simplify( Commutator(LL, Lz) );

(6.3.20)

______________________________________________________________

2. Using tensor notation to represent the quantum operator components of , show that (see the exercises of Chap VI in Cohen-Tannoudji).

For this purpose, set the dimension of spacetime to 3 and its signature to Euclidean, so that "spacetime" tensors are actually 3-D space tensors. To follow textbook notation, use also lowercaselatin tensor indices (see Setup).

>	Setup(dimension = 3, signature = `+`, spacetimeindices = lowercaselatin, quiet);

(6.3.21)

To set Commutator rules for r and p using tensor notation and have Simplify tackling their products using Einstein's summation convention for repeated indices, Define r and p as tensors of this 3-D Euclidean space.

>	Define(r, p);

${p, r, Physics:-D_[mu], Physics:-Dgamma[mu], Physics:-Psigma[mu], Physics:-Ricci[mu, nu], Physics:-Riemann[mu, nu, alpha, beta], Physics:-Weyl[mu, nu, alpha, beta], Physics:-d_[mu], Physics:-g_[mu, nu], Physics:-Christoffel[mu, nu, alpha], Physics:-Einstein[mu, nu], Physics:-KroneckerDelta[mu, nu], Physics:-LeviCivita[alpha, mu, nu]}$

(6.3.22)

Now set the related Commutator rules for the algebra in tensor notation; in doing so, erase previous settings for quantum operators and algebra rules (by using the redo option of Setup; this erasure of previous definitions is not necessary in this example, but it is sometimes desired).

>	Setup(redo, quantumoperators = '{L, L_, r, p}', algebrarules = {%Commutator(r[i], p[j]) = I*kd_[i,j], %Commutator(r[i], r[j]) = 0, %Commutator(p[i], p[j]) = 0});

$[algebrarules = {%Commutator(p[i], p[j]) = 0, %Commutator(r[i], p[j]) = I*Physics:-KroneckerDelta[i, j], %Commutator(r[i], r[j]) = 0}, quantumoperators = {L, L_, p, r}]$

(6.3.23)

Note in the above how compact the algebra rules are when written in tensor notation: you have only three instead of fifteen Commutator rules. Verify how these algebra rules work: