Classical Mechanics with Maple

Section 2.1: Mass Particles in Cartesian, Polar and Natural Coordinates

Harald Kammerer
maple@jademountain.de

Initialisation

> restart;

> libname:="C:/mylib/m6dynlib","C:/mylib/m6dynfig",libname;

> with(plots):with(plottools):with(linalg):with(dynamics);with(figures_chapter_2);

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

2.1.1 Cartesian Coordinates

This course deals with the kinetics of mass particles and rigid bodies. A body with finite mass and negligible dimensions is called a mass particle or simply a particle . A rigid body is a body with fixed distances among all of its points. Mass particles have negligible dimensions. This means that only translation and not rotation is defined for particles. The state of a mass particle under motion is unambiguously defined by its position at time t . The position at time t is defined by the vector r ( t ) which starts at a fixed point O of the inertial system of the particle. Accordingly, the position of the particle at time t+ D t (call: "Delta t") is defined by r ( t+ D t ) (see the animation in Fig. 1). To see the animation, click on the diagram and hit "play" in the menu bar.

> display(Fig_2_1(),insequence=true,scaling=constrained,axes=none,title="Figure 1");

The velocity v ( t ) is defined by

>

Notice that in this expressions r ( t ) and v ( t ) are vectors.

We can write this also as

>

In the last equation, the vector r ( t+ D t ) - r ( t ) points in the direction of the tangent of the curve of the motion when the timestep D t decreases and becomes infinitesimal. So the vector v ( t ) points in this direction, too. To see what happens in this situation click the picture for Figure 1 above and press "play" in the menu bar.

For two different times t and t+ D t we get for the velocities v ( t ) and v ( t+ D t ) (see Fig. 2).

> display(Fig_2_2(),scaling=constrained,axes=none, title="Figure 2");

The acceleration is defined as

>

Notice that in this expresions v ( t ) and a ( t ) are vectors.

We can write this as

>

A cartesian coordinate system is a coordinate system in which the locations of points in space are expressed by reference to three planes, called coordinate planes, no two of which are parallel. We describe this system by three orthogonal axes. The direction in space of every axis is described by the unit vectors e x, e y and e z

> ex:=vector(3,[1,0,0]);

> ey:=vector(3,[0,1,0]);

> ez:=vector(3,[0,0,1]);

The position vector

> r(t):=x(t)*ex+y(t)*ey+z(t)*ez;

with the coordinates x ( t ) , y ( t ) , z ( t ) describes the position of the particle at time t .

The derivative of a variable with respect to time is usually denoted by a dot above the variable. Here it is denoted by an appended letter d. So we write for the coordinates of the velocity

>

>

>

and for the velocity vector

>

Written in the form of a vector we get

> evalm(v(t));

And for the acceleration we write

>

>

>

>

> evalm(a(t));

In general the vector a ( t ) does not have the direction of the tangent of the curve of the motion. We illustrate this in the following example.

Example: Horizontal Throwing of a Ball

Consider a man standing on a small hill. We state a coordinate system with the x-axis in the horizontal direction, the y-axis in the vertical and the z-axis orthogonal to the picture plane. The origin of the system is at the head of the man (see Fig. 3).

> display(Fig_2_3(10,0,0,0,0,15,0),scaling=constrained,axes=none, title="Figure 3");

Now the man throws a ball in the horizontal direction, starting at the origin with velocity v0. Air resistance is not taken into account. The velocity in the horizontal direction then is constant and the acceleration is consequently 0 up to the moment when the ball touches the ground. In the vertical direction, the acceleration is caused by gravity. There is no motion orthogonal to the picture plane .

So we write for the accelerations in the x-, y- and z-direction

>

>

>

where g is the acceleration of gravity = 9.81 .

We get for the acceleration vector

> a(t):=evalm(ax(t)*ex+ay(t)*ey+az(t)*ez);

At all times, the acceleration vector points in the vertical direction. Because the velocity at the beginning points in the horizontal direction we see already now that the directions of the two vectors are not the same.

Integration with respect to time t yields for the velocity vector

>

>

>

From this integration we get the three constants cx , cy and cz . We know the values of the components of the velocity in every direction at the time t=0 . These are the three initial conditions.

> ic_x:=subs(t=0,vx(t))=v0;

> ic_y:=subs(t=0,vy(t))=0;

> ic_z:=subs(t=0,vz(t))=0;

Solving these conditions yields the integration constants

> sol_1:=solve({ic_x,ic_y,ic_z},{cx,cy,cz});

> assign(sol_1);

So we get for the velocity vector

> v(t):=evalm(vx(t)*ex+vy(t)*ey+vz(t)*ez);

One more integration yields the position of the ball at time t

>

>

>

Now we have the three constants cxx , cyy and czz . We know the position of the ball at the beginning of the motion, because we defined this point as the origin. This information yields three initial conditions.

> ic_xx:=subs(t=0,x(t))=0;

> ic_yy:=subs(t=0,y(t))=0;

> ic_zz:=subs(t=0,z(t))=0;

Solving these conditions yields the integration constants

> sol_2:=solve({ic_xx,ic_yy,ic_zz},{cxx,cyy,czz});

> assign(sol_2);

So we get for the position vector

> r(t):=evalm(x(t)*ex+y(t)*ey+z(t)*ez);

With this result the problem is solved. But now we want to make some additional considerations.

Let h be the vertical distance between the starting point of the ball and the ground. Then we can calculate the time period the ball needs on its way up to the moment when it touches the ground from the relation

> eq:=eval(y(t))=h;

The solution of this equation yields the desired time

> Tend:=solve(eq,t);

We get two solutions with the same absolute value but different signs. Only the positive solution makes sense, so we use a little trick (note: here it is not clearly defined whether g and h are positive)

> Tend:=sqrt(Tend[1]**2);

At last we want to calculate the horizontal distance that the ball travels. We get this by substituting the time Tend in the expression for x ( t )

> xend:=subs(t=Tend,eval(x(t)));

To show the result we choose some concrete values for the height h and the velocity v0 .

> h:=10:

> v0:=10:

The acceleration of gravity is

> g:=9.81:

Fig. 4 shows the result. To see what happens click the picture and press "play" in the menu bar. (Notice that the man's arm actually moves =) )

> display(Fig_2_3(h,v0,eval(x(t)),eval(y(t)),Tend,eval(xend),t),insequence=true,scaling=constrained,axes=none, title="Figure 4");

In Fig 5 we see the velocity vector and the acceleration vector together with the trajectory of the motion at different time steps.

> display(Fig_2_4(h,v0,eval(x(t)),eval(y(t)),Tend,eval(xend),t),insequence=false,scaling=constrained,axes=none, title="Figure 5");

To try the example for different initial velocities and different heights, change the values h and v0 above.

2.1.2 Polar Coordinates

Now we consider a coordinate system with rotated initial vectors e r( t ) and e f ( t ) (in the figures, this is written as ephi) against the vectors e x and e y by the vector f( t ) (phi(t) in the figures). The directions of e r and e f are changing. See the presentation in Fig. 6 to get the following relation.

> display(Fig_2_5(),scaling=constrained,axes=none, title="Figure 6");

The relations between the different systems of initial vectors are

>

>

Or written as vectors

> evalm(er(t));

> evalm(ephi(t));

The initial vectors e x and e y are not dependent on time. So we get by the derivative with respect to time for vectors e r and e f

>

>

>

Or written as vectors

> evalm(erd(t));

> evalm(ephid(t));

as the derivative of f ( t ) with respect to time t (spoken as: "e r dot" and ephid: "e f dot").

We can calculate this result by using Maple when we calculate the derivatives for every coordinate of the vectors

> erd(t):=vector(3,[diff(evalm(er(t))[1],t),diff(evalm(er(t))[2],t),diff(evalm(er(t))[3],t)]);

> ephid(t):=vector(3,[diff(evalm(ephi(t))[1],t),diff(evalm(ephi(t))[2],t),diff(evalm(ephi(t))[3],t)]);

Because we consider only plane motions, there are no changes of the position in the direction orthogonal to the plane. So the z -components of the velocity vectors are 0.

For clarity of the exposition that follows, we need to unassign the vectors e r( t ) and e f ( t ) and their derivatives.

> unassign('er(t)','ephi(t)','erd(t)','ephid(t)');

Comparison with the vectors e r( t ) and e f ( t ) yields the relation

>

>

for the velocity.

Now we consider a mass particle on a trajectory by use of polar coordinates. We use the example with the man throwing a ball (see Fig. 4). Fig. 7 shows the situation at an arbitrary time.

> Tend := sqrt(2)*sqrt(1/g*h);

> display(Fig_2_6(10,10,eval(x(t)),eval(y(t)),Tend,eval(xend),t),insequence=false,scaling=constrained,axes=none, title ="Figure 7");

The position of the ball is described by the vector r ( t ).

>

with the distance R(t) between the origin and the current position of the ball and the initial vector e r( t ). This vector is rotated from the fixed x-axes by the angle f ( t ). The initial vector e f ( t ) is choosen orthogonal to e r( t ) in the direction of increasing f. We get for the velocity

>

and with the result above for the derivative of e r( t )

> v(t):=subs(diff(er(t),t)=erd(t),v(t));

For the acceleration we get

> a(t):=diff(v(t),t);

and with the corresponding substitutions

> a(t):=subs({diff(er(t),t)=erd(t),diff(ephi(t),t)=ephid(t)},a(t));

At last we get

> a(t):=collect(a(t),{er(t),ephi(t)});

Now we collect the components of the velocity and the acceleration in the directions of e f ( t ) and e r( t )

> vr(t):=coeff(v(t),er(t));

> vphi(t):=coeff(v(t),ephi(t));

> ar(t):=coeff(a(t),er(t));

> aphi(t):=coeff(a(t),ephi(t));

These relations are included in the dynamics package. See polar_vel , polar_acc to know how to use this functions.

Example: Horizontal Throwing of a Ball (continuation)

As an example we consider again the man who throws a ball from the top of a hill.

But first we unassign the variables to describe the results in a general manner

> unassign('h','v0','g');

Above we had calculated the position vector of the ball

> r(t) := vector([v0*t, 1/2*g*t^2, 0]);

So the distance between the origin and the position of the ball is

>

and for the angle between the x -axis and the position vector we get from Fig. 7

>

Using the relations above we get for the velocity in the form

> v_polar(t):=polar_vel(R(t),phi(t),t);

and for the acceleration in the form

> a_polar(t):=polar_acc(R(t),phi(t),t);

Of course this result is much more complicated than the result in cartesian coordinates. But this example shows how to calculate the velocity and the acceleration in polar coordinates from the given position vector. And with the functions polar_vel and polar_acc the calculation is easy.

With the same concrete values for the height h and the velocity v0 as above

> h:=10:

> v0:=10:

and the acceleration of gravity

> g:=9.81:

we get for the velocity

> vr(t):=eval(v_polar(t)[1]);

> vphi(t):=eval(v_polar(t)[2]);

and for the acceleration

> ar(t):=eval(a_polar(t)[1]);

> aphi(t):=eval(a_polar(t)[2]);

Fig. 8 shows the distance R ( t ) between the origin and the position of the ball, Fig. 9 the angle f ( t ) between the position vector r ( t ) and the x -axis as a function of t

> Plot1p:=plot(R(t),t=0..Tend, color=green, thickness=3, legend="R(t)"):

> Plot2p:=plot(phi(t),t=0..Tend, 0..Pi/2, color=blue, thickness=3, legend="phi(t)"):

> display({Plot1p}, title="Figure 8");

> display({Plot2p}, title="Figure 9");

The result of the calculation is shown in Fig. 10 (velocity) and Fig. 11 (acceleration).

> Plot1v:=plot(vr(t),t=0..Tend, color=green, thickness=3, legend="velocity in r-direction"):

> Plot2v:=plot(vphi(t),t=0..Tend, color=blue, thickness=3, legend="velocity in phi-direction"):

> display({Plot1v,Plot2v}, title="Figure 10");

> Plot1a:=plot(ar(t),t=0..Tend, color=green, thickness=3, legend="acceleration in r-direction"):

> Plot2a:=plot(aphi(t),t=0..Tend, color=blue, thickness=3, legend="acceleration in phi-direction"):

> display({Plot1a,Plot2a}, title="Figure 11");

Remember that the origin of the f ( t )-r( t )-coordinate system is at the starting point of the motion. At time Tend the ball touches the ground, and the motion ends. You can see what happens when the hill on which the man stands is higher or when he throws the ball faster by changing the constants above.

The description of the motion in this example in polar coordinates is not particularly useful, but it helps to understand the relations between cartesian and polar coordinates.

2.1.3 Natural Coordinates

The last method for describing motion in this section is to use natural coordinates . The trajectory of the motion should be given in this case. The position of the partricle is given by the arc length s . Starting point is chosen at s=0 . At the time t the position of the particle is s ( t ). The tangential inital vector e t( t ) points in the tangential direction of increasing s . The normal initial vector e n( t ) points orthogonal to the tangential inital vector in the direction of the center of the circle (see Fig. 12).

> display(Fig_2_7()[1],scaling=constrained,axes=none, title="Figure 12");

Now we consider two states of the motion, first at time t and second at time t + D t (in the following figures D is written "Delta"). The positions can be described by the arc length s ( t ) and s ( t+ D ) or as above by the position vectors r ( t ) and r ( t+ D t ). The direction the particle is moving this time period is D s (Delta_s in the figure). This situation is shown in Fig 13.

> display(Fig_2_7()[2],scaling=constrained,axes=none, title="Figure 13");

The velocity points in the direction of the trajectory of the motion, as was shown above. That means that the velocity at every time t points in the direction of e t( t ) and has no component in the direction of e n( t ). With the derivative of s with respect to time t

>

we get for the velocity

>

Next we want to derive the acceleration of the motion of the particle. The derivative of the velocity with respect to time is

>

The question is at this point what happens after a change of the initial vector e t( t ) in time. In Fig. 14 the velocity vectors v ( t ) and v ( t + D t ) are shown.

> display(Fig_2_7()[3],scaling=constrained,axes=none, title="Figure 14");

The radius of curvature for both considered time points are drawn in Fig. 15.

> display(Fig_2_7()[4],scaling=constrained,axes=none, title="Figure 15");

We relocate the vector of the velocity at time t parallel to the point s ( t+ D ) (Fig. 16)

> display(Fig_2_7()[5],scaling=constrained,axes=none, title="Figure 16");

Next we replace the vectors of the velocity v ( t ) and v ( t+ D t ) by the initial vectors e t( t ) and e t( t+ D t ). Additionally we include the initial vector e n( t ) (see Fig. 17).

> display(Fig_2_7()[6],scaling=constrained,axes=none, title="Figure 17");

The changing of the tangential initial vector while the motion from s ( t ) to s ( t+ D t ) is given by e t( t+ D t )- e t( t ) (see Fig. 18).

> display(Fig_2_7()[7],scaling=constrained,axes=none, title="Figure 18");

Now we consider the situation for D t -> 0. The triangles (in Fig. 19 painted cyan and red) are similar (that means that the angles in both triangles are equal).

> display(Fig_2_7()[8],scaling=constrained,axes=none, title="Figure 19");

So we can see the relation

>

Notice that the angle between e t( t ) and e t( t+ D t ) and between r ( t+ D t ) and r ( t+ D t ) for infinitesimal D t is infinitesimal, too. The arc D s is then approximately straight.

The absolute value of all initial vectors, especially the vector et( t ) is 1.

>

The changing of the tangential initial vector e t( t ) is defined by

>

For D t -> 0 the direction of e t( t+ D t )- e t( t ) is the direction of e n( t ). So we get

>

Together with the relations REL_1 and REL_2 we get

>

And with

>

we get

> RES:=diff(et(t),t) = diff(s(t),t)/rho(t)*en(t);

Using RES we obtain for the acceleration vector

> a(t):=subs(RES,a(t));

The relations for the velocity and the acceleration vector are included in the dynamics package. See natural_vel , natural_acc to see how to use these functions.

Example: Horizontal Throw of a Ball (continuation)

As an example we consider again the man who throws the ball.

First we unassign the variables to describe the results in a general manner

> unassign('h','v0','g');

Above we had calculated for the position vector of the ball

> r(t) := vector([v0*t, 1/2*g*t^2, 0]);

The position of the ball is in natural coordinates described by the arc length. For a given function f( x ) the arc length between the points s1 and s2 is given by

Additionally we need the radius of curvature of the trajectory for every position of the ball and so for every time t.

This is given by

In the example we know the function of the trajectory as a function of time.

With

> x(t):=r(t)[1];

and

> y(t):=r(t)[2];

we get

> y(x):=subs(t=solve(x(t)=x,t),y(t));

Remember that the coordinate y points downwards.

For the arc length of the trajectory we get

>

And the relation for the radius of curvature yields

>

Now we replace the variable x by x ( t )

> s(t):=subs(xi=x(t),s(xi));

> rho(t):=subs(x=x(t),rho(x));

We get for the velocity vector in natural coordinates

> v_natural(t):=natural_vel(s(t),t);

and for the acceleration

> a_natural(t):=natural_acc(s(t),rho(t),t);

We don't look whether this expressions can be simplified but we use the same values as above for the example

> h:=10:

> v0:=10:

> g:=9.81:

With this values we get for the trajectory

> s(t):=eval(s(t));

and the curvature

> rho(t):=eval(rho(t));

Then we get for the velocity

> vt(t):=simplify(eval(v_natural(t)[1]));

> vn(t):=eval(v_natural(t)[2]);

and for the acceleration

> at(t):=simplify(eval(a_natural(t)[1]));

> an(t):=simplify(eval(a_natural(t)[2]));

Fig. 20 shows the trajectory s ( t ) of the motion of the ball, Fig.9 the radius of curvature r ( t ) in dependence of the time t

> Plot1n:=plot(s(t),t=0..Tend, color=green, thickness=3, legend="s(t)"):

> Plot2n:=plot(rho(t),t=0..Tend, color=blue, thickness=3, legend="rho(t)"):

> display({Plot1n}, title="Figure 20");

> display({Plot2n}, title="Figure 21");

This result of the calculation is shown in Fig. 22 (velocity) and Fig. 23 (acceleration).

> Plot1nv:=plot(vt(t),t=0..Tend, color=green, thickness=3, legend="velocity in t-direction"):

> Plot2nv:=plot(vn(t),t=0..Tend, color=blue, thickness=3, legend="velocity in n-direction"):

> display({Plot1nv,Plot2nv}, title="Figure 22");

> Plot1na:=plot(at(t),t=0..Tend, color=green, thickness=3, legend="acceleration in t-direction"):

> Plot2na:=plot(an(t),t=0..Tend, color=blue, thickness=3, legend="acceleration in n-direction"):

> display({Plot1na,Plot2na}, title="Figure 23");

At time Tend the ball touches the ground and the motion ends. You can see what happens when the hill on which the man stands is higher or when he throws the ball faster by changing the constants above.

As with the description of the motion in this example in polar coordinates, it makes no sense in general to describe the solution in natural coordinates, but it helps to understand the relations between cartesian and natural coordinates. Natural coordinates are useful if the trajectory of the motion is prescribed.