Example Worksheet - Maple Help
For the best experience, we recommend viewing online help using Google Chrome or Mozilla Firefox.

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : Physics : Example Worksheet

Examples in Physics 

Vectors and Analytical Geometry 

Vectorial equation of a line 

 

Problem 

  • Derive the vectorial equation of a straight line which passes through two generic points, A and B.
 

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

  • Generate a 3D plot of this line.
 

 

Solution  

The vectorial equation of a line is the equation which is satisfied by the position vector of any point on this line. We want to express this equation in terms of the position vectors of the points A and B. The parametric equation can then be constructed from the vectorial equation by taking the coefficients of the unit vectors. 

 

> restart; -1; with(Physics); -1; with(Physics[Physics:-Vectors]); -1
 

1. To construct the vectorial equation of the line, let A_; and B_; be vectors pointing to A and B. The vector difference A_; - B_; is parallel to the line we want to represent. 

> Physics:-Vectors:-`+`(A_, `+`(`-`(B_)));
 

`+`(A_, `-`(B_)) (1)
 

Now let r_; be the position vector of any point of this line; the vector Physics:-Vectors:-`+`(r_, `+`(`-`(A_))); is also parallel to the line so their cross product is equal to 0 (you can use &x, or the operator from the palette of Common Symbols) 

> Eq := Typesetting:-delayCrossProduct(Physics:-Vectors:-`+`(r_, `+`(`-`(A_))), `+`(A_, `-`(B_))) = 0;
 

Typesetting:-mprintslash([Eq := Physics:-Vectors:-`&x`(`+`(r_, `-`(A_)), `+`(A_, `-`(B_))) = 0], [Physics:-Vectors:-`&x`(`+`(r_, `-`(A_)), `+`(A_, `-`(B_))) = 0]) (2)
 

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

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

> A_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(2, _i), `+`(`-`(Physics:-`*`(3, _j)))), Physics:-`*`(4, _k));
 

Typesetting:-mprintslash([A_ := `+`(`*`(2, `*`(_i)), `-`(`*`(3, `*`(_j))), `*`(4, `*`(_k)))], [`+`(`*`(2, `*`(_i)), `-`(`*`(3, `*`(_j))), `*`(4, `*`(_k)))]) (3)
 

> B_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(`+`(`-`(Physics:-`*`(4, _i))), Physics:-`*`(2, _j)), `+`(`-`(_k)));
 

Typesetting:-mprintslash([B_ := `+`(`-`(`*`(4, `*`(_i))), `*`(2, `*`(_j)), `-`(_k))], [`+`(`-`(`*`(4, `*`(_i))), `*`(2, `*`(_j)), `-`(_k))]) (4)
 

Regardless of the values of A_; and B_; , the position vector of a point is 

> r_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x, _i), Physics:-`*`(y, _j)), Physics:-`*`(z, _k));
 

Typesetting:-mprintslash([r_ := `+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)), `*`(_k, `*`(z)))], [`+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)), `*`(_k, `*`(z)))]) (5)
 

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

> Eq;
 

`+`(`*`(`+`(`*`(5, `*`(y)), `-`(5), `*`(5, `*`(z))), `*`(_i)), `*`(`+`(`*`(6, `*`(z)), `-`(14), `-`(`*`(5, `*`(x)))), `*`(_j)), `*`(`+`(`-`(`*`(5, `*`(x))), `-`(8), `-`(`*`(6, `*`(y)))), `*`(_k))) = 0 (6)
 

Since the unit vectors [_i, _j, _k; ] are independent, their Coefficients in the vectorial equation constitute the parametric equations; you can get those coefficients by taking the scalar product 

> Typesetting:-delayDotProduct([_i, _j, _k], Eq);
 

[`+`(`*`(5, `*`(y)), `-`(5), `*`(5, `*`(z))) = 0, `+`(`*`(6, `*`(z)), `-`(14), `-`(`*`(5, `*`(x)))) = 0, `+`(`-`(`*`(5, `*`(x))), `-`(8), `-`(`*`(6, `*`(y)))) = 0] (7)
 

We now have a system of three equations in terms of three unknowns x, y, and z, but it still represents the equation of a line, so two of the variables can be expressed in terms of the third one, which acts as the parameter of the parametric equations. Thus, by solving the system you get the following parametric equations 

> solve([`+`(`*`(5, `*`(y)), `-`(5), `*`(5, `*`(z))) = 0, `+`(`*`(6, `*`(z)), `-`(14), `-`(`*`(5, `*`(x)))) = 0, `+`(`-`(`*`(5, `*`(x))), `-`(8), `-`(`*`(6, `*`(y)))) = 0], [x, y, z]);
 

 

Warning, solve may be ignoring assumptions on the input variables.
[[x = `+`(`-`(`/`(14, 5)), `*`(`/`(6, 5), `*`(z))), y = `+`(1, `-`(z)), z = z]] (8)
 

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 select the curve's parameter: 

> select(evalb, [[x = `+`(`-`(`/`(14, 5)), `*`(`/`(6, 5), `*`(z))), y = `+`(1, `-`(z)), z = z]][1]);
 

[z = z] (9)
 

Now construct the appropriate input for the spacecurve command, say for z from -4 to 4 

> input := [op(map(rhs, [[x = `+`(`-`(`/`(14, 5)), `*`(`/`(6, 5), `*`(z))), y = `+`(1, `-`(z)), z = z]][1]))], lhs([z = z][1]) = -4 .. 4;
 

Typesetting:-mprintslash([input := [`+`(`-`(`/`(14, 5)), `*`(`/`(6, 5), `*`(z))), `+`(1, `-`(z)), z], z = -4 .. 4], [[`+`(`-`(`/`(14, 5)), `*`(`/`(6, 5), `*`(z))), `+`(1, `-`(z)), z], z = -4 .. 4]) (10)
 

Some options for improving the visualization are 

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

Typesetting:-mprintslash([opts := axes = boxed, scaling = constrained, orientation = [-130, 70], labels = [x, y, z]], [axes = boxed, scaling = constrained, orientation = [-130, 70], labels = [x, y, z]... (11)
 

> plots[spacecurve](input, opts);
 

Plot_2d
 

Vectorial equation of a plane 

 

Problem 

  • Derive the vectorial equation of a plane which passes through three generic points A, B, and C.
 

  • Choose three concrete points A, B, and C on this plane and plot it.
 

 

Solution 

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

> restart; -1; with(Physics[Physics:-Vectors]); -1
 

1. To construct this equation, let A_; , B_; , and C_; be the vectors pointing to A, B, and C, and let r_; be the position vector of any point on this plane. The differences Physics:-Vectors:-`+`(A_, `+`(`-`(C_))); and Physics:-Vectors:-`+`(A_, `+`(`-`(B_))); are vectors parallel to the plane we want to represent, and also the difference between r_;  and any of A_, B_; , or C_; is a vector parallel to the plane. So the equation of the plane can be obtained by first taking the cross product of any difference involving A_, B_; , and C_; to construct a vector perpendicular to the plane (you can use &x, or the operator from the palette of Common Symbols) 

> G_ := Typesetting:-delayCrossProduct(Physics:-Vectors:-`+`(A_, `+`(`-`(B_))), Physics:-Vectors:-`+`(A_, `+`(`-`(C_))));
 

Typesetting:-mprintslash([G_ := Physics:-Vectors:-`&x`(`+`(A_, `-`(B_)), `+`(A_, `-`(C_)))], [Physics:-Vectors:-`&x`(`+`(A_, `-`(B_)), `+`(A_, `-`(C_)))]) (12)
 

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

> Eq := Typesetting:-delayDotProduct(Physics:-Vectors:-`+`(r_, `+`(`-`(A_))), G_) = 0;
 

Typesetting:-mprintslash([Eq := Physics:-Vectors:-`.`(`+`(r_, `-`(A_)), Physics:-Vectors:-`&x`(`+`(A_, `-`(B_)), `+`(A_, `-`(C_)))) = 0], [Physics:-Vectors:-`.`(`+`(r_, `-`(A_)), Physics:-Vectors:-`&x... (13)
 

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

> A_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(2, _i), `+`(`-`(Physics:-`*`(3, _j)))), Physics:-`*`(4, _k));
 

Typesetting:-mprintslash([A_ := `+`(`*`(2, `*`(_i)), `-`(`*`(3, `*`(_j))), `*`(4, `*`(_k)))], [`+`(`*`(2, `*`(_i)), `-`(`*`(3, `*`(_j))), `*`(4, `*`(_k)))]) (14)
 

> B_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(5, _i), Physics:-`*`(4, _j)), `+`(`-`(Physics:-`*`(7, _k))));
 

Typesetting:-mprintslash([B_ := `+`(`*`(5, `*`(_i)), `*`(4, `*`(_j)), `-`(`*`(7, `*`(_k))))], [`+`(`*`(5, `*`(_i)), `*`(4, `*`(_j)), `-`(`*`(7, `*`(_k))))]) (15)
 

> C_ := Physics:-Vectors:-`+`(Physics:-`*`(30, Physics:-Vectors:-`^`(4, -1), A_), Physics:-`*`(90, Physics:-Vectors:-`^`(7, -1), B_));
 

Typesetting:-mprintslash([C_ := `+`(`*`(`/`(555, 7), `*`(_i)), `*`(`/`(405, 14), `*`(_j)), `-`(`*`(60, `*`(_k))))], [`+`(`*`(`/`(555, 7), `*`(_i)), `*`(`/`(405, 14), `*`(_j)), `-`(`*`(60, `*`(_k))))]) (16)
 

For r_; we always have 

> r_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x, _i), Physics:-`*`(y, _j)), Physics:-`*`(z, _k));
 

Typesetting:-mprintslash([r_ := `+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)), `*`(_k, `*`(z)))], [`+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)), `*`(_k, `*`(z)))]) (17)
 

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

> Eq;
 

`+`(`-`(`*`(`/`(1355, 14), `*`(x))), `-`(`*`(`/`(4607, 7), `*`(y))), `-`(`*`(`/`(6233, 14), `*`(z)))) = 0 (18)
 

To verify that the surface represented by this equation contains the points A, B, and C, you can substitute the values of the coordinates of these points and check if the equation is satisfied. For example, these are the coordinates of the first point: 

> Typesetting:-delayDotProduct([_i, _j, _k], A_);
 

[2, -3, 4] (19)
 

The following shows that, for these points A_, B_, C_; , the equation is satisfied: 

> for P in [A_, B_, C_] do eval(Eq, [x = Typesetting:-delayDotProduct(_i, P), y = Typesetting:-delayDotProduct(_j, P), z = Typesetting:-delayDotProduct(_k, P)]) end do;
 

 

 

0 = 0
0 = 0
0 = 0 (20)
 

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 with the command implicitplot3d. 

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

Typesetting:-mprintslash([opts := axes = boxed, scaling = constrained, orientation = [125, 65], style = surface], [axes = boxed, scaling = constrained, orientation = [125, 65], style = surface]) (21)
 

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

Plot_2d
 

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 on it. 

 

> restart; -1; with(Physics[Physics:-Vectors]); -1
 

Let r_; represent the position vector of any point on the plane, A_; be a vector pointing to the center of the sphere, and B_; be a vector pointing to the point B where the plane is tangent to the sphere. So the difference Physics:-Vectors:-`+`(r_, `+`(`-`(B_))); is a vector on the plane, and the difference Physics:-Vectors:-`+`(B_, `+`(`-`(A_))); 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 := Typesetting:-delayDotProduct(Physics:-Vectors:-`+`(r_, `+`(`-`(B_))), Physics:-Vectors:-`+`(B_, `+`(`-`(A_)))) = 0;
 

Typesetting:-mprintslash([Eq := Physics:-Vectors:-`.`(`+`(r_, `-`(B_)), `+`(B_, `-`(A_))) = 0], [Physics:-Vectors:-`.`(`+`(r_, `-`(B_)), `+`(B_, `-`(A_))) = 0]) (22)
 

This is already the vectorial equation of the tangent plane, only it is not yet expressed in terms of the radius of the sphere. Now since Physics:-Vectors:-`+`(B_, `+`(`-`(A_))); 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 := Physics:-Vectors:-Norm(Physics:-Vectors:-`+`(B_, `+`(`-`(A_)))) = a;
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (23)
 

Expand both expressions to use them together. 

> expand(Eq);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (24)
 

> map(proc (u) options operator, arrow; Physics:-Vectors:-`^`(u, 2) end proc, key);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msup(Typesetting:-mrow(Typesetting:-mi( (25)
 

> expand(`*`(`^`(Physics:-Vectors:-Norm(`+`(B_, `-`(A_))), 2)) = 1764);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msup(Typesetting:-mrow(Typesetting:-mi( (26)
 

So simplify one with respect to the other one, eliminating Physics:-Vectors:-Norm(B_); (see simplify/siderels) 

> simplify(`+`(Physics:-Vectors:-`.`(r_, B_), `-`(Physics:-Vectors:-`.`(r_, A_)), `-`(`*`(`^`(Physics:-Vectors:-Norm(B_), 2))), Physics:-Vectors:-`.`(B_, A_)) = 0, {`+`(`*`(`^`(Physics:-Vectors:-Norm(B_...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msup(Typesetting:-mrow(Typesetting:-mi( (27)
 

After collecting terms in the above, the requested vectorial equation can be more compactly rewritten as Typesetting:-delayDotProduct(Physics:-Vectors:-`+`(B_, `+`(`-`(A_))), Physics:-Vectors:-`+`(r_, `+`(`-`(A_)))) = Physics:-Vectors:-`^`(a, 2); .  

> (rhs = lhs)(expand(Typesetting:-delayDotProduct(Physics:-Vectors:-`+`(B_, `+`(`-`(A_))), Physics:-Vectors:-`+`(r_, `+`(`-`(A_)))) = Physics:-Vectors:-`^`(a, 2)));
 

Typesetting:-mrow(Typesetting:-mn( (28)
 

> subs(1764 = `+`(Physics:-Vectors:-`.`(r_, B_), `-`(Physics:-Vectors:-`.`(B_, A_)), `-`(Physics:-Vectors:-`.`(r_, A_)), `*`(`^`(Physics:-Vectors:-Norm(A_), 2))), `+`(`*`(`^`(Physics:-Vectors:-Norm(A_),...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msup(Typesetting:-mrow(Typesetting:-mi( (29)
 

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 derived equation and plot the sphere and 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 r_ = R_(u, v, w); 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 the equation as  Physics:-`*`(Physics:-Vectors:-`^`(d, 3), r_) = Physics:-`*`(Typesetting:-delayDotProduct(Physics:-Vectors:-diff(r_, u), Typesetting:-delayCrossProduct(Physics:-Vectors:-diff(r_, v), Physics:-Vectors:.... 

> restart; -1; with(Physics:-Vectors); -1
 

We want this volume element to be expressed in spherical coordinates (r, phi, theta; ); we can always choose these coordinates themselves as parameters u, v, w. Therefore, we are interested in the explicit form of (note the use of  %diff, the inert form diff, and for the cross product you can use &x, or the operator from the palette of Common Symbols) 

> answer := Typesetting:-delayDotProduct(%diff(r_, r), Typesetting:-delayCrossProduct(%diff(r_, theta), %diff(r_, phi)));
 

Typesetting:-mprintslash([answer := Physics:-Vectors:-`.`(%diff(r_, r), Physics:-Vectors:-`&x`(%diff(r_, theta), %diff(r_, phi)))], [Physics:-Vectors:-`.`(%diff(r_, r), Physics:-Vectors:-`&x`(%diff(r_... (30)
 

The first step is to write the vectorial equation r_ = R_(r, phi, theta); for a sphere of radius r; that is the equation satisfied by the position vector r_; of any point of the sphere. In spherical coordinates, with the origin of the reference system at the center of the sphere, this vectorial equation has its simplest form where r_; points to any point of the sphere, r is the radial coordinate (constant over the sphere) and _r;  is the radial unit vector. So 

> r_ := Physics:-`*`(r, _r);
 

Typesetting:-mprintslash([r_ := `*`(r, `*`(_r))], [`*`(r, `*`(_r))]) (31)
 

From which the value of 

> answer;
 

Typesetting:-mprintslash([Physics:-Vectors:-`.`(%diff(`*`(r, `*`(_r)), r), Physics:-Vectors:-`&x`(%diff(`*`(r, `*`(_r)), theta), %diff(`*`(r, `*`(_r)), phi)))], [Physics:-Vectors:-`.`(%diff(`*`(r, `*`... (32)
 

can be computed directly, using the value command 

> value(answer);
 

`*`(`^`(r, 2), `*`(sin(theta))) (33)
 

Alternatively, one could compute this result one step at a time by making it explicit that _r; depends on varphi; and For this purpose, change the basis in the vectorial equation to the Cartesian basis  (_i, _j, _k; ), where all the unit vectors are constant and so the partial derivatives can be performed directly. 

> r_ := Physics:-Vectors:-ChangeBasis(r_, 1);
 

Typesetting:-mprintslash([r_ := `+`(`*`(r, `*`(sin(theta), `*`(cos(phi), `*`(_i)))), `*`(r, `*`(sin(phi), `*`(sin(theta), `*`(_j)))), `*`(r, `*`(cos(theta), `*`(_k))))], [`+`(`*`(r, `*`(sin(theta), `*... (34)
 

So, the answer introduced in (30) becomes: 

> answer;
 

Typesetting:-mprintslash([Physics:-Vectors:-`.`(%diff(`+`(`*`(r, `*`(sin(theta), `*`(cos(phi), `*`(_i)))), `*`(r, `*`(sin(phi), `*`(sin(theta), `*`(_j)))), `*`(r, `*`(cos(theta), `*`(_k)))), r), Physi...
Typesetting:-mprintslash([Physics:-Vectors:-`.`(%diff(`+`(`*`(r, `*`(sin(theta), `*`(cos(phi), `*`(_i)))), `*`(r, `*`(sin(phi), `*`(sin(theta), `*`(_j)))), `*`(r, `*`(cos(theta), `*`(_k)))), r), Physi...
(35)
 

> value(answer);
 

`*`(`^`(r, 2), `*`(sin(theta))) (36)
 

Hence the volume element is  

>
 

Parametrization of Curves, Surfaces and Volumes 

 

Consider the following C; representing a curve in space 

> restart; -1; with(Physics:-Vectors)
 

[`&x`, `+`, `.`, Assume, ChangeBasis, ChangeCoordinates, CompactDisplay, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, ParametrizeCurve, ParametrizeSurface,...
[`&x`, `+`, `.`, Assume, ChangeBasis, ChangeCoordinates, CompactDisplay, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, ParametrizeCurve, ParametrizeSurface,...
(37)
 

> C := {y = Physics:-Vectors:-`^`(x, 2), z = 0};
 

Typesetting:-mprintslash([C := {y = `*`(`^`(x, 2)), z = 0}], [{y = `*`(`^`(x, 2)), z = 0}]) (38)
 

The parametric equations for this curve are 

> Physics:-Vectors:-ParametrizeCurve(C);
 

[x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t (39)
 

The right-hand sides of equations above are the components of the position vector `#mover(mi( in cartesian coordinates, from where a vectorial form of these equations is 

> Physics:-Vectors:-ParametrizeCurve(C, output = vector);
 

`+`(`*`(_j, `*`(`^`(t, 2))), `*`(_i, `*`(t))), t (40)
 

The curve C; can also be passed in vector form 

> C := Physics:-Vectors:-`+`(Physics:-`*`(x, _i), Physics:-`*`(Physics:-Vectors:-`^`(x, 2), _j));
 

Typesetting:-mprintslash([C := `+`(`*`(_j, `*`(`^`(x, 2))), `*`(_i, `*`(x)))], [`+`(`*`(_j, `*`(`^`(x, 2))), `*`(_i, `*`(x)))]) (41)
 

> Physics:-Vectors:-ParametrizeCurve(C);
 

[x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t (42)
 

The equations of circle of radius a; on the x, y; plane can be written as 

> C := {z = 0, Physics:-Vectors:-`+`(Physics:-Vectors:-`^`(x, 2), Physics:-Vectors:-`^`(y, 2)) = Physics:-Vectors:-`^`(a, 2)};
 

Typesetting:-mprintslash([C := {z = 0, `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = 1764}], [{z = 0, `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = 1764}]) (43)
 

> Physics:-Vectors:-ParametrizeCurve(C);
 

[x(phi) = `+`(`*`(42, `*`(cos(phi)))), y(phi) = `+`(`*`(42, `*`(sin(phi)))), z(phi) = 0], phi (44)
 

In vector notation, 

> Physics:-Vectors:-ParametrizeCurve(C, output = vector);
 

`+`(`*`(42, `*`(cos(phi), `*`(_i))), `*`(42, `*`(sin(phi), `*`(_j)))), phi (45)
 

Changing the basis of this vector, for example to cylindrical, we get 

> Physics:-Vectors:-ChangeBasis((`+`(`*`(42, `*`(cos(phi), `*`(_i))), `*`(42, `*`(sin(phi), `*`(_j)))), phi)[1], cylindrical);
 

`+`(`*`(42, `*`(_rho))) (46)
 

The same result can be obtained by specifying the basis 

> Physics:-Vectors:-ParametrizeCurve(C, output = vector, basis = cylindrical);
 

`+`(`*`(42, `*`(_rho))), phi (47)
 

To parametrize surfaces and volumes you can use ParametrizeSurface and ParametrizeVolume, which work in the same way as ParametrizeCurve, only with two and three parameters respectively. 

> C__2 := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`^`(x, 2), Physics:-Vectors:-`^`(y, 2)), Physics:-Vectors:-`^`(z, 2)) = Physics:-Vectors:-`^`(a, 2);
 

Typesetting:-mprintslash([C__2 := `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = 1764], [`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = 1764]) (48)
 

> Physics:-Vectors:-ParametrizeSurface(C__2);
 

[x(theta, phi) = `+`(`*`(42, `*`(sin(theta), `*`(cos(phi))))), y(theta, phi) = `+`(`*`(42, `*`(sin(theta), `*`(sin(phi))))), z(theta, phi) = `+`(`*`(42, `*`(cos(theta))))], [theta, phi] (49)
 

> C__3 := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`^`(x, 2), Physics:-Vectors:-`^`(y, 2)), Physics:-Vectors:-`^`(z, 2)) = Physics:-Vectors:-`^`(r, 2);
 

Typesetting:-mprintslash([C__3 := `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = `*`(`^`(r, 2))], [`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = `*`(`^`(r, 2))]) (50)
 

> Physics:-Vectors:-ParametrizeVolume(C__3);
 

[x(r, theta, phi) = `*`(r, `*`(sin(theta), `*`(cos(phi)))), y(r, theta, phi) = `*`(r, `*`(sin(theta), `*`(sin(phi)))), z(r, theta, phi) = `*`(r, `*`(cos(theta)))], [r, theta, phi] (51)
 

> Physics:-Vectors:-ParametrizeVolume(C__3, basis = spherical, output = vector);
 

`*`(r, `*`(_r)), [r, theta, phi] (52)
 

> Physics:-Vectors:-ParametrizeVolume(C__3, basis = cylindrical, output = vector);
 

`+`(`*`(_k, `*`(z)), `*`(_rho, `*`(rho))), [rho, phi, z] (53)
 

 

>
 

Integral Vector Calculus 

 

Integral vector calculus is about computing path, surface and volume vector integrals. Those are integrals where the integrand is a scalar or vector function, and the computation is done from any description (algebraic, parametric, vectorial) of the region of integration - a path, surface of volume. 

 

There are three kinds of line or path integrals: 

Physics:-Vectors:-int(F, r_ = A_ .. B_)[path = C];  

 

 

where `#mover(mi( and `#mover(mi( are points in space, the limits of integration, that belong to the curve C; over which the integral is performed. In the first integral, F; is a scalar function and the result of the integration is thus a vector. In the second and third integrals the integrand `#mover(mi( is a vector function, so that the dot product `#mrow(mover(mi( is a scalar and so is the integral, while in the third one `#mrow(mover(mi( is a vector and so is its integration over the region C; . 

 

Likewise, there are three kinds of surface integrals 

Physics:-Vectors:-int(F, S_)[surface = C, parameters = [u = a .. b, v = c .. d]];  

 

 

  • and two types of volume integrals
 

Physics:-Vectors:-int(F, V)[volume = C, parameters = [u = a .. b, v = c .. d, w = f .. g]];  

Physics:-Vectors:-int(F_, V)[volume = C, parameters = [u = a .. b, v = c .. d, w = f .. g]];  

 

The line element `#mrow(mi( in path integrals is expressed in terms of a parameter t; as 

`#mrow(mi( 

Using indexed notation for derivatives `#mrow(mscripts(mi(, the surface element `#mrow(mi( in surface integrals is expressed in terms of parameters as 

Physics:-`*`(`#mrow(mi( 

and the volume element `#mrow(mi( in volume integrals as 

Physics:-`*`(`#mrow(mi( 

 

The integrals in the three cases are computed by first expressing the integrand and the integration element in terms of the parameters using the parametric equations derived with ParametrizeCurve, ParametrizeSurface and ParametrizeVolume, then performing the vector product operations, then the integration.  

 

Vectorial Path integrals 

> restart; -1; with(Physics:-Vectors); 1
 

[`&x`, `+`, `.`, Assume, ChangeBasis, ChangeCoordinates, CompactDisplay, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, ParametrizeCurve, ParametrizeSurface,...
[`&x`, `+`, `.`, Assume, ChangeBasis, ChangeCoordinates, CompactDisplay, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, ParametrizeCurve, ParametrizeSurface,...
(54)
 

Consider the following scalar function F; , path C; and integration limits `#mover(mi( and `#mover(mi( 

> F := Physics:-`*`(x, y);
 

Typesetting:-mprintslash([F := `*`(x, `*`(y))], [`*`(x, `*`(y))]) (55)
 

> C := {y = Physics:-Vectors:-`^`(x, 2), z = 0};
 

Typesetting:-mprintslash([C := {y = `*`(`^`(x, 2)), z = 0}], [{y = `*`(`^`(x, 2)), z = 0}]) (56)
 

> A_ := [0, 0, 0];
 

Typesetting:-mprintslash([A_ := [0, 0, 0]], [[0, 0, 0]]) (57)
 

> B_ := [1, 1, 0];
 

Typesetting:-mprintslash([B_ := [1, 1, 0]], [[1, 1, 0]]) (58)
 

The line or path integral, shown here in inert form on the left-hand side and active, computed to the end on the right-hand side, is 

> (Int = Physics:-Vectors:-int)(F, r_ = A_ .. B_, path = C);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (59)
 

The output on the right-hand side is a vector. Within int, to perform this integration the curve C = {y = Physics:-Vectors:-`^`(x, 2), z = 0}; is first parametrized using ParametrizeCurve 

> Physics:-Vectors:-ParametrizeCurve(C);
 

[x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t (60)
 

The output above is a sequence, first the parametric equations as an ordered list (order [x, y, z]; ) then the parameter, in this case t; . For the formulation of the integral to make sense, the limits of integration `#mover(mi( and `#mover(mi( must belong to this curve, i.e. satisfy the parametric equations for some value of the parameter, in this case t = 0; and t = 1;  

> A_ = Eval(([x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t)[1], t = 0);
 

[0, 0, 0] = Eval([x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t = 0) (61)
 

> B_ = Eval(([x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t)[1], t = 1);
 

[1, 1, 0] = Eval([x(t) = t, y(t) = `*`(`^`(t, 2)), z(t) = 0], t = 1) (62)
 

To see the integral after being parametrized and before performing the integration you can use the option inert 

> (Int = Physics:-Vectors:-int)(F, r_ = A_ .. B_, path = C, inert);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (63)
 

Since there is a relation one-to-one between the integration limits `#mover(mi( and `#mover(mi( and the parameter's range, instead of indicating `#mover(mi( and `#mover(mi( you can also indicate the range of t; itself, getting the same result (7) 

> (Int = Physics:-Vectors:-int)(F, r_, path = C, parameter = {t = 0 .. 1});
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (64)
 

When the parameter's range is passed, you can also use a shortcut notation, passing the second argument as an equation r_ = C, making more explicit that r_ is the vector representation of the region of integration C; , so the line element `#mrow(mi( is the differential of the parametric equations of C. In this case indicating path = C is redundant and can be omitted. 

> (Int = Physics:-Vectors:-int)(F, r_ = C, parameter = {t = 0 .. 1});
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (65)
 

The integration path and the limits of integration can be expressed in vector notation as well 

> C_ := Physics:-Vectors:-`+`(Physics:-`*`(x, _i), Physics:-`*`(Physics:-Vectors:-`^`(x, 2), _j));
 

Typesetting:-mprintslash([C_ := `+`(`*`(_j, `*`(`^`(x, 2))), `*`(_i, `*`(x)))], [`+`(`*`(_j, `*`(`^`(x, 2))), `*`(_i, `*`(x)))]) (66)
 

> A_ := 0;
 

Typesetting:-mprintslash([A_ := 0], [0]) (67)
 

> B_ := Physics:-Vectors:-`+`(_i, _j);
 

Typesetting:-mprintslash([B_ := `+`(_i, _j)], [`+`(_i, _j)]) (68)
 

> (Int = Physics:-Vectors:-int)(F, r_ = A_ .. B_, path = C_);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (69)
 

Vectorial Surface and Volume integrals 

The case of surface and volume integrals is analogous to that of line (path) integrals, but for two things: instead of one, there are two or three parameters, and instead of indicating integration limits, it is required that you indicate the parameter's ranges. 

 

The following C; represents the surface of a sphere of radius a; centered at the origin 

> Physics:-Vectors:-Assume(`<`(0, a));
 

Error, (in Physics:-Assume) {} is a set without a property; it is not a valid assumption
 

> C__2 := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`^`(x, 2), Physics:-Vectors:-`^`(y, 2)), Physics:-Vectors:-`^`(z, 2)) = Physics:-Vectors:-`^`(a, 2);
 

Typesetting:-mprintslash([C__2 := `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = 1764], [`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = 1764]) (70)
 

> Physics:-Vectors:-ParametrizeSurface(C__2);
 

[x(theta, phi) = `+`(`*`(42, `*`(sin(theta), `*`(cos(phi))))), y(theta, phi) = `+`(`*`(42, `*`(sin(theta), `*`(sin(phi))))), z(theta, phi) = `+`(`*`(42, `*`(cos(theta))))], [theta, phi] (71)
 

From the definition of the vectorial surface element as 

Physics:-`*`(`#mrow(mi( 

`#mrow(mi( is a vector perpendicular to the surface of the sphere with magnitude equal to the radius a; . Hence, the vectorial surface integral of `#mrow(mi( over the upper half of the sphere should be a vector perpendicular to the plane x, y; in the `#mover(mi( direction, with modulus equal to half the surface of a sphere specified by C[2]; . Also, checking the form of the parametrization returned by ParametrizeSurface is useful to understand what the ranges for the parameters theta; and phi; need to be in order to represent the desired region. 

> (Int = Physics:-Vectors:-int)(1, S_, surface = C__2, parameters = [theta = 0 .. Physics:-`*`(Pi, Physics:-Vectors:-`^`(2, -1)), phi = 0 .. Physics:-`*`(2, Pi)]);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (72)
 

From symmetry considerations, the integral of `#mrow(mi( over the lower half of the sphere should have the same magnitude but opposite direction (`+`(`-`(`#mover(mi(), from where the integral over the whole sphere, so for theta = 0 .. Pi; , should be equal to 0 

> (Int = Physics:-Vectors:-int)(1, S_, surface = C__2, parameters = [theta = 0 .. Pi, phi = 0 .. Physics:-`*`(2, Pi)]);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (73)
 

From this example we see that to get the area computing the integral of `#mrow(mi( over the whole surface it is necessary to take as integrand the modulus of `#mrow(mi(, that is its scalar product with a unit vector parallel to it. By definition of surface element, Physics:-`*`(`#mrow(mi(, so that the unit vector _n_ = Physics:-`*`(d, S_, Physics:-Vectors:-`^`(Physics:-Vectors:-Norm(Physics:-`*`(d, S_)), -1)); is given by 

> Physics:-Vectors:-CompactDisplay(r_(theta, phi));
 

Typesetting:-mprintslash([`*`(r_(theta, phi), `*`(`will now be displayed as`, `*`(r_)))], [`*`(r_(theta, phi), `*`(`will now be displayed as`, `*`(r_)))]) (74)
 

> _n_ := Physics:-`*`(Physics:-Vectors:-`&x`(Physics:-Vectors:-diff(r_(theta, phi), theta), Physics:-Vectors:-diff(r_(theta, phi), phi)), Physics:-Vectors:-`^`(Physics:-Vectors:-Norm(Physics:-Vectors:-`...
 

Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (75)
 

where `#mover(mi( is the vectorial form of the parametric equations of the surface. It is easy to see this vector is the radial unit vector _r; . For that purpose you can use the option output = vector; of ParametrizeSurface to get the vectorial form of the parametric equations for C[2];  

> Physics:-Vectors:-ParametrizeSurface(C__2, output = vector);
 

`+`(`*`(42, `*`(sin(theta), `*`(cos(phi), `*`(_i)))), `*`(42, `*`(sin(theta), `*`(sin(phi), `*`(_j)))), `*`(42, `*`(cos(theta), `*`(_k)))), [theta, phi] (76)
 

Introduce this value of `#mover(mi( into the expression for _n_; and change basis to spherical 

> _n_ := eval(_n_, r_(theta, phi) = (`+`(`*`(42, `*`(sin(theta), `*`(cos(phi), `*`(_i)))), `*`(42, `*`(sin(theta), `*`(sin(phi), `*`(_j)))), `*`(42, `*`(cos(theta), `*`(_k)))), [theta, phi])[1]);
 

Typesetting:-mprintslash([_n_ := `+`(`/`(`*`(`/`(1, 882), `*`(`+`(`*`(1764, `*`(`^`(sin(theta), 2), `*`(cos(phi), `*`(_i)))), `*`(1764, `*`(`^`(sin(theta), 2), `*`(sin(phi), `*`(_j)))), `*`(882, `*`(s... (77)
 

> _n_ := Physics:-Vectors:-ChangeBasis(_n_, spherical);
 

Typesetting:-mprintslash([_n_ := _r], [_r]) (78)
 

Putting all together, the area of a sphere of radius a; is given by this closed surface integral 

> (Int = Physics:-Vectors:-int)(_n_, S_, surface = C__2, parameters = [theta = 0 .. Pi, phi = 0 .. Physics:-`*`(2, Pi)]);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (79)
 

In the case of the volume of a sphere, an algebraic representation of the region is 

> C__3 := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`^`(x, 2), Physics:-Vectors:-`^`(y, 2)), Physics:-Vectors:-`^`(z, 2)) = Physics:-Vectors:-`^`(r, 2);
 

Typesetting:-mprintslash([C__3 := `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = `*`(`^`(r, 2))], [`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2))) = `*`(`^`(r, 2))]) (80)
 

from where the volume of a sphere of radius a; is equal to the integral of the corresponding Physics:-`*`(`#mrow(mi( with parameters r = 0 .. a; , theta = 0 .. Pi; and phi = 0 .. Physics:-`*`(2, Pi);  

> (Int = Physics:-Vectors:-int)(1, V, volume = C__3, parameters = [r = 0 .. a, theta = 0 .. Pi, phi = 0 .. Physics:-`*`(2, Pi)]);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (81)
 

Using ParametrizeVolume one can also see that the vectorial representation of this region is just the position vector `#mover(mi( written in spherical coordinates and basis 

> Physics:-Vectors:-ParametrizeVolume(C__3, basis = spherical, output = vector);
 

`*`(r, `*`(_r)), [r, theta, phi] (82)
 

>
 

Vectors in Spherical Coordinates using Tensor Notation (advanced) 

 

The following is a topic that appears frequently in formulations: given a 3D vector in spherical (or any curvilinear) coordinates, how do you represent and relate, in simple terms, the vector and the corresponding vectorial operations Gradient, Divergence, Curl and Laplacian using tensor notation?  

 

I. The line element in spherical coordinates and the scale-factors 

 

Start by setting the spacetime to be 3-dimensional, Euclidean, and use Cartesian coordinates 

> restart; -1; with(Physics); -1; with(Physics[Physics:-Vectors]); -1
 

 

> Physics:-Vectors:-Setup(dimension = 3, coordinates = cartesian, metric = Euclidean, spacetimeindices = lowercaselatin); -1
 

 

 

 

 

 

 

Physics:-`*`(`Systems of spacetime coordinates are:`, {X = (x, y, z)});
_______________________________________________________;
Physics:-`*`(`The Euclidean metric in coordinates `, [x, y, z]);
_______________________________________________________;
Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(
_______________________________________________________; (83)
 

 

In vector calculus,the line element dr_; is at the root of everything, and in Cartesian coordinates it has the simple form 

> dr_ = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(dx, _i), Physics:-`*`(dy, _j)), Physics:-`*`(dz, _k));
 

dr_ = `+`(`*`(_i, `*`(dx)), `*`(_j, `*`(dy)), `*`(_k, `*`(dz))) (84)
 

To compute the line element  dr_; in spherical coordinates, the starting point is the transformation  

> tr := `~`[`=`]([X], ` $`, Physics:-Vectors:-ChangeCoordinates([X], spherical));
 

Typesetting:-mprintslash([tr := [x = `*`(r, `*`(sin(theta), `*`(cos(phi)))), y = `*`(r, `*`(sin(theta), `*`(sin(phi)))), z = `*`(r, `*`(cos(theta)))]], [[x = `*`(r, `*`(sin(theta), `*`(cos(phi)))), y ... (85)
 

> Physics:-Coordinates(S = [r, theta, phi]); -1
 

Physics:-`*`(`Systems of spacetime coordinates are:`, {S = (r, theta, phi), X = (x, y, z)}); (86)
 

Since in Physics:-`*`(dr_ = `+`(`*`(_i, `*`(dx)), `*`(_j, `*`(dy)), `*`(_k, `*`(dz))), [dx, dy, dz]); are just symbols with no relationship to [x, y, z]; , start transforming these differentials using the chain rule, computing the Jacobian of the transformation tr := `~`[`=`]([X], ` $`, Physics:-Vectors:-ChangeCoordinates([X], spherical)); . In this Jacobian J, the first line is [Physics:-`*`(`∂x`, Physics:-Vectors:-`^`(`∂r`, -1), dr), Physics:-`*`(Physics:-`*`(`∂x`, Physics:-Vectors:-`^`(Physics:-`*`(`∂`, theta), -1), d), theta), Physics:-...  

> J := VectorCalculus:-Jacobian(map(rhs, [x = `*`(r, `*`(sin(theta), `*`(cos(phi)))), y = `*`(r, `*`(sin(theta), `*`(sin(phi)))), z = `*`(r, `*`(cos(theta)))]), [S]); -1
 

So in matrix notation, 

> Vector([dx, dy, dz]) = Typesetting:-delayDotProduct(J, Vector([dr, dtheta, dphi]));
 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mi( (87)
 

To complete the computation of  dr_; in spherical coordinates we can now use ChangeBasis, provided that we next substitute Vector([dx, dy, dz]) = Typesetting:-delayDotProduct(J, Vector([dr, dtheta, dphi])); in the result, expressing the abstract objects [dx, dy, dz]; in terms of [dr, dtheta, dphi]; . 

 

In two steps: 

> lhs(dr_ = `+`(`*`(_i, `*`(dx)), `*`(_j, `*`(dy)), `*`(_k, `*`(dz)))) = Physics:-Vectors:-ChangeBasis(rhs(dr_ = `+`(`*`(_i, `*`(dx)), `*`(_j, `*`(dy)), `*`(_k, `*`(dz)))), spherical);
 

dr_ = `+`(`*`(`+`(`*`(dx, `*`(sin(theta), `*`(cos(phi)))), `*`(dy, `*`(sin(theta), `*`(sin(phi)))), `*`(dz, `*`(cos(theta)))), `*`(_r)), `*`(`+`(`*`(dx, `*`(cos(phi), `*`(cos(theta)))), `*`(dy, `*`(si...
dr_ = `+`(`*`(`+`(`*`(dx, `*`(sin(theta), `*`(cos(phi)))), `*`(dy, `*`(sin(theta), `*`(sin(phi)))), `*`(dz, `*`(cos(theta)))), `*`(_r)), `*`(`+`(`*`(dx, `*`(cos(phi), `*`(cos(theta)))), `*`(dy, `*`(si...
(88)
 

The line element 

>
 

dr_ = `+`(`*`(_phi, `*`(dphi, `*`(r, `*`(sin(theta))))), `*`(_theta, `*`(dtheta, `*`(r))), `*`(_r, `*`(dr))) (89)
 

This result is important: it gives us the so-called scale factors, the key that connects 3D vectors with the related covariant and contravariant tensors in curvilinear coordinates. The scale factors are computed from by taking the scalar product with each of the unit vectors [_r, _theta, _phi]; , then taking the coefficients of the differentials [dr, dtheta, dphi]; (just substitute them by the number 1) 

> h := subs(`~`[`=`]([dr, dtheta, dphi], 1), [seq(Typesetting:-delayDotProduct(rhs(dr_ = `+`(`*`(_phi, `*`(dphi, `*`(r, `*`(sin(theta))))), `*`(_theta, `*`(dtheta, `*`(r))), `*`(_r, `*`(dr)))), q), q = ...
 

Typesetting:-mprintslash([h := [1, r, `*`(r, `*`(sin(theta)))]], [[1, r, `*`(r, `*`(sin(theta)))]]) (90)
 

The scale factors are relevant because the components of the 3D vector and the corresponding tensor are not the same in curvilinear coordinates. For instance, representing the differential of the coordinates as the tensor Physics:-Vectors:-`^`(dS, j) = [dr, dtheta, dphi]; , we see that the corresponding vector, the line element in spherical coordinates dS_; , is not constructed by directly equating its components to the components of Physics:-Vectors:-`^`(dS, j) = [dr, dtheta, dphi]; , so   

`<>`(dS_, Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(dr, _r), Physics:-`*`(dtheta, _theta)), Physics:-`*`(dphi, _phi)));  

The vector dS_; is constructed multiplying these contravariant components [dr, dtheta, dphi]; by the scaling factors 

dS_ = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(h__r, dr, _r), Physics:-`*`(`h__θ`, dtheta, _theta)), Physics:-`*`(`h__φ`, dphi, _phi));  

This rule applies in general. The vectorial components of a 3D vector in an orthogonal system (curvilinear or not) are always expressed in terms of the contravariant components Physics:-Vectors:-`^`(A, j); the same way we did in the line above with the line element, using the scale-factors h__j; , so that 

A_ = Sum(Physics:-`*`(h[j], Physics:-Vectors:-`^`(A, j), _e__j_), j = 1 .. 3);  

where on the right-hand side we see the contravariant components Physics:-Vectors:-`^`(A, j); and the scale-factors h[j]; . Because the system is orthogonal, each vector component `#msub(mi( satisfies 

A__j = Physics:-`*`(h[j], A[`~j`]);  

The scale-factors h[j]; do not constitute a tensor, so on the right-hand side we do not sum over j.  Also, from  

Physics:-Vectors:-Norm(A_) = Physics:-`*`(A[j], A[`~j`]);  

it follows that, 

A__j = Physics:-`*`(A__j, Physics:-Vectors:-`^`(h__j, -1));  

where on the right-hand side we now have the covariant tensor components A__j; . 

 

  • This relationship between the components of a 3D vector and the contravariant and covariant components of a tensor representing the vector is key for translating vector-component to corresponding tensor-component formulas.
 

 

II. Transformation of contravariant and covariant tensors (can skip) 

 

Define two representations for the same tensor: A__c; will represent A in Cartesian coordinates, while A__s; will represent A in spherical coordinates. 

> Physics:-Define(A__c[j], A__s[j], quiet);
 

Typesetting:-mprintslash([{A__c[j], A__s[j], Physics:-Dgamma[b], Physics:-Psigma[b], Physics:-d_[b], Physics:-g_[b, c], Physics:-LeviCivita[b, c, d], S[b], X[b]}], [{A__c[j], A__s[j], Physics:-Dgamma[... (91)
 

 

Transformation rule for a contravariant tensor 

 

We know, by definition, that the transformation rule for the components of a contravariant tensor is `#mrow(msup(mi(, the same as the rule for the differential of the coordinates. Thus, the transformation rule from to computed using TransformCoordinates, should return the same relation as Vector([dx, dy, dz]) = Typesetting:-delayDotProduct(J, Vector([dr, dtheta, dphi])); . However, the application of the command requires one to pay attention because, just like in Vector([dx, dy, dz]) = Typesetting:-delayDotProduct(J, Vector([dr, dtheta, dphi])); , we want to isolate the Cartesian (not the spherical) components. This is like performing a reversed transformation. So we will use 

 

 

where on the left-hand side we isolate the three components of A in Cartesian coordinates, and on the right-hand side we transform the spherical components from spherical S = (r, theta, phi); (4th argument) to Cartesian X = (x, y, z); (3rd argument), which according to the 5th bullet of TransformCoordinates will result in a transformation expressed in terms of the old coordinates (here the spherical [S]; ). Expand things to compare with Vector([dx, dy, dz]) = Typesetting:-delayDotProduct(J, Vector([dr, dtheta, dphi]));  

 

> Vector[column](Physics:-TensorArray(A__c[`~j`])) = Physics:-TransformCoordinates(tr, A__s[`~j`], [X], [S], simplifier = expand);
 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-msub(Typesetting:-mi( (92)
 

We see that the transformation rule for a contravariant vector is, indeed, like the transformation Vector([dx, dy, dz]) = Typesetting:-delayDotProduct(J, Vector([dr, dtheta, dphi])); for the differential of the coordinates. 

 

Transformation rule for a covariant tensor 

 

For the transformation rule for the components of a covariant tensor A__c[j]; , we know by definition that it is `#mrow(msub(mi(, and so the same transformation rule for the gradient [Physics:-d_[x], Physics:-d_[y], Physics:-d_[z]]; , where Physics:-d_[x] = (proc (u) options operator, arrow; Physics:-Vectors:-diff(u, x) end proc); and so on. We can directly experiment with this by changing variables in the differential operators [Physics:-d_[x], Physics:-d_[y], Physics:-d_[z]];  

> Physics:-d_[x] = PDEtools:-dchange(tr, proc (u) options operator, arrow; Physics:-Vectors:-diff(u, x) end proc, simplify);
 

Physics:-d_[x] = (proc (u) options operator, arrow; `/`(`*`(`+`(`*`(diff(u, r), `*`(sin(theta), `*`(cos(phi), `*`(r)))), `*`(diff(u, theta), `*`(cos(phi), `*`(cos(theta)))), `-`(`*`(diff(u, phi), `*`(... (93)
 

This result, and the equivalent ones replacing x by y or z in the input above, can be computed in one go, in matrix and simplified form, using the Jacobian of the transformation computed after Physics:-Coordinates(S = [r, theta, phi]); -1. Now we need to take the transpose of the inverse of J (because we are now transforming the components of the gradient [Physics:-d_[x], Physics:-d_[y], Physics:-d_[z]]; ) 

> H := simplify(Physics:-Vectors:-`^`(Physics:-Vectors:-`^`(J, -1), %T)); -1
 

> Vector([Physics:-d_[x], Physics:-d_[y], Physics:-d_[z]]) = Typesetting:-delayDotProduct(H, Vector([Physics:-d_[r], Physics:-d_[theta], Physics:-d_[phi]]));
 

rtable(1 .. 3, [Physics:-d_[x], Physics:-d_[y], Physics:-d_[z]], subtype = Vector[column]) = rtable(1 .. 3, [`+`(`*`(sin(theta), `*`(cos(phi), `*`(Physics:-d_[r]))), `/`(`*`(cos(phi), `*`(cos(theta), ... (94)
 

The corresponding transformation equations relating the tensors A__c; and A__s; in Cartesian and spherical coordinates is computed with TransformCoordinates as in Vector[column](Physics:-TensorArray(A__c[`~j`])) = Physics:-TransformCoordinates(tr, A__s[`~j`], [X], [S], simplifier = expand); , just lower the indices on the left and right hand sides, i.e., remove the tilde ~ 

> Vector[column](Physics:-TensorArray(A__c[j])) = Physics:-TransformCoordinates(tr, A__s[j], [X], [S], simplifier = expand);
 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-msub(Typesetting:-mi( (95)
 

We see that the transformation rule for a covariant vector A__c[j]; is, indeed, just like the transformation rule H := simplify(Physics:-Vectors:-`^`(Physics:-Vectors:-`^`(J, -1), %T)); -1 for the gradient.  

 

Side Note: once the computation of these transformation rules is understood, we can also perform the inverse of Vector[column](Physics:-TensorArray(A__c[j])) = Physics:-TransformCoordinates(tr, A__s[j], [X], [S], simplifier = expand); as follows 

> Vector[column](Physics:-TensorArray(A__s[j])) = Physics:-TransformCoordinates(tr, A__c[j], [S], [X], simplifier = expand);
 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-msub(Typesetting:-mi( (96)
 

III. Deriving the transformation rule for the Gradient using TransformCoordinates (can skip) 

 

Turn the CompactDisplay notation for derivatives ON, so that the differentiation variable is displayed as an index: 

> ON;
 


The gradient of a function
f, in Cartesian and spherical coordinates respectively, is given by 

> (%Nabla = Physics:-Vectors:-Nabla)(f(X));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (97)
 

> (%Nabla = Physics:-Vectors:-Nabla)(f(S));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (98)
 

What we want now is to depart from (%Nabla = Physics:-Vectors:-Nabla)(f(X)); in Cartesian coordinates and obtain (%Nabla = Physics:-Vectors:-Nabla)(f(S)); in spherical coordinates using the transformation rule for a covariant tensor we previously computed with TransformCoordinates in Vector[column](Physics:-TensorArray(A__c[j])) = Physics:-TransformCoordinates(tr, A__s[j], [X], [S], simplifier = expand); . (An equivalent derivation which is simpler and has less steps is done in Sec. IV.) 

 

Start by changing the vector basis in the gradient (%Nabla = Physics:-Vectors:-Nabla)(f(X));  

> lhs(%Nabla(f(X)) = `+`(`*`(diff(f(X), x), `*`(_i)), `*`(diff(f(X), y), `*`(_j)), `*`(diff(f(X), z), `*`(_k)))) = Physics:-Vectors:-ChangeBasis(rhs(%Nabla(f(X)) = `+`(`*`(diff(f(X), x), `*`(_i)), `*`(d...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(99)
 

We can see that in this result the coefficients of [_r, _theta, _phi]; are the three lines in the right-hand side of Vector[column](Physics:-TensorArray(A__s[j])) = Physics:-TransformCoordinates(tr, A__c[j], [S], [X], simplifier = expand); after replacing the covariant components A__j; by the derivatives of f with respect to the jth coordinate, shown here using indexed notation due to CompactDisplay being ON 

> `~`[`=`]([A__s[1], A__s[2], A__s[3]], [Physics:-Vectors:-diff(f(S), r), Physics:-Vectors:-diff(f(S), theta), Physics:-Vectors:-diff(f(S), phi)]);
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (100)
 

> `~`[`=`]([A__c[1], A__c[2], A__c[3]], [Physics:-Vectors:-diff(f(X), x), Physics:-Vectors:-diff(f(X), y), Physics:-Vectors:-diff(f(X), z)]);
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (101)
 

So since Vector[column](Physics:-TensorArray(A__c[j])) = Physics:-TransformCoordinates(tr, A__s[j], [X], [S], simplifier = expand); is the inverse of Vector[column](Physics:-TensorArray(A__s[j])) = Physics:-TransformCoordinates(tr, A__c[j], [S], [X], simplifier = expand); , replace A by `∂`(f); in Vector[column](Physics:-TensorArray(A__c[j])) = Physics:-TransformCoordinates(tr, A__s[j], [X], [S], simplifier = expand); (the formula computed using TransformCoordinates), then insert the result in lhs(%Nabla(f(X)) = `+`(`*`(diff(f(X), x), `*`(_i)), `*`(diff(f(X), y), `*`(_j)), `*`(diff(f(X), z), `*`(_k)))) = Physics:-Vectors:-ChangeBasis(rhs(%Nabla(f(X)) = `+`(`*`(diff(f(X), x), `*`(_i)), `*`(d... to relate the gradient in Cartesian and spherical coordinates. We expect to arrive at the formula for the gradient in spherical coordinates (%Nabla = Physics:-Vectors:-Nabla)(f(S)); .  

>
 

Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-... (102)
 

>
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(103)
 

Simplifying, we arrive at (%Nabla = Physics:-Vectors:-Nabla)(f(S));  

> (lhs = `@`(`@`(expand, simplify), rhs))(%Nabla(f(X)) = `+`(`*`(`+`(`*`(`+`(`*`(sin(theta), `*`(cos(phi), `*`(diff(f(S), r)))), `/`(`*`(cos(theta), `*`(cos(phi), `*`(diff(f(S), theta)))), `*`(r)), `-`(...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (104)
 

> (lhs = `@`(simplify, rhs))(%Nabla(f(S)) = `+`(`*`(diff(f(S), r), `*`(_r)), `/`(`*`(diff(f(S), theta), `*`(_theta)), `*`(r)), `/`(`*`(diff(f(S), phi), `*`(_phi)), `*`(r, `*`(sin(theta))))));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (105)
 

IV. Deriving the transformation rule for the Divergence, Curl, Gradient and Laplacian, using Covariant derivatives 

 

  • The Divergence
 

 

Introducing the vector A in spherical coordinates, its Divergence is given by 

> A__s_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(A__r(S), _r), Physics:-`*`(`A__θ`(S), _theta)), Physics:-`*`(`A__φ`(S), _phi));
 

Typesetting:-mprintslash([A__s_ := `+`(`*`(A__r(S), `*`(_r)), `*`(`A__θ`(S), `*`(_theta)), `*`(`A__φ`(S), `*`(_phi)))], [`+`(`*`(A__r(S), `*`(_r)), `*`(`A__θ`(S), `*`(_theta)), `*`(`A_... (106)
 

> Physics:-Vectors:-CompactDisplay(%);
 

 

 

Typesetting:-mprintslash([`*`(A__r(r, theta, phi), `*`(`will now be displayed as`, `*`(A__r)))], [`*`(A__r(r, theta, phi), `*`(`will now be displayed as`, `*`(A__r)))])
Typesetting:-mprintslash([`*`(`A__φ`(r, theta, phi), `*`(`will now be displayed as`, `*`(`A__φ`)))], [`*`(`A__φ`(r, theta, phi), `*`(`will now be displayed as`, `*`(`A__φ`)))])
Typesetting:-mprintslash([`*`(`A__θ`(r, theta, phi), `*`(`will now be displayed as`, `*`(`A__θ`)))], [`*`(`A__θ`(r, theta, phi), `*`(`will now be displayed as`, `*`(`A__θ`)))]) (107)
 

> %Divergence(%A__s_) = Physics:-Vectors:-Divergence(A__s_);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (108)
 

We want to see how this result, %Divergence(%A__s_) = Physics:-Vectors:-Divergence(A__s_); , can be obtained departing from a tensorial representation of the object, this time the covariant derivative , and then using TransformCoordinates. For that purpose, first transform the coordinates and the metric introducing nonzero Christoffel symbols 

> Physics:-TransformCoordinates(tr, Physics:-g_[j, k], [S], setmetric); -1
 

 

 

 

 

_______________________________________________________;
Physics:-`*`(`Coordinates: `[r, theta, phi], `. Signature: `(`+ + +`));
_______________________________________________________;
Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(
_______________________________________________________; (109)
 

Side Note: despite having nonzero Christoffel symbols the space still has no curvature as all the components of the Riemann tensor are equal to zero 

> Physics:-Riemann[nonzero];
 

Physics:-Riemann[b, c, d, e] = {} (110)
 

Now consider the divergence of the contravariant tensor, computed in tensor notation 

> Physics:-Vectors:-CompactDisplay(A__s(S));
 

Typesetting:-mprintslash([`*`(A__s(r, theta, phi), `*`(`will now be displayed as`, `*`(A__s)))], [`*`(A__s(r, theta, phi), `*`(`will now be displayed as`, `*`(A__s)))]) (111)
 

> Physics:-D_[j](A__s[`~j`](S));
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (112)
 

Side Note: the covariant derivative expressed using the D_ operator can be rewritten in terms of the non-covariant d_ and Christoffel symbols as follows 

> Physics:-D_[j](A__s[`~j`](S), [S]) = convert(Physics:-D_[j](A__s[`~j`](S), [S]), Physics:-d_);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (113)
 

Summing over the repeated indices in Physics:-D_[j](A__s[`~j`](S)); , we have 

> %D_[j](%A__s[`~j`]) = Physics:-SumOverRepeatedIndices(Physics:-D_[j](A__s[`~j`](S), [S]));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (114)
 

How is this related to the expression of the Typesetting:-delayDotProduct(Physics:-Vectors:-Nabla, A__s_); in %Divergence(%A__s_) = Physics:-Vectors:-Divergence(A__s_); ? The answer is in the relationship established at the end of Sec I between the components of the tensor and the components of the vector A__s_; , namely that the vector components are obtained by multiplying the contravariant tensor components with the scale-factors h__j; . So in the above we need to substitute the contravariant with the vector components A__j; divided by the scale-factors 

> [seq(A__s[Physics:-Library:-Contravariant(j)](S) = Physics:-`*`(Physics:-Vectors:-Component(A__s_, j), Physics:-Vectors:-`^`(h[j], -1)), j = 1 .. 3)];
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (115)
 

> subs[eval]([A__s[`~1`](S) = A__r(S), A__s[`~2`](S) = `/`(`*`(`A__θ`(S)), `*`(r)), A__s[`~3`](S) = `/`(`*`(`A__φ`(S)), `*`(r, `*`(sin(theta))))], %D_[j](%A__s[`~j`]) = `+`(diff(A__s[`~1`](S),...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (116)
 

> (lhs = `@`(simplify, rhs))(%D_[j](%A__s[`~j`]) = `+`(diff(A__r(S), r), `/`(`*`(diff(`A__θ`(S), theta)), `*`(r)), `/`(`*`(diff(`A__φ`(S), phi)), `*`(r, `*`(sin(theta)))), `/`(`*`(2, `*`(A__r(...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (117)
 

Comparing with %Divergence(%A__s_) = Physics:-Vectors:-Divergence(A__s_); , we see these two expressions are the same: 

> (lhs = `@`(simplify, rhs))(%Divergence(%A__s_) = `+`(`/`(`*`(`+`(`*`(diff(A__r(S), r), `*`(r)), `*`(2, `*`(A__r(S))))), `*`(r)), `/`(`*`(`+`(`*`(diff(`A__θ`(S), theta), `*`(sin(theta))), `*`(`A_...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (118)
 

  • The Curl
 

 

The Curl of the the vector A__s_; in spherical coordinates is given by  

> %Curl(%A__s_) = Physics:-Vectors:-Curl(A__s_);
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (119)
 

 

One could think that the expression for the Curl in tensor notation is the same as in a non-curvilinear system 

 

But in a curvilinear system Physics:-LeviCivita[i, j, k]; is not a tensor, so instead we need to use the non-Galilean form , where is the determinant of the metric. Moreover, since the expression has one free covariant index (the first one), to compare with the vectorial formula (lhs = `@`(simplify, rhs))(%Divergence(%A__s_) = `+`(`/`(`*`(`+`(`*`(diff(A__r(S), r), `*`(r)), `*`(2, `*`(A__r(S))))), `*`(r)), `/`(`*`(`+`(`*`(diff(`A__θ`(S), theta), `*`(sin(theta))), `*`(`A_... this index also needs to be rewritten as a vector component, as discussed at the end of Sec. I, by using  

A__j = Physics:-`*`(A__j, Physics:-Vectors:-`^`(h__j, -1));  

The formula %Curl(%A__s_) = Physics:-Vectors:-Curl(A__s_); for the vectorial Curl is thus expressed using tensor notation as  

> Physics:-Vectors:-Setup(levicivita = nongalilean);
 

[levicivita = nongalilean] (120)
 

> %Curl(%A__s_) = Physics:-`*`(Physics:-LeviCivita[i, j, k], Physics:-Vectors:-`^`(%h[i], -1), Physics:-D_[`~j`](A__s[`~k`](S)));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (121)
 

after which we replace the contravariant tensor components with the vector components Physics:-`*`(A__k, Physics:-Vectors:-`^`(h__k, -1)); using [seq(A__s[Physics:-Library:-Contravariant(j)](S) = Physics:-`*`(Physics:-Vectors:-Component(A__s_, j), Physics:-Vectors:-`^`(h[j], -1)), j = 1 .. 3)]; . Proceeding the same way we did with the Divergence, expand this expression. We could use TensorArray, but Library:-TensorComponents places a comma between components which makes things more readable in this case 

> lhs(%Curl(%A__s_) = `/`(`*`(Physics:-LeviCivita[i, j, k], `*`(Physics:-D_[`~j`](A__s[`~k`](S), [S]))), `*`(%h[i]))) = Physics:-Library:-TensorComponents(rhs(%Curl(%A__s_) = `/`(`*`(Physics:-LeviCivita...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(122)
 

Now replace the components of the tensor with the components of the 3D vector A__s_; using [seq(A__s[Physics:-Library:-Contravariant(j)](S) = Physics:-`*`(Physics:-Vectors:-Component(A__s_, j), Physics:-Vectors:-`^`(h[j], -1)), j = 1 .. 3)];  

> lhs(%Curl(%A__s_) = [`/`(`*`(`+`(`*`(2, `*`(cos(theta), `*`(A__s[`~3`](S)))), `*`(diff(A__s[`~3`](S), theta), `*`(sin(theta))), `-`(`*`(csc(theta), `*`(diff(A__s[`~2`](S), phi)))))), `*`(%h[1])), `/`(...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(123)
 

> (lhs = `@`(simplify, rhs))(%Curl(%A__s_) = [`/`(`*`(`+`(`/`(`*`(2, `*`(cos(theta), `*`(`A__φ`(S)))), `*`(r, `*`(sin(theta)))), `*`(`+`(`/`(`*`(diff(`A__φ`(S), theta)), `*`(r, `*`(sin(theta))))...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (124)
 

We see that these are exactly the components of the Curl %Curl(%A__s_) = Physics:-Vectors:-Curl(A__s_);  

> (lhs = `@`(`@`(Physics:-Vectors:-Collect, simplify), rhs))(%Curl(%A__s_) = `+`(`/`(`*`(`+`(`*`(diff(`A__φ`(S), theta), `*`(sin(theta))), `*`(`A__φ`(S), `*`(cos(theta))), `-`(diff(`A__θ`(...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (125)
 

  • The Gradient
 

 

Once the problem is fully understood, it is easy to redo the computations of Sec.III for the Gradient, this time using tensor notation and the covariant derivative. In tensor notation, the components of the Gradient are given by the components of the right-hand side 

> Typesetting:-delayGradient(f(S)) = Physics:-`*`(Physics:-D_[j](f(S)), Physics:-Vectors:-`^`(h[j], -1));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi( (126)
 

where on the left-hand side we have the vectorial Nabla differential operator and on the right-hand side, since f(S); is a scalar, the covariant derivative Physics:-D_[j](f); becomes the standard derivative Physics:-d_[j](f); . 

> lhs(`+`(`*`(diff(f(S), r), `*`(_r)), `/`(`*`(diff(f(S), theta), `*`(_theta)), `*`(r)), `/`(`*`(diff(f(S), phi), `*`(_phi)), `*`(r, `*`(sin(theta))))) = `/`(`*`(Physics:-d_[j](f(S), [S])), `*`([1, r, `...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi( (127)
 

The above is the expected result (%Nabla = Physics:-Vectors:-Nabla)(f(S));  

> %Nabla(f(S)) = `+`(`*`(diff(f(S), r), `*`(_r)), `/`(`*`(diff(f(S), theta), `*`(_theta)), `*`(r)), `/`(`*`(diff(f(S), phi), `*`(_phi)), `*`(r, `*`(sin(theta)))));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (128)
 

  • The Laplacian
 

 

Likewise we can compute the Laplacian directly as 

> %Laplacian(f(S)) = Physics:-D_[j](Physics:-D_[j](f(S)));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (129)
 

In this case there are no free indices nor tensor components to be rewritten as vector components, so there is no need for scale-factors. Summing over the repeated indices, 

> Physics:-SumOverRepeatedIndices(%Laplacian(f(S)) = Physics:-D_[j](Physics:-d_[`~j`](f(S), [S]), [S]));
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (130)
 

Evaluating the  Laplacian on the left-hand side, 

> value(%Laplacian(f(S)) = `+`(Physics:-dAlembertian(f(S), [S]), `/`(`*`(2, `*`(diff(f(S), r))), `*`(r)), `/`(`*`(cot(theta), `*`(diff(f(S), theta))), `*`(`^`(r, 2)))));
 

Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msup(Typesetting:-mrow(Typesetting:-mi( (131)
 

On the right-hand side we see the dAlembertian, in curvilinear coordinates; rewrite it using standard diff derivatives and expand both sides of the equation for comparison 

> expand(convert(`/`(`*`(`+`(`*`(`^`(csc(theta), 2), `*`(diff(diff(f(S), phi), phi))), `*`(diff(diff(f(S), r), r), `*`(`^`(r, 2))), `*`(cot(theta), `*`(diff(f(S), theta))), `*`(2, `*`(diff(f(S), r), `*`...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (132)
 

This is an identity, the left and right hand sides are equal: 

> evalb(`+`(diff(diff(f(S), r), r), `/`(`*`(diff(diff(f(S), theta), theta)), `*`(`^`(r, 2))), `/`(`*`(`^`(csc(theta), 2), `*`(diff(diff(f(S), phi), phi))), `*`(`^`(r, 2))), `/`(`*`(2, `*`(diff(f(S), r))...
 

true (133)
 

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 alpha; with the horizontal plane and its projection BC over this plane forms an angle beta; with the vertical plane. The cable BC is on the same vertical plane as the bar. Determine the reactions of the planes at A and B as well as the tensions on the cables. 

Image 

Solution 

There are two equations that contain information about the state of equilibrium of a system. The first one, assuming that the sum of the forces acting on the body are equal to zero, states that the center of mass of the body is not accelerated. The second one, assuming that the sum of the moments of the forces acting on the body (that is, the total torque) is zero, states that the rotation of the body around its center of masses unchanging (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 overall problem. There is no friction so it is also clear that the reactions R_[A]; and R_[B]; are perpendicular to the planes, as shown in the figure, and the tensions T_[A]; and T_[B]; on the cables have direction AD and BC, respectively. 

 

The steps to solve this problem are: 

  • Determine each force acting on the bar and its application point
 

  • Equate the sum of the forces to zero.
 

  • Equate the sum of the moments Typesetting:-delayCrossProduct(r_, F_); to zero.
 

  • Solve these two vectorial equations for R_[A], R_[B]; , T_[A]; , and T_[B]; , 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; -1; with(Physics[Physics:-Vectors]); -1
 

The forces acting on the bar are its weight w_; and the reactions and tensions  R_[A], R_[B]; , T_[A]; , and T_[B]; . So the two equilibrium equations are: 

> eq[1] := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(w_, R_[A]), R_[B]), T_[A]), T_[B]) = 0;
 

Typesetting:-mprintslash([eq[1] := `+`(w_, R_[A], R_[B], T_[A], T_[B]) = 0], [`+`(w_, R_[A], R_[B], T_[A], T_[B]) = 0]) (134)
 

> eq[2] := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Typesetting:-delayCrossProduct(r_[w], w_), Typesetting:-delayCrossProduct(r_[A], R_[A])), Typesetting:-...
 

Typesetting:-mprintslash([eq[2] := `+`(Physics:-Vectors:-`&x`(r_[w], w_), Physics:-Vectors:-`&x`(r_[A], R_[A]), Physics:-Vectors:-`&x`(r_[B], R_[B]), Physics:-Vectors:-`&x`(r_[A], T_[A]), Physics:-Vec... (135)
 

where, in the input above, to enter the cross products you can use &x or the operator from the palette of Common Symbols. 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 R_[B]; such that the y axis in the direction of R_[A]; , and the x axis in the remaining direction, anti-parallel to T_[A]; . With these selections, the vectors entering eq[1]; and eq[2]; are projected as follows:  

> R_[B] := Physics:-`*`(abs(R[B]), _k);
 

Typesetting:-mprintslash([R_[B] := `*`(abs(R[B]), `*`(_k))], [`*`(abs(R[B]), `*`(_k))]) (136)
 

where for simplicity we use abs to represent Norm, so that abs(R[B]); represents the norm of R_[B]; to be determined, 

> r_[B] := 0;
 

Typesetting:-mprintslash([r_[B] := 0], [0]) (137)
 

 

> R_[A] := Physics:-`*`(abs(R[A]), _j);
 

Typesetting:-mprintslash([R_[A] := `*`(abs(R[A]), `*`(_j))], [`*`(abs(R[A]), `*`(_j))]) (138)
 

and where abs(R[A]); is to be determined. This reaction R_[A]; is applied to the bar at A, represented by r_[A]; ; its component along the x axis is obtained by projecting the segment BA onto the horizontal plane Physics:-`*`(L, cos(alpha)); , resulting in BC, and then into the intersection of the two planes. So: 

> r_[A] := Physics:-Vectors:-`+`(Physics:-`*`(L, cos(alpha), Physics:-Vectors:-`+`(Physics:-`*`(cos(beta), _i), Physics:-`*`(sin(Physics:-Vectors:-`+`(Physics:-`*`(2, Pi), `+`(`-`(beta)))), _j))), Physi...
 

Typesetting:-mprintslash([r_[A] := `+`(`*`(L, `*`(cos(alpha), `*`(cos(beta), `*`(_i)))), `-`(`*`(L, `*`(cos(alpha), `*`(sin(beta), `*`(_j))))), `*`(L, `*`(sin(alpha), `*`(_k))))], [`+`(`*`(L, `*`(cos(... (139)
 

For the other vectors we have 

> T_[A] := `+`(`-`(Physics:-`*`(abs(T[A]), _i)));
 

Typesetting:-mprintslash([T_[A] := `+`(`-`(`*`(abs(T[A]), `*`(_i))))], [`+`(`-`(`*`(abs(T[A]), `*`(_i))))]) (140)
 

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

Typesetting:-mprintslash([T_[B] := `+`(`*`(abs(T[B]), `*`(cos(beta), `*`(_i))), `-`(`*`(abs(T[B]), `*`(sin(beta), `*`(_j)))))], [`+`(`*`(abs(T[B]), `*`(cos(beta), `*`(_i))), `-`(`*`(abs(T[B]), `*`(sin... (141)
 

where abs(T[A]); and abs(T[B]); are to be determined, 

> w_ := `+`(`-`(Physics:-`*`(abs(w), _k)));
 

Typesetting:-mprintslash([w_ := `+`(`-`(`*`(abs(w), `*`(_k))))], [`+`(`-`(`*`(abs(w), `*`(_k))))]) (142)
 

> r_[w] := Physics:-`*`(r_[A], Physics:-Vectors:-`^`(2, -1));
 

Typesetting:-mprintslash([r_[w] := `+`(`*`(`/`(1, 2), `*`(L, `*`(cos(alpha), `*`(cos(beta), `*`(_i))))), `-`(`*`(`/`(1, 2), `*`(L, `*`(cos(alpha), `*`(sin(beta), `*`(_j)))))), `*`(`/`(1, 2), `*`(L, `*... (143)
 

The two equilibrium equations now appear as 

> eq[1];
 

`+`(`-`(`*`(abs(w), `*`(_k))), `*`(abs(R[A]), `*`(_j)), `*`(abs(R[B]), `*`(_k)), `-`(`*`(abs(T[A]), `*`(_i))), `*`(abs(T[B]), `*`(cos(beta), `*`(_i))), `-`(`*`(abs(T[B]), `*`(sin(beta), `*`(_j))))) = ... (144)
 

> 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...
`+`(`*`(`/`(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...
(145)
 

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

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

Typesetting:-mprintslash([Eq[1, 2, 3] := `+`(`-`(abs(T[A])), `*`(abs(T[B]), `*`(cos(beta)))) = 0, `+`(abs(R[A]), `-`(`*`(abs(T[B]), `*`(sin(beta))))) = 0, `+`(`-`(abs(w)), abs(R[B])) = 0], [`+`(`-`(ab... (146)
 

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

Typesetting:-mprintslash([Eq[4, 5, 6] := `+`(`*`(`/`(1, 2), `*`(L, `*`(cos(alpha), `*`(sin(beta), `*`(abs(w)))))), `-`(`*`(L, `*`(sin(alpha), `*`(abs(R[A])))))) = 0, `+`(`*`(`/`(1, 2), `*`(L, `*`(cos(...
Typesetting:-mprintslash([Eq[4, 5, 6] := `+`(`*`(`/`(1, 2), `*`(L, `*`(cos(alpha), `*`(sin(beta), `*`(abs(w)))))), `-`(`*`(L, `*`(sin(alpha), `*`(abs(R[A])))))) = 0, `+`(`*`(`/`(1, 2), `*`(L, `*`(cos(...
(147)
 

So the system of equations to be solved is 

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

Typesetting:-mprintslash([sys := {`+`(`*`(`/`(1, 2), `*`(L, `*`(cos(alpha), `*`(cos(beta), `*`(abs(w)))))), `-`(`*`(L, `*`(sin(alpha), `*`(abs(T[A])))))) = 0, `+`(`*`(L, `*`(cos(alpha), `*`(cos(beta),...
Typesetting:-mprintslash([sys := {`+`(`*`(`/`(1, 2), `*`(L, `*`(cos(alpha), `*`(cos(beta), `*`(abs(w)))))), `-`(`*`(L, `*`(sin(alpha), `*`(abs(T[A])))))) = 0, `+`(`*`(L, `*`(cos(alpha), `*`(cos(beta),...
(148)
 

The unknowns are 

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

Typesetting:-mprintslash([var := {abs(R[A]), abs(R[B]), abs(T[A]), abs(T[B])}], [{abs(R[A]), abs(R[B]), abs(T[A]), abs(T[B])}]) (149)
 

and the solution is 

> solve(sys, var);
 

{abs(R[A]) = `+`(`/`(`*`(`/`(1, 2), `*`(cos(alpha), `*`(abs(w), `*`(sin(beta))))), `*`(sin(alpha)))), abs(R[B]) = abs(w), abs(T[A]) = `+`(`/`(`*`(`/`(1, 2), `*`(cos(alpha), `*`(abs(w), `*`(cos(beta)))... (150)
 

Work 

 

Problem 

A particle is submitted to a force whose Cartesian components are given byF[z] = Physics:-`*`(c, z); . Calculate the work done by this force when moving the particle along a straight line from the origin to a point (x[0], y[0], z[0]; ). 

 

Solution 

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

 

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

 

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

  • Express x, y, and z entering the components of F_; in terms of the components of r_(t); , after which the scalar product r_; becomes Typesetting:-delayDotProduct(F_(t), Physics:-Vectors:-diff(r_(t), t)); dt and the work can be expressed as Int(Typesetting:-delayDotProduct(F_, v_), t = 0 .. t); .
 

  • Compute the integral.
 

> restart; -1; with(Physics[Physics:-Vectors]); -1
 

The work to be calculated is given by 

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

Typesetting:-mrow(Typesetting:-mi( (151)
 

The force acting on the particle is 

> F_ := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(a, Physics:-Vectors:-`^`(x, 3)), Physics:-`*`(Physics:-`*`(b, x), Physics:-Vecto...
 

Typesetting:-mprintslash([F_ := `+`(`*`(`+`(`*`(b, `*`(x, `*`(`^`(y, 2)))), `*`(42, `*`(`^`(x, 3))), `*`(c, `*`(z))), `*`(_i)), `*`(`+`(`*`(b, `*`(x, `*`(`^`(y, 2)))), `*`(42, `*`(`^`(x, 3)))), `*`(_j... (152)
 

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

> r_[0] := Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x__0, _i), Physics:-`*`(y__0, _j)), Physics:-`*`(z__0, _k));
 

Typesetting:-mprintslash([r_[0] := `+`(`*`(_i, `*`(x__0)), `*`(_j, `*`(y__0)), `*`(_k, `*`(z__0)))], [`+`(`*`(_i, `*`(x__0)), `*`(_j, `*`(y__0)), `*`(_k, `*`(z__0)))]) (153)
 

In such a case, r_(t); can be written directly as the product of the parameter t and the vector representing r_[0]; . 

> r_ := Physics:-`*`(t, r_[0]);
 

Typesetting:-mprintslash([r_ := `*`(t, `*`(`+`(`*`(_i, `*`(x__0)), `*`(_j, `*`(y__0)), `*`(_k, `*`(z__0)))))], [`*`(t, `*`(`+`(`*`(_i, `*`(x__0)), `*`(_j, `*`(y__0)), `*`(_k, `*`(z__0)))))]) (154)
 

So t = 0 at the origin, and t = 1 when r_ = r_[0]; . To construct the scalar product r_; , first express F_; 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  r_; . 

> x = Physics:-Vectors:-Component(r_, 1), y = Physics:-Vectors:-Component(r_, 2), z = Physics:-Vectors:-Component(r_, 3);
 

x = `*`(t, `*`(x__0)), y = `*`(t, `*`(y__0)), z = `*`(t, `*`(z__0)) (155)
 

> F_ := eval(F_, [x = `*`(t, `*`(x__0)), y = `*`(t, `*`(y__0)), z = `*`(t, `*`(z__0))]);
 

Typesetting:-mprintslash([F_ := `+`(`*`(`+`(`*`(b, `*`(`^`(t, 3), `*`(x__0, `*`(`^`(y__0, 2))))), `*`(42, `*`(`^`(t, 3), `*`(`^`(x__0, 3)))), `*`(c, `*`(t, `*`(z__0)))), `*`(_i)), `*`(`+`(`*`(b, `*`(`... (156)
 

Now compute v_ = Physics:-Vectors:-diff(r_(t), t); . 

> v_ := Physics:-Vectors:-diff(r_, t);
 

Typesetting:-mprintslash([v_ := `+`(`*`(_i, `*`(x__0)), `*`(_j, `*`(y__0)), `*`(_k, `*`(z__0)))], [`+`(`*`(_i, `*`(x__0)), `*`(_j, `*`(y__0)), `*`(_k, `*`(z__0)))]) (157)
 

So the line integral for the work W is now 

> W;
 

Typesetting:-mrow(Typesetting:-mi( (158)
 

When the particle is at r_[0]; we have t = 1, so the value of W to move the particle to r_[0]; is 

> eval(W, t = 1);
 

Typesetting:-mrow(Typesetting:-mi( (159)
 

> value(Int(`+`(`*`(b, `*`(`^`(t, 3), `*`(`^`(x__0, 2), `*`(`^`(y__0, 2))))), `*`(b, `*`(`^`(t, 3), `*`(x__0, `*`(`^`(y__0, 3))))), `*`(42, `*`(`^`(t, 3), `*`(`^`(x__0, 4)))), `*`(42, `*`(`^`(t, 3), `*`...
 

`+`(`*`(`/`(1, 4), `*`(b, `*`(`^`(x__0, 2), `*`(`^`(y__0, 2))))), `*`(`/`(1, 4), `*`(b, `*`(x__0, `*`(`^`(y__0, 3))))), `*`(`/`(21, 2), `*`(`^`(x__0, 4))), `*`(`/`(21, 2), `*`(`^`(x__0, 3), `*`(y__0))... (160)
 

Lagrangian for a pendulum 

 

Problem 

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

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

b) oscillates horizontally on the plane of the pendulum according to x = cos(Physics:-`*`(omega, t)); . 

 

Solution 

 

a) A figure for this part of the problem is 

Image 

The Lagrangian is defined as  

> restart; -1; with(Physics[Physics:-Vectors]); -1
 

> L := Physics:-Vectors:-`+`(T, `+`(`-`(U)));
 

Typesetting:-mprintslash([L := `+`(T, `-`(U))], [`+`(T, `-`(U))]) (161)
 

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

> U := `+`(`-`(Physics:-`*`(m, g, y)));
 

Typesetting:-mprintslash([U := `+`(`-`(`*`(m, `*`(g, `*`(y)))))], [`+`(`-`(`*`(m, `*`(g, `*`(y)))))]) (162)
 

where g is the gravitational constant and the y axis is chosen to be along the vertical, pointing downwards, resulting in the gravitational force F_[g] = Physics:-`*`(mg, _j); . The kinetic energy is: 

> T := Physics:-`*`(Physics:-Vectors:-`^`(2, -1), m, Typesetting:-delayDotProduct(v_, v_));
 

Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn( (163)
 

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

> r_ := Physics:-Vectors:-`+`(Physics:-`*`(x, _i), Physics:-`*`(y, _j));
 

Typesetting:-mprintslash([r_ := `+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)))], [`+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)))]) (164)
 

must be determined. Choosing the x axis along the horizontal 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 = Physics:-Vectors:-`+`(Physics:-`*`(a, cos(Physics:-`*`(omega, t))), Physics:-`*`(l, sin(phi(t)))), y = Physics:-Vectors:-`+`(`+`(`-`(Physics:-`*`(a, sin(Physics:-`*`(omega...
 

Typesetting:-mprintslash([parametric_equations := [x = `+`(`*`(42, `*`(cos(`*`(omega, `*`(t))))), `*`(l, `*`(sin(phi(t))))), y = `+`(`-`(`*`(42, `*`(sin(`*`(omega, `*`(t)))))), `*`(l, `*`(cos(phi(t)))... (165)
 

> r_ := eval(r_, parametric_equations);
 

Typesetting:-mprintslash([r_ := `+`(`*`(`+`(`*`(42, `*`(cos(`*`(omega, `*`(t))))), `*`(l, `*`(sin(phi(t))))), `*`(_i)), `*`(`+`(`-`(`*`(42, `*`(sin(`*`(omega, `*`(t)))))), `*`(l, `*`(cos(phi(t))))), `... (166)
 

> v_ := Physics:-Vectors:-diff(r_, t);
 

Typesetting:-mprintslash([v_ := `+`(`*`(`+`(`-`(`*`(42, `*`(omega, `*`(sin(`*`(t, `*`(omega))))))), `*`(l, `*`(diff(phi(t), t), `*`(cos(phi(t)))))), `*`(_i)), `*`(`+`(`-`(`*`(42, `*`(omega, `*`(cos(`*... (167)
 

> T;
 

Typesetting:-mprintslash([`+`(`*`(`/`(1, 2), `*`(m, `*`(`+`(`*`(`^`(`+`(`-`(`*`(42, `*`(omega, `*`(sin(`*`(t, `*`(omega))))))), `*`(l, `*`(diff(phi(t), t), `*`(cos(phi(t)))))), 2)), `*`(`^`(`+`(`-`(`*... (168)
 

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

> combine(T, trig);
 

Typesetting:-mprintslash([`+`(`*`(`/`(1, 2), `*`(m, `*`(`+`(`*`(`^`(diff(phi(t), t), 2), `*`(`^`(l, 2))), `-`(`*`(84, `*`(diff(phi(t), t), `*`(l, `*`(omega, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(phi(t))... (169)
 

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

> U := eval(U, parametric_equations);
 

Typesetting:-mprintslash([U := `+`(`-`(`*`(m, `*`(g, `*`(`+`(`-`(`*`(42, `*`(sin(`*`(t, `*`(omega)))))), `*`(l, `*`(cos(phi(t))))))))))], [`+`(`-`(`*`(m, `*`(g, `*`(`+`(`-`(`*`(42, `*`(sin(`*`(t, `*`(... (170)
 

So the requested Lagrangian is 

> L := combine(L, trig);
 

Typesetting:-mprintslash([L := `+`(`*`(`/`(1, 2), `*`(m, `*`(`+`(`*`(`^`(diff(phi(t), t), 2), `*`(`^`(l, 2))), `-`(`*`(84, `*`(diff(phi(t), t), `*`(l, `*`(omega, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(ph... (171)
 

> L := frontend(expand, [L]);
 

Typesetting:-mprintslash([L := `+`(`-`(`*`(42, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(phi(t)))), `*`(diff(phi(t), t), `*`(l, `*`(m, `*`(omega))))))), `*`(`/`(1, 2), `*`(`^`(diff(phi(t), t), 2), `*`(`^`(l... (172)
 

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 two terms arePhysics:-`*`(m, Physics:-Vectors:-`^`(a, 2), Physics:-Vectors:-`^`(omega, 2), Physics:-Vectors:-`^`(2, -1)); and so 

> L := subs(Physics:-Vectors:-`^`(omega, 2) = 0, sin(Physics:-`*`(omega, t)) = 0, L);
 

Typesetting:-mprintslash([L := `+`(`-`(`*`(42, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(phi(t)))), `*`(diff(phi(t), t), `*`(l, `*`(m, `*`(omega))))))), `*`(`/`(1, 2), `*`(`^`(diff(phi(t), t), 2), `*`(`^`(l... (173)
 

_________________________________________________________________ 

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

 

> restart; -1; with(Physics[Physics:-Vectors]); -1
 

> L := Physics:-Vectors:-`+`(T, `+`(`-`(U)));
 

Typesetting:-mprintslash([L := `+`(T, `-`(U))], [`+`(T, `-`(U))]) (174)
 

> U := `+`(`-`(Physics:-`*`(m, g, y)));
 

Typesetting:-mprintslash([U := `+`(`-`(`*`(m, `*`(g, `*`(y)))))], [`+`(`-`(`*`(m, `*`(g, `*`(y)))))]) (175)
 

> T := Physics:-`*`(Physics:-Vectors:-`^`(2, -1), m, Typesetting:-delayDotProduct(v_, v_));
 

Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn( (176)
 

> r_ := Physics:-Vectors:-`+`(Physics:-`*`(x, _i), Physics:-`*`(y, _j));
 

Typesetting:-mprintslash([r_ := `+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)))], [`+`(`*`(_i, `*`(x)), `*`(_j, `*`(y)))]) (177)
 

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

> y = Physics:-`*`(l, cos(phi(t)));
 

y = `*`(l, `*`(cos(phi(t)))) (178)
 

So the parametric equations in this case are 

> parametric_equations := [x = Physics:-Vectors:-`+`(Physics:-`*`(a, cos(Physics:-`*`(omega, t))), Physics:-`*`(l, sin(phi(t)))), y = `*`(l, `*`(cos(phi(t))))];
 

Typesetting:-mprintslash([parametric_equations := [x = `+`(`*`(42, `*`(cos(`*`(omega, `*`(t))))), `*`(l, `*`(sin(phi(t))))), y = `*`(l, `*`(cos(phi(t))))]], [[x = `+`(`*`(42, `*`(cos(`*`(omega, `*`(t)... (179)
 

> r_ := eval(r_, parametric_equations);
 

Typesetting:-mprintslash([r_ := `+`(`*`(`+`(`*`(42, `*`(cos(`*`(omega, `*`(t))))), `*`(l, `*`(sin(phi(t))))), `*`(_i)), `*`(l, `*`(cos(phi(t)), `*`(_j))))], [`+`(`*`(`+`(`*`(42, `*`(cos(`*`(omega, `*`... (180)
 

> v_ := Physics:-Vectors:-diff(r_, t);
 

Typesetting:-mprintslash([v_ := `+`(`*`(`+`(`-`(`*`(42, `*`(omega, `*`(sin(`*`(t, `*`(omega))))))), `*`(l, `*`(diff(phi(t), t), `*`(cos(phi(t)))))), `*`(_i)), `-`(`*`(l, `*`(diff(phi(t), t), `*`(sin(p... (181)
 

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

> U := eval(U, parametric_equations);
 

Typesetting:-mprintslash([U := `+`(`-`(`*`(m, `*`(g, `*`(l, `*`(cos(phi(t))))))))], [`+`(`-`(`*`(m, `*`(g, `*`(l, `*`(cos(phi(t))))))))]) (182)
 

So the requested Lagrangian is 

> L;
 

Typesetting:-mprintslash([`+`(`*`(`/`(1, 2), `*`(m, `*`(`+`(`*`(`^`(`+`(`-`(`*`(42, `*`(omega, `*`(sin(`*`(t, `*`(omega))))))), `*`(l, `*`(diff(phi(t), t), `*`(cos(phi(t)))))), 2)), `*`(`^`(l, 2), `*`... (183)
 

> L := combine(L, trig);
 

Typesetting:-mprintslash([L := `+`(`*`(`/`(1, 2), `*`(m, `*`(`+`(`*`(`^`(diff(phi(t), t), 2), `*`(`^`(l, 2))), `-`(`*`(42, `*`(diff(phi(t), t), `*`(l, `*`(omega, `*`(sin(`+`(`*`(t, `*`(omega)), phi(t)... (184)
 

> L := frontend(expand, [L]);
 

Typesetting:-mprintslash([L := `+`(`-`(`*`(21, `*`(sin(`+`(`*`(t, `*`(omega)), phi(t))), `*`(diff(phi(t), t), `*`(l, `*`(m, `*`(omega))))))), `-`(`*`(21, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(phi(t)))),...
Typesetting:-mprintslash([L := `+`(`-`(`*`(21, `*`(sin(`+`(`*`(t, `*`(omega)), phi(t))), `*`(diff(phi(t), t), `*`(l, `*`(m, `*`(omega))))))), `-`(`*`(21, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(phi(t)))),...
(185)
 

The terms in L that can be expressed as total derivatives can be discarded are those involving Physics:-Vectors:-`^`(omega, 2); and cos(Physics:-`*`(Physics:-`*`(2, omega), t)); , so the Lagrangian is 

> L := subs(Physics:-Vectors:-`^`(omega, 2) = 0, cos(Physics:-`*`(Physics:-`*`(2, omega), t)) = 0, L);
 

Typesetting:-mprintslash([L := `+`(`-`(`*`(21, `*`(sin(`+`(`*`(t, `*`(omega)), phi(t))), `*`(diff(phi(t), t), `*`(l, `*`(m, `*`(omega))))))), `-`(`*`(21, `*`(sin(`+`(`*`(t, `*`(omega)), `-`(phi(t)))),... (186)
 

 

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 m[1]; in the extremes of the base and a mass m[2]; in the third vertex. The distance between the two masses m[1]; is equal to a, and the height of the triangle is equal to h. 

 

Solution 

> restart; -1; with(Physics, Physics:-KroneckerDelta); -1; with(Physics:-Vectors); -1
 

The general formula for the Inertia tensor is given by (note the use of the abbreviation kd_ for KroneckerDelta) 

> InertiaTensor := %sum(Physics:-`*`(m[k], Physics:-Vectors:-`+`(Physics:-`*`(Physics:-Vectors:-`^`(Physics:-Vectors:-Norm(r_[k]), 2), kd_[i, j]), `+`(`-`(Physics:-`*`(Physics:-Vectors:-Component(r_[k],...
 

Typesetting:-mrow(Typesetting:-mi( (187)
 

where N is the number of particles, m[k]; is the mass of each particle, r_[k]; is its position in a reference system with the origin at the "center of mass", r_[k][i]; is the component in the ith direction of the position vector associated to the kth particle in the reference system, and delta[i, j]; is the Kronecker delta, part of Physics. Set this definition of the InertiaTensor to be an indexing function for the InertiaTensor matrix. 

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

Typesetting:-mprintslash([IT := proc (i, j) options operator, arrow; %sum(`*`(m[k], `*`(`+`(`*`(`^`(Physics:-Vectors:-Norm(r_[k]), 2), `*`(Physics:-KroneckerDelta[i, j])), `-`(`*`(Physics:-Vectors:-Co... (188)
 

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

> IT(1, 1);
 

Typesetting:-mrow(Typesetting:-mi( (189)
 

and outside of the diagonal we have 

> IT(1, 2);
 

Typesetting:-mrow(Typesetting:-mi( (190)
 

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

> N := 3;
 

Typesetting:-mprintslash([N := 3], [3]) (191)
 

Hence the matrix is 

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

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mi( (192)
 

Two of the masses are equal 

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

Typesetting:-mprintslash([m[3] := m[1]], [m[1]]) (193)
 

Now choose any system of reference (not at the center of mass) where we are going to project the position vectors R_[k]; of each atom as well as the center of mass R_[CM]; . The vectors r_[k]; entering the definition of the Inertia tensor InertiaTensor := %sum(Physics:-`*`(m[k], Physics:-Vectors:-`+`(Physics:-`*`(Physics:-Vectors:-`^`(Physics:-Vectors:-Norm(r_[k]), 2), kd_[i, j]), `+`(`-`(Physics:-`*`(Physics:-Vectors:-Component(r_[k],... are related to R_[k]; and R_[CM]; by 

> position := r_[k] = Physics:-Vectors:-`+`(R_[k], `+`(`-`(R_[CM])));
 

Typesetting:-mprintslash([position := r_[k] = `+`(R_[k], `-`(R_[CM]))], [r_[k] = `+`(R_[k], `-`(R_[CM]))]) (194)
 

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

> R_[1] := `+`(`-`(Physics:-`*`(a, Physics:-Vectors:-`^`(2, -1), _i)));
 

Typesetting:-mprintslash([R_[1] := `+`(`-`(`*`(21, `*`(_i))))], [`+`(`-`(`*`(21, `*`(_i))))]) (195)
 

> R_[2] := Physics:-`*`(h, _k);
 

Typesetting:-mprintslash([R_[2] := `*`(h, `*`(_k))], [`*`(h, `*`(_k))]) (196)
 

> R_[3] := Physics:-`*`(a, Physics:-Vectors:-`^`(2, -1), _i);
 

Typesetting:-mprintslash([R_[3] := `+`(`*`(21, `*`(_i)))], [`+`(`*`(21, `*`(_i)))]) (197)
 

Indicate the real objects of this problem so that simplification steps further below can take that into account 

> Physics:-Vectors:-Setup(real = {a, h, m[1], m[2], m[3]});
 

 

Physics:-`*`(`* Partial match of '`, `real' against keyword '`, `realobjects' `);
Error, (in Physics:-Setup) invalid input: Physics:-Setup expects value for keyword parameter realobjects to be of type {"not given", Physics:-name, function, list({Physics:-name, function}), set({Physics:-name, function})}, but received {42, h, m[1], m[2]}
 

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

> R_[CM] := Physics:-`*`(%sum(Physics:-`*`(m[k], R_[k]), k = 1 .. N), Physics:-Vectors:-`^`(%sum(m[k], k = 1 .. N), -1));
 

Typesetting:-mprintslash([R_[CM] := `/`(`*`(%sum(`*`(m[k], `*`(R_[k])), k = 1 .. 3)), `*`(%sum(m[k], k = 1 .. 3)))], [`/`(`*`(%sum(`*`(m[k], `*`(R_[k])), k = 1 .. 3)), `*`(%sum(m[k], k = 1 .. 3)))]) (198)
 

Evaluating these sums, we have the value of R_[CM]; : 

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

Typesetting:-mprintslash([R_[CM] := `/`(`*`(m[2], `*`(h, `*`(_k))), `*`(`+`(`*`(2, `*`(m[1])), m[2])))], [`/`(`*`(m[2], `*`(h, `*`(_k))), `*`(`+`(`*`(2, `*`(m[1])), m[2])))]) (199)
 

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

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

r_[1] = `+`(`-`(`*`(21, `*`(_i))), `-`(`/`(`*`(m[2], `*`(h, `*`(_k))), `*`(`+`(`*`(2, `*`(m[1])), m[2]))))), r_[2] = `+`(`*`(h, `*`(_k)), `-`(`/`(`*`(m[2], `*`(h, `*`(_k))), `*`(`+`(`*`(2, `*`(m[1])),... (200)
 

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

> IT_answer := simplify(eval(value(IT_Matrix), [r_[1] = `+`(`-`(`*`(21, `*`(_i))), `-`(`/`(`*`(m[2], `*`(h, `*`(_k))), `*`(`+`(`*`(2, `*`(m[1])), m[2]))))), r_[2] = `+`(`*`(h, `*`(_k)), `-`(`/`(`*`(m[2]...
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mfrac(Typesetting:-mrow(Typesetting:-mn( (201)
 

The equations of motion in curvilinear coordinates, tensor notation and Coriolis force 

 

The formulation of the equations of motion of a particle in Cartesian coordinates using vector notation is simple. However, depending on the problem, it may be advantageous to use curvilinear coordinates and also tensor notation, for example when describing the motion of a particle as seen from a non-inertial system of reference (e.g. a rotating planet, like earth). When the particle's movement is observed from one such rotating referential, we also see "acceleration" that is not due to any force but to the fact that the referential itself is accelerated. The material below discusses and formulates these topics, and derives the expression for the Coriolis and centripetal force in cylindrical coordinates as seen from a rotating system of reference. 

 

Cartesian coordinates 

 

Vector notation 

 

Generally speaking, the equations of motion of a particle are easy to formulate: the position vector is a function of time, the velocity is its first derivative, and the acceleration is its second derivative. For instance, in Cartesian coordinates 

> restart; 1; with(Physics); -1; with(Physics:-Vectors); -1
 

> r_(t) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x(t), _i), Physics:-`*`(y(t), _j)), Physics:-`*`(z(t), _k));
 

r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k))) (202)
 

> Physics:-Vectors:-diff(r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k))), t);
 

Typesetting:-mprintslash([diff(r_(t), t) = `+`(`*`(diff(x(t), t), `*`(_i)), `*`(diff(y(t), t), `*`(_j)), `*`(diff(z(t), t), `*`(_k)))], [diff(r_(t), t) = `+`(`*`(diff(x(t), t), `*`(_i)), `*`(diff(y(t)... (203)
 

> Physics:-Vectors:-diff(diff(r_(t), t) = `+`(`*`(diff(x(t), t), `*`(_i)), `*`(diff(y(t), t), `*`(_j)), `*`(diff(z(t), t), `*`(_k))), t);
 

Typesetting:-mprintslash([diff(r_(t), `$`(t, 2)) = `+`(`*`(diff(x(t), `$`(t, 2)), `*`(_i)), `*`(diff(y(t), `$`(t, 2)), `*`(_j)), `*`(diff(z(t), `$`(t, 2)), `*`(_k)))], [diff(diff(r_(t), t), t) = `+`(`... (204)
 

Newton's 2nd law, that in an inertial system of reference when there is force there is acceleration and vice versa, is then given by 

> F_(t) = Physics:-`*`(m, lhs(diff(diff(r_(t), t), t) = `+`(`*`(diff(diff(x(t), t), t), `*`(_i)), `*`(diff(diff(y(t), t), t), `*`(_j)), `*`(diff(diff(z(t), t), t), `*`(_k)))));
 

Typesetting:-mprintslash([F_(t) = `*`(m, `*`(diff(r_(t), `$`(t, 2))))], [F_(t) = `*`(m, `*`(diff(diff(r_(t), t), t)))]) (205)
 

where F_(t) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(F__x(t), _i), Physics:-`*`(F__y(t), _j)), Physics:-`*`(F__z(t), _k)); represents the total force acting on the particle. This vectorial equation represents three second order differential equations which, given initial conditions, can be integrated to arrive at a closed form expression for r_(t); . 

 

Tensor notation 

 

In Cartesian coordinates, the tensorial form of the equations F_(t) = Physics:-`*`(m, lhs(diff(diff(r_(t), t), t) = `+`(`*`(diff(diff(x(t), t), t), `*`(_i)), `*`(diff(diff(y(t), t), t), `*`(_j)), `*`(diff(diff(z(t), t), t), `*`(_k))))); is also straightforward. In a flat spacetime - Galilean system of reference - the three space coordinates x, y, z; form a 3D tensor, and so do its first and second derivatives. Set the spacetime to be 3-dimensional and Euclidean, and use lowercaselatin indices for the corresponding tensors 

> Physics:-Vectors:-Setup(coordinates = cartesian, metric = Euclidean, dimension = 3, spacetimeindices = lowercaselatin); -1
 

 

 

 

 

 

 

Physics:-`*`(`Systems of spacetime coordinates are:`, {X = (x, y, z)});
_______________________________________________________;
Physics:-`*`(`The Euclidean metric in coordinates `, [x, y, z]);
_______________________________________________________;
Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(
_______________________________________________________; (206)
 

The position, velocity, and acceleration vectors are expressed in tensor notation as in r_(t) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x(t), _i), Physics:-`*`(y(t), _j)), Physics:-`*`(z(t), _k)); , Physics:-Vectors:-diff(r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k))), t); and Physics:-Vectors:-diff(diff(r_(t), t) = `+`(`*`(diff(x(t), t), `*`(_i)), `*`(diff(y(t), t), `*`(_j)), `*`(diff(z(t), t), `*`(_k))), t);  

> X[j](t);
 

(X)[j](t) (207)
 

> Physics:-Vectors:-diff((X)[j](t), t);
 

Typesetting:-mprintslash([diff(X[j](t), t)], [Physics:-Vectors:-diff((Physics:-SpaceTimeVector[j](X))(t), t)]) (208)
 

> Physics:-Vectors:-diff(Physics:-Vectors:-diff((Physics:-SpaceTimeVector[j](X))(t), t), t);
 

Typesetting:-mprintslash([diff(X[j](t), `$`(t, 2))], [Physics:-Vectors:-diff(Physics:-Vectors:-diff((Physics:-SpaceTimeVector[j](X))(t), t), t)]) (209)
 

Setting a tensor F__j(t); to represent the three Cartesian components of the force 

> Physics:-Define(F[j] = [F__x(t), F__y(t), F__z(t)]);
 

 

`Defined objects with tensor properties`;
Typesetting:-mprintslash([{Physics:-Dgamma[b], F[j], Physics:-Psigma[b], Physics:-d_[b], Physics:-g_[b, c], Physics:-LeviCivita[b, c, d], X[b]}], [{Physics:-Dgamma[b], F[j], Physics:-Psigma[b], Physic... (210)
 

Newton's 2nd law F_(t) = Physics:-`*`(m, lhs(diff(diff(r_(t), t), t) = `+`(`*`(diff(diff(x(t), t), t), `*`(_i)), `*`(diff(diff(y(t), t), t), `*`(_j)), `*`(diff(diff(z(t), t), t), `*`(_k))))); , now expressed in tensorial notation, is given by 

> F[j] = Physics:-`*`(m, Physics:-Vectors:-diff(Physics:-Vectors:-diff((Physics:-SpaceTimeVector[j](X))(t), t), t));
 

Typesetting:-mprintslash([F[j] = `*`(m, `*`(diff(X[j](t), `$`(t, 2))))], [F[j] = `*`(m, `*`(diff(diff((Physics:-SpaceTimeVector[j](X))(t), t), t)))]) (211)
 

The three differential equations represented by this tensorial form of F_(t) = Physics:-`*`(m, lhs(diff(diff(r_(t), t), t) = `+`(`*`(diff(diff(x(t), t), t), `*`(_i)), `*`(diff(diff(y(t), t), t), `*`(_j)), `*`(diff(diff(z(t), t), t), `*`(_k))))); are as expected 

> Physics:-TensorArray(F[j] = `*`(m, `*`(diff(diff((Physics:-SpaceTimeVector[j](X))(t), t), t))), output = setofequations);
 

Typesetting:-mprintslash([{F__x(t) = `*`(m, `*`(diff(x(t), `$`(t, 2)))), F__y(t) = `*`(m, `*`(diff(y(t), `$`(t, 2)))), F__z(t) = `*`(m, `*`(diff(z(t), `$`(t, 2))))}], [{F__x(t) = `*`(m, `*`(diff(diff(... (212)
 

Things are straightforward in Cartesian coordinates because the components of the line element dr_ = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(dx, _i), Physics:-`*`(dy, _j)), Physics:-`*`(dz, _k)); are exactly the components of the tensor Physics:-d_(X[j]);  

> Physics:-TensorArray(Physics:-d_(X[j]));
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mi( (213)
 

and so the form factors (see related MaplePrimes post) are all equal to 1. 

 

In the case of no external forces `and`(F_(t) = 0, 0 = F[j]); , and the equations of motion, whose solution is the trajectory, can be formulated as the path of minimal length between two points, a geodesic. In this case, as the spacetime is flat the geometry is Euclidean, as can be seen in the 3x3 identity metric matrix Physics:-Vectors:-Setup(coordinates = cartesian, metric = Euclidean, dimension = 3, spacetimeindices = lowercaselatin); -1; the previously mentioned two points lie on a plane, and the geodesic is a straight line. The differential equations of this geodesic are thus the equations of motion Physics:-TensorArray(F[j] = `*`(m, `*`(diff(diff((Physics:-SpaceTimeVector[j](X))(t), t), t))), output = setofequations); with  F_(t) = 0; , and can be computed using Geodesics 

> Physics:-Geodesics(t);
 

Typesetting:-mprintslash([[diff(z(t), `$`(t, 2)) = 0, diff(y(t), `$`(t, 2)) = 0, diff(x(t), `$`(t, 2)) = 0]], [[diff(diff(z(t), t), t) = 0, diff(diff(y(t), t), t) = 0, diff(diff(x(t), t), t) = 0]]) (214)
 

Curvilinear coordinates 

 

Vector notation 

 

In the case of curvilinear coordinates, for example with cylindrical or spherical variables, the form of these equations is obtained by performing a change of coordinates. 

> tr := `~`[`=`]([X], ` $`, Physics:-Vectors:-ChangeCoordinates([X], cylindrical));
 

Typesetting:-mprintslash([tr := [x = `*`(rho, `*`(cos(phi))), y = `*`(rho, `*`(sin(phi))), z = z]], [[x = `*`(rho, `*`(cos(phi))), y = `*`(rho, `*`(sin(phi))), z = z]]) (215)
 

This change of coordinates leaves the z axis and its corresponding unit vector _k; unchanged. By changing the basis and coordinates used to represent the position vector r_(t) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x(t), _i), Physics:-`*`(y(t), _j)), Physics:-`*`(z(t), _k)); , the equation becomes 

> r_(t) = Physics:-Vectors:-ChangeBasis(rhs(r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k)))), cylindrical, alsocomponents);
 

r_(t) = `+`(`*`(z(t), `*`(_k)), `*`(rho(t), `*`(_rho(t)))) (216)
 

where, since in r_(t) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(x(t), _i), Physics:-`*`(y(t), _j)), Physics:-`*`(z(t), _k)); the coordinates x, y, z; depend on t, in r_(t) = Physics:-Vectors:-ChangeBasis(rhs(r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k)))), cylindrical, alsocomponents); not just rho(t); and z(t); but also the unit vector _rho(t); depend on t. The velocity is computed as usual, by differentiating 

> Physics:-Vectors:-diff(r_(t) = `+`(`*`(z(t), `*`(_k)), `*`(rho(t), `*`(_rho(t)))), t);
 

Typesetting:-mprintslash([diff(r_(t), t) = `+`(`*`(diff(z(t), t), `*`(_k)), `*`(diff(rho(t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(phi(t), t), `*`(_phi(t)))))], [diff(r_(t), t) = `+`(`*`(diff(z(t), ... (217)
 

The second term is due to the dependency of _rho; on the coordinate phi; as well as the chain rule `and`(Physics:-Vectors:-diff(_rho(t), t) = Physics:-`*`(Physics:-Vectors:-diff(_rho, phi), Physics:-Vectors:-diff(phi(t), t)), Physics:-`*`(Physics:-Vectors:-diff(_rho, phi), Physics:-Vectors:-diff(ph.... The dependency of curvilinear unit vectors on the coordinates is automatically taken into account when differentiating due to the Setup setting geometricdifferentiation = true.  

 

For the acceleration, 

> Physics:-Vectors:-diff(diff(r_(t), t) = `+`(`*`(diff(z(t), t), `*`(_k)), `*`(diff(rho(t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(phi(t), t), `*`(_phi(t))))), t);
 

Typesetting:-mprintslash([diff(r_(t), `$`(t, 2)) = `+`(`*`(`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi... (218)
 

where the term involving Physics:-Vectors:-`^`(Physics:-Vectors:-diff(phi(t), t), 2); comes from differentiating _phi(t); in Physics:-Vectors:-diff(r_(t) = `+`(`*`(z(t), `*`(_k)), `*`(rho(t), `*`(_rho(t)))), t); while taking into account the dependency of _phi; on the coordinate This result can also be obtained by directly changing the variables in the acceleration Physics:-Vectors:-diff(r_(t), t, t); in Physics:-Vectors:-diff(diff(r_(t), t) = `+`(`*`(diff(x(t), t), `*`(_i)), `*`(diff(y(t), t), `*`(_j)), `*`(diff(z(t), t), `*`(_k))), t);  

> lhs(diff(diff(r_(t), t), t) = `+`(`*`(diff(diff(x(t), t), t), `*`(_i)), `*`(diff(diff(y(t), t), t), `*`(_j)), `*`(diff(diff(z(t), t), t), `*`(_k)))) = Physics:-Vectors:-ChangeBasis(rhs(diff(diff(r_(t)...
 

Typesetting:-mprintslash([diff(r_(t), `$`(t, 2)) = `+`(`*`(`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi... (219)
 

 

Newton's 2nd law becomes 

> F_(t) = Physics:-`*`(m, rhs(diff(diff(r_(t), t), t) = `+`(`*`(`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff...
 

Typesetting:-mprintslash([F_(t) = `*`(m, `*`(`+`(`*`(`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t... (220)
 

In the absence of external forces, by equating the vector components to 0, that is the coefficients of the unit vectors of the acceleration Physics:-Vectors:-diff(r_(t), t, t); , we get the system of differential equations for rho(t), phi(t), z(t); whose solution is the trajectory of the particle in cylindrical coordinates 

> `~`[`=`](Typesetting:-delayDotProduct([_rho(t), _phi(t), _k], rhs(diff(diff(r_(t), t), t) = `+`(`*`(`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+...
 

Typesetting:-mprintslash([[`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))) = 0, `+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), `*`(diff(phi(t), `$`(t,... (221)
 

> solve([`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))) = 0, `+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), `*`(diff(diff(phi(t), t), t)))) = 0, diff(...
 

Typesetting:-mprintslash([{diff(phi(t), `$`(t, 2)) = `+`(`-`(`/`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t))))), diff(rho(t), `$`(t, 2)) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2)))... (222)
 

In this formulation solve([`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))) = 0, `+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), `*`(diff(diff(phi(t), t), t)))) = 0, diff(... with F_(t) = 0; , although Physics:-Vectors:-diff(z(t), t, t) = 0; (no acceleration in the _k; direction) is naturally expected, the same cannot be said for the other two equations for Physics:-Vectors:-diff(phi(t), t, t); and Physics:-Vectors:-diff(rho(t), t, t); , as discussed below under Coriolis and Centripetal forces. The key observation at this point is that the right-hand sides of both unexpected equations involve Physics:-Vectors:-diff(phi(t), t); , rotation around the z axis. 

 

Tensor notation 

 

Both F_(t) = Physics:-`*`(m, rhs(diff(diff(r_(t), t), t) = `+`(`*`(`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff... and solve([`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))) = 0, `+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), `*`(diff(diff(phi(t), t), t)))) = 0, diff(... return the same equations when using tensor notation. When doing so, one can transform the position (168), velocity (169), and acceleration (170) tensors, but since these tensors are expressed as functions of a parameter (the time), it is simpler to only transform the underlying metric via TransformCoordinates. This has the advantage that all the geometrical subtleties of curvilinear coordinates, like scale factors and dependency of unit vectors on curvilinear coordinates, are automatically and succinctly encoded in the metric: 

> Physics:-TransformCoordinates(tr, Physics:-g_[j, k], [X], setmetric); -1
 

 

 

 

 

_______________________________________________________;
Physics:-`*`(`Coordinates: `[X], `. Signature: `(`+ + +`));
_______________________________________________________;
Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(
_______________________________________________________; (223)
 

The computation of geodesics assumes that the coordinates X; depend on a parameter. That parameter is passed as the first argument to Geodesics 

> Physics:-Geodesics(t);
 

Typesetting:-mprintslash([[diff(rho(t), `$`(t, 2)) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2))), diff(phi(t), `$`(t, 2)) = `+`(`-`(`/`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t)))))... (224)
 

These equations of motion are the same as the equations in solve([`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))) = 0, `+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), `*`(diff(diff(phi(t), t), t)))) = 0, diff(... which were computed using standard vector notation while differentiating and taking into account the dependency of curvilinear unit vectors on the curvilinear coordinates in Physics:-Vectors:-diff(r_(t) = `+`(`*`(z(t), `*`(_k)), `*`(rho(t), `*`(_rho(t)))), t); and Physics:-Vectors:-diff(diff(r_(t), t) = `+`(`*`(diff(z(t), t), `*`(_k)), `*`(diff(rho(t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(phi(t), t), `*`(_phi(t))))), t); .  One of the interesting features of computing with tensors is that, as previously mentioned, the geometrical subtleties of curvilinear coordinates are automatically encoded in the metric Physics:-TransformCoordinates(tr, Physics:-g_[j, k], [X], setmetric); -1.  

 

To understand how the geodesic equations are computed in one go in Physics:-Geodesics(t); , one can perform the calculation in steps:  

  • Make rho; a function of t; directly in the metric
 

  • Compute - not the final form of the equations Physics:-Geodesics(t); - but the intermediate form which expresses the geodesic equation using tensor notation, in terms of Christoffel symbols
 

  • Compute the components of this tensorial equation for the geodesic (using TensorArray)
 

 

For step 1, we have 

> subs(rho = rho(t), Physics:-g_[]);
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (225)
 

Set this metric where `≡`(rho, rho(t));  

> Physics:-Vectors:-Setup(Physics:-g_[b, c] = Matrix(%id = 36893628832830130164)); -1
 

 

 

 

 

_______________________________________________________;
Physics:-`*`(`Coordinates: `[X], `. Signature: `(`+ + +`));
_______________________________________________________;
Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi(
_______________________________________________________; (226)
 

For step 2, the geodesic equations in tensor notation with the coordinates depending on the time t are computed by passing the optional argument tensornotation 

> Physics:-Geodesics(t, tensornotation);
 

Typesetting:-mprintslash([`+`(diff(X[`~b`](t), `$`(t, 2)), `*`(Physics:-Christoffel[`~b`, c, d], `*`(diff(X[`~c`](t), t), `*`(diff(X[`~d`](t), t))))) = 0], [`+`(diff(diff((Physics:-SpaceTimeVector[`~b... (227)
 

Step 3: compute the components of this tensorial equation 

> Physics:-TensorArray(`+`(diff(diff((Physics:-SpaceTimeVector[`~b`](X))(t), t), t), `*`(Physics:-Christoffel[`~b`, c, d], `*`(diff((Physics:-SpaceTimeVector[`~c`](X))(t), t), `*`(diff((Physics:-SpaceTi...
 

Typesetting:-mprintslash([[`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))) = 0, `+`(diff(phi(t), `$`(t, 2)), `/`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rh... (228)
 

This is the expected result, the same as Physics:-Geodesics(t); .  

 

Having the tensorial equation Physics:-Geodesics(t, tensornotation); is also useful for formulating the equations of motion in the presence of force. For this purpose, redefine the contravariant tensor Physics:-Vectors:-`^`(F, j); to represent the force in the cylindrical basis 

> Physics:-Define(F[`~j`] = [`F__ρ`(t), `F__φ`(t), F__z(t)]);
 

 

`Defined objects with tensor properties`;
Typesetting:-mprintslash([{Physics:-D_[b], Physics:-Dgamma[b], F[`~j`], Physics:-Psigma[b], Physics:-Ricci[b, c], Physics:-Riemann[b, c, d, e], Physics:-Weyl[b, c, d, e], Physics:-d_[b], Physics:-g_[b... (229)
 

 

Newton's 2nd law F_(t) = Physics:-`*`(m, rhs(diff(diff(r_(t), t), t) = `+`(`*`(`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff... 

> F_(t) = `*`(m, `*`(`+`(`*`(`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), `*`(dif...
 

Typesetting:-mprintslash([F_(t) = `*`(m, `*`(`+`(`*`(`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t... (230)
 

now using tensorial notation, becomes 

> F[`~a`] = Physics:-`*`(m, lhs(`+`(diff(diff((Physics:-SpaceTimeVector[`~b`](X))(t), t), t), `*`(Physics:-Christoffel[`~b`, c, d], `*`(diff((Physics:-SpaceTimeVector[`~c`](X))(t), t), `*`(diff((Physics...
 

Typesetting:-mprintslash([F[`~a`] = `*`(m, `*`(`+`(diff(X[`~b`](t), `$`(t, 2)), `*`(Physics:-Christoffel[`~b`, c, d], `*`(diff(X[`~c`](t), t), `*`(diff(X[`~d`](t), t)))))))], [F[`~a`] = `*`(m, `*`(`+`... (231)
 

> Physics:-TensorArray(F[`~a`] = `*`(m, `*`(`+`(diff(diff((Physics:-SpaceTimeVector[`~b`](X))(t), t), t), `*`(Physics:-Christoffel[`~b`, c, d], `*`(diff((Physics:-SpaceTimeVector[`~c`](X))(t), t), `*`(d...
 

Error, (in Physics:-TensorArray) free indices on both sides of the equation are different: found {~a} on the left-hand side and {~b} on the right-hand side
 

where we recall (see related MaplePrimes post) that, to obtain the vector components entering F_(t); from these tensor components F[`~a`]; , we need to multiply the latter by the scale factors 1, rho(t), 1; , and so the component of F_(t); in the direction of _phi; is equal to Physics:-`*`(rho(t), m, Physics:-Vectors:-`+`(Physics:-Vectors:-diff(phi(t), t, t), Physics:-`*`(2, Physics:-Vectors:-diff(rho(t), t), Physics:-Vectors:-diff(phi(t), t), Physics:-Vectors:-`^`(rho(t), .... 

Coriolis force and centripetal force 

 

After changing variables the position vector r_(t); in r_(t) = Physics:-Vectors:-ChangeBasis(rhs(r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k)))), cylindrical, alsocomponents); is expressed as 

r_(t) = Physics:-Vectors:-`+`(Physics:-`*`(z(t), _k), Physics:-`*`(_rho(t), rho(t)));  

A distinction needs to be made here as to whether or not the unit vector _rho; depends on the time t; , with the former being the general case. When _rho; is a constant, the value of the coordinate phi; (the angle between _rho; and the x axis) does not change; there is no rotation around the z axis. On the other hand, when `≡`(_rho, _rho(t)); , the orientation of _rho; , and thus the coordinate phi; , changes with time, meaning that either the force F_(t); acting on the particle has a component in the _phi; direction which produces rotation around the z axis or the system of reference itself is rotating around the z axis. 

 

Likewise, the expression r_(t) = Physics:-Vectors:-ChangeBasis(rhs(r_(t) = `+`(`*`(x(t), `*`(_i)), `*`(y(t), `*`(_j)), `*`(z(t), `*`(_k)))), cylindrical, alsocomponents);  can also represent the position vector measured in the original Galilean (inertial) system of reference, where a force F_(t); is producing rotation around the z axis, as well as the position of the particle measured in a rotating non-inertial system of reference. It follows that the transformation tr := `~`[`=`]([X], ` $`, Physics:-Vectors:-ChangeCoordinates([X], cylindrical)); can also be interpreted in two different ways: representing the position of the particle in the original inertial system of reference, or representing a transformation from an inertial system of reference to a rotating non-inertial system of reference. 

 

This equivalence between the trajectory of a particle subject to an external force as observed in an inertial system of reference, and the same trajectory with no external force observed from a non-inertial accelerated system of reference, is actually at the root of the formulation of general relativity, and is also well known in classical mechanics. The (apparent) forces observed in the rotating non-inertial system of reference, due only to its acceleration, are called Coriolis and centripetal forces. 

 

The equations  

Physics:-Vectors:-diff(rho(t), t, t) = Physics:-`*`(Physics:-Vectors:-`^`(Physics:-Vectors:-diff(phi(t), t), 2), rho(t)); ,  Physics:-Vectors:-diff(phi(t), t, t) = `+`(`-`(Physics:-`*`(2, Physics:-Vectors:-diff(phi(t), t), Physics:-Vectors:-diff(rho(t), t), Physics:-Vectors:-`^`(rho(t), -1))));  

that appeared in Physics:-TensorArray(`+`(diff(diff((Physics:-SpaceTimeVector[`~b`](X))(t), t), t), `*`(Physics:-Christoffel[`~b`, c, d], `*`(diff((Physics:-SpaceTimeVector[`~c`](X))(t), t), `*`(diff((Physics:-SpaceTi... when in the inertial system of reference `and`(F_(t) = Physics:-`*`(m, Physics:-Vectors:-diff(r_(t), t, t)), Physics:-`*`(m, Physics:-Vectors:-diff(r_(t), t, t)) = 0); , are related to the Coriolis and centripetal forces in the non-inertial referential. To see this relation, following the presentation in [1], introduce a vector omega_; representing the rotation of the non-inertial referential around the z axis. When the non-inertial system rotates clockwise in the inertial system of reference, in the non-inertial system phi; increases in value in the anti-clockwise direction. 

> omega_ = `+`(`-`(Physics:-`*`(Physics:-Vectors:-diff(phi(t), t), _k)));
 

Typesetting:-mprintslash([omega_ = `+`(`-`(`*`(diff(phi(t), t), `*`(_k))))], [omega_ = `+`(`-`(`*`(diff(phi(t), t), `*`(_k))))]) (232)
 

According to [1], (39.7), the acceleration Physics:-Vectors:-diff(r_(t), t, t); in the inertial system is expressed in terms of the quantities in the non-inertial rotating system by the sum of the following three vectorial terms. First, the components of the acceleration a_(t); measured in the non-inertial system, which are given by the second derivatives of the coordinates rho(t), phi(t), z(t); multiplied by the scale factors, which in this case are given by 1, rho(t), 1; (see this post in MaplePrimes) 

> a_(t) = Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(Physics:-`*`(Physics:-Vectors:-diff(rho(t), t, t), _rho(t)), Physics:-`*`(rho(t), Physics:-Vectors:-diff(phi(t), t, t), _phi(t))), Physics:-`*`(Phys...
 

Typesetting:-mprintslash([a_(t) = `+`(`*`(diff(rho(t), `$`(t, 2)), `*`(_rho(t))), `*`(rho(t), `*`(diff(phi(t), `$`(t, 2)), `*`(_phi(t)))), `*`(diff(z(t), `$`(t, 2)), `*`(_k)))], [a_(t) = `+`(`*`(diff(... (233)
 

Second, the Coriolis force divided by the mass, which is by definition given by  

> Typesetting:-delayCrossProduct(Physics:-`*`(2, diff(r_(t), t) = `+`(`*`(diff(z(t), t), `*`(_k)), `*`(diff(rho(t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(phi(t), t), `*`(_phi(t)))))), omega_ = `+`(`-`...
 

Typesetting:-mprintslash([`+`(`*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_)))) = `+`(`-`(`*`(2, `*`(rho(t), `*`(`^`(diff(phi(t), t), 2), `*`(_rho(t)))))), `*`(2, `*`(diff(rho(t), t), `*`(d... (234)
 

Third, the centripetal force divided by the mass, which is defined by 

> Typesetting:-delayCrossProduct(omega_ = `+`(`-`(`*`(diff(phi(t), t), `*`(_k)))), Typesetting:-delayCrossProduct(r_(t) = `+`(`*`(z(t), `*`(_k)), `*`(rho(t), `*`(_rho(t)))), omega_ = `+`(`-`(`*`(diff(ph...
 

Typesetting:-mprintslash([Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2), `*`(_rho(t))))], [Physics:-Vectors:-`&x`(omega_, Physics:-Vec... (235)
 

Adding these three terms, 

> Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(a_(t) = `+`(`*`(diff(diff(rho(t), t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(diff(phi(t), t), t), `*`(_phi(t)))), `*`(diff(diff(z(t), t), t), `*`(_k))), `+...
 

Typesetting:-mprintslash([`+`(a_(t), `*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_))), Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_))) = `+`(`*`(`+`(diff(rho(t), `$`(... (236)
 

where the right-hand side of Physics:-Vectors:-`+`(Physics:-Vectors:-`+`(a_(t) = `+`(`*`(diff(diff(rho(t), t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(diff(phi(t), t), t), `*`(_phi(t)))), `*`(diff(diff(z(t), t), t), `*`(_k))), `+... is the same result computed in lhs(diff(diff(r_(t), t), t) = `+`(`*`(diff(diff(x(t), t), t), `*`(_i)), `*`(diff(diff(y(t), t), t), `*`(_j)), `*`(diff(diff(z(t), t), t), `*`(_k)))) = Physics:-Vectors:-ChangeBasis(rhs(diff(diff(r_(t)... 

> diff(diff(r_(t), t), t) = `+`(`*`(`+`(diff(diff(rho(t), t), t), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t), ...
 

Typesetting:-mprintslash([diff(r_(t), `$`(t, 2)) = `+`(`*`(`+`(diff(rho(t), `$`(t, 2)), `-`(`*`(rho(t), `*`(`^`(diff(phi(t), t), 2))))), `*`(_rho(t))), `*`(`+`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi... (237)
 

So that 

> Physics:-Vectors:-diff(r_(t), t, t) = lhs(`+`(a_(t), `*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_))), Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_))) = `+`(`*`(`+`(d...
 

Typesetting:-mprintslash([diff(r_(t), `$`(t, 2)) = `+`(a_(t), `*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_))), Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)))], [dif... (238)
 

 

From this result, in the absence of external forces, Physics:-Vectors:-diff(r_(t), t, t) = 0; and the acceleration a_(t); measured in the rotating system is given by the sum of the Coriolis and centripetal accelerations 

> isolate(rhs(diff(diff(r_(t), t), t) = `+`(a_(t), `*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_))), Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)))), a_(t));
 

Typesetting:-mprintslash([a_(t) = `+`(`-`(`*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_)))), `-`(Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_))))], [a_(t) = `+`(`-`(`... (239)
 

In conclusion: in the absence of external forces, the acceleration of a particle observed in a rotating (non-inertial) system of reference is not zero.  

 

Expressing this equation isolate(rhs(diff(diff(r_(t), t), t) = `+`(a_(t), `*`(2, `*`(Physics:-Vectors:-`&x`(diff(r_(t), t), omega_))), Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)))), a_(t)); in terms of rho(t), phi(t), z(t); we get 

> subs(a_(t) = `+`(`*`(diff(diff(rho(t), t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(diff(phi(t), t), t), `*`(_phi(t)))), `*`(diff(diff(z(t), t), t), `*`(_k))), `+`(`-`(`+`(`*`(2, `*`(Physics:-Vectors:-...
 

Typesetting:-mprintslash([`+`(`*`(diff(rho(t), `$`(t, 2)), `*`(_rho(t))), `*`(rho(t), `*`(diff(phi(t), `$`(t, 2)), `*`(_phi(t)))), `*`(diff(z(t), `$`(t, 2)), `*`(_k))) = `+`(`*`(rho(t), `*`(`^`(diff(p... (240)
 

resulting in the three equations 

> Typesetting:-delayDotProduct(`+`(`*`(diff(diff(rho(t), t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(diff(phi(t), t), t), `*`(_phi(t)))), `*`(diff(diff(z(t), t), t), `*`(_k))) = `+`(`*`(rho(t), `*`(`^`(...
 

Typesetting:-mprintslash([diff(rho(t), `$`(t, 2)) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2)))], [diff(diff(rho(t), t), t) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2)))]) (241)
 

> Typesetting:-delayDotProduct(`+`(`*`(diff(diff(rho(t), t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(diff(phi(t), t), t), `*`(_phi(t)))), `*`(diff(diff(z(t), t), t), `*`(_k))) = `+`(`*`(rho(t), `*`(`^`(...
 

Typesetting:-mprintslash([`*`(rho(t), `*`(diff(phi(t), `$`(t, 2)))) = `+`(`-`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t))))))], [`*`(rho(t), `*`(diff(diff(phi(t), t), t))) = `+`(`-`(`*`(2, `*`(di... (242)
 

> Typesetting:-delayDotProduct(`+`(`*`(diff(diff(rho(t), t), t), `*`(_rho(t))), `*`(rho(t), `*`(diff(diff(phi(t), t), t), `*`(_phi(t)))), `*`(diff(diff(z(t), t), t), `*`(_k))) = `+`(`*`(rho(t), `*`(`^`(...
 

Typesetting:-mprintslash([diff(z(t), `$`(t, 2)) = 0], [diff(diff(z(t), t), t) = 0]) (243)
 

which are the same equations returned by Geodesics in Physics:-Geodesics(t);  

> [diff(diff(rho(t), t), t) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2))), diff(diff(phi(t), t), t) = `+`(`-`(`/`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t))))), diff(diff(z(t), t), t)...
 

Typesetting:-mprintslash([[diff(rho(t), `$`(t, 2)) = `*`(rho(t), `*`(`^`(diff(phi(t), t), 2))), diff(phi(t), `$`(t, 2)) = `+`(`-`(`/`(`*`(2, `*`(diff(rho(t), t), `*`(diff(phi(t), t)))), `*`(rho(t)))))... (244)
 

References 

[1] L.D. Landau, E.M. Lifchitz, Mechanics, Course of Theoretical Physics, Volume 1, third edition, Elsevier.