Work and Force

We investigate the relationship between force and energy (work). For energy-conserving systems one can formulate the equations of motion by making use of scalar quantities that represent the kinetic and potential energy for the particle in the force field. For non-conservative forces the work integral has to be calculated for the individual particle trajectories using a line integral. Both cases are considered in this worksheet.

We begin with the conservative case. Take the example of a two-centre electrostatic potential as experienced in a diatomic molecule by the valence electron (e.g., the H2+ molecule, cf. R.L. Greene, Classical Mechanics with Maple (Springer 1994), pp. 88-93, and M.L. de Jong, Introduction to Computational Physics (Addison-Wesley 1991), chapter 10). The problem of a spaceship, or asteroid, or small planet moving in the field of two nearly stationary stars would be calculated in a completely analogous way. Two protons are kept at a fixed distance 2 d along the z axis, one located at (0,0,- d ), the other at (0,0, d ). The electron is located at ( x , y , z ) and being of opposite charge experiences an attractive force. It is convenient to express the electron distance from either proton in terms of its coordinate ( x , y , z ) using coordinates, r 1 and r 2 respectively.

> restart;

> r1:=sqrt(x^2+y^2+(z-d)^2);

> r2:=sqrt(x^2+y^2+(z+d)^2);

The electrotatic potential energy in MKSA (SI) units Usi is then given in terms of the charges for the protons and the electron, Q and q respectively, the dielectric constant epsilon[0] and the two distances as:

> f[0]:=q*Q/(4*Pi*epsilon[0]);

> Usi:=-f[0]*(1/r1+1/r2);

For plotting as well as calculations it is more convenient to measure energy in units in which the dielectric constant, the charge units and a characteristic length scale are absorbed into a single constant. The characteristic length scale can be chosen as d , (or as 2 d ), i.e., the equilibrium separation of the nuclei in the molecule. In these units we simply have:

> U:=unapply(subs(d=1,-1/r1-1/r2),x,y,z);

There are various ways to visualize this function. It is possible to make cuts such that a 2D graph conveys some information. Two of the coordinates have to be fixed for this purpose, and keep in mind that all distances are measured in terms of d which has been set to unity, i.e., all coordinates are measured as multiples of the half-separation between the nuclei.

We take cuts along the x axis ( x =0) with fixed y (treated as a parameter y 0) and display as a function of z . We restrict the range of the ordinate to avoid a huge scaling from the singularity in the potential at (0,0,- d ) and (0,0, d ) respectively.

> P1:=plot(U(0,0,z),z=-3..3,-5..0,color=blue):

> P2:=plot(U(0,0.5,z),z=-3..3,color=red):

> P3:=plot(U(0,1,z),z=-3..3,color=green):

> with(plots):

> display({P1,P2,P3});

It is straightforward to read off from this diagram the z component of the force for the electron when it is located at (0, y 0, z ) for the three values of y 0 (0, 0.5 and 1) represented by the three curves above. Due to the superposition principle we can treat these potentials as 1D potentials for the motion with fixed x 0 and y 0, and thus the force is simply given as the negative of the derivative w.r.t. z :

> P1f:=plot(-diff(U(0,0,z),z),z=-3..3,-5..5,color=blue,discont=true):

> P2f:=plot(-diff(U(0,0.5,z),z),z=-3..3,color=red):

> P3f:=plot(-diff(U(0,1,z),z),z=-3..3,color=green):

> display({P1f,P2f,P3f});

We added an option in the graph for x 0=0 to avoid a vertical line showing up at the location of the singularity (incorrect joining of the asymptotic behaviour done by Maple without the option). Note that the sign of the force changes, and that the location of the equilibrium points (intercept with zero) depends on the value of y 0. For y 0=0 the behaviour is qualitatively different as a singularity appears.

How can the full 3D behaviour be displayed? The correct approach would be to switch to cylindrical (polar) coordinates and make use of axial symmetry w.r.t. the z axis. We prefer to keep Cartesian coordinates however. Thus we keep the constraint x 0=0, and display the potential in the y - z plane using a surface plot and a contour diagram.

Since we called the plots package already, contourplot is available..

> plot3d(U(0,y,z),y=-3..3,z=-3..3,view=-3..0,style=patchcontour,shading=zhue,axes=boxed,numpoints=3000,orientation=[-30,15]);

The surface plot contains contours drawn at various heights. We can turn it into a contour diagram by simply choosing the right projection (choice of theta/phi angles). Alternatively we use the contourplot command from the plots package. We avoid the singularity by picking a value of x 0=0.1:

> contourplot(U(0.1,y,z),z=-3..3,y=-3..3,axes=boxed,scaling=CONSTRAINED,grid=[40,40],contours=[-4,-3,-2,-1.5,-1],filled=true);

The force acting on the particle is given as the negative of the gradient. We need to load the linalg package and obtain:

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> F:=grad(-U(x,y,z),[x,y,z]);

Note that one could not simply evaluate F:=-grad(U...) as this would amount to the multiplication of a vector by a scalar (namely by -1). Maple requires an evalm wrapped around such a command to be able to evaluate. The evalm procedure is used for matrix and vector multiplications.

To confirm our previous result for the z component of the force we graph:

> plot(subs(x=0,y=0.5,F[3]),z=-3..3,-5..5,color=red);

We can demonstrate that the force field is curl-free:

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

We can demonstrate by example that the line integral is independent of the chosen path. We choose two paths in the y - z plane for which the potential was shown in the contour diagram. Let us calculate the work to reach the location [0,1,1] from the origin [0,0,0] in two different ways. Once through a straight-line connection, and once through an arc of a circle with radius sqrt(2)/2 around the point [0,0.5,0.5].

To set up the line integral to calculate the work done by the force on the particle as it traverses the trajectory r ( t ) we use a parameter representation of the path. If the path is differentiable with respect to its parameter (e.g., time in physics), then the line integral can be reduced to a definite Riemann integral.

W = int ( F . d r ) = int( F ( r ( t )) . d r /d t , t = t0 .. t1)

For the straight-line connection we have:

> r_in:=[0,0,0]; r_fin:=[0,1,1];

> rSL:=evalm(r_in+t*(r_fin-r_in));

> rSL[1], rSL[2], rSL[3];

The above does not work properly without the call to evalm, and the operation below does not work without the map construct, once rSL is a vector.

> vSL:=map(diff,rSL,t);

> whattype(rSL);

In this parameter representation of the straight line the endpoints are realized for t =0 and t =1 respectively. It is important to evaluate the force along the trajectory, i.e., to substitute for x , y , z in the force assignment the values of x ( t ), y ( t ), z ( t ), so that the definite integration makes sense, and the force evaluations happen at the correct positions as a function of t , namely along the chosen trajectory. This is the most important aspect of the line integral to be understood by the beginning physicist (and also the most common source of error).

The op(F) rather than F in the substitution is important. It is needed for the substitution to take place in all components:

> subs(x=rSL[1],y=rSL[2],z=rSL[3],op(F));

We use the dotprod procedure from the linalg package to calculate the dot product between F ( r ( t )) and d r /d t , and then integrate over t :

> W1:=int(dotprod(subs(x=rSL[1],y=rSL[2],z=rSL[3],op(F)),vSL),t=0..1);

We can check whether this result agrees with the difference between the potential energies evaluated at the two endpoints respectively:

> U(op(r_fin))-U(op(r_in));

> evalf(%,4); evalf(W1,4);

To see how the first of the two numbers came about look at the contributions to the potential difference:

> evalf(U(op(r_fin)),4); evalf(U(op(r_in)),4);

The amount of work to move the electron from location [0,0,0] to [0,1,1] is negative, which means that on average the force was acting more against the motion.

Try simpler trajectories for which the force vector is aligned with the position vector, such as moving from [0,0,0] to [0,0,1- c ], i.e., from the unstable equilibrium point towards one of the protons, e.g., the one located at z =1. Also see what happens as your path moves across the singularity...

The potential energy at the initial position r_in=[0,0,0] is more attractive than at r_fin=[0,1,1] . The difference in potential energy between the final and initial points of the trajectory is therefore positive. It costs energy to move from [0,0,0] to [0,1,1]. This is consistent with the previous finding that the work is negative: we have to invest energy to move the electron from [0,0,0] to [0,1,1]. The electron can trade kinetic energy to achieve this goal.

Now we repeat the calculation for a circular trajectory:

> rC:=evalm(sqrt(2)/2*[0,cos(t),sin(t)]+[0,1/2,1/2]);

> rC[1];

> plot([rC[2],rC[3],t=Pi/4..5/4*Pi],scaling=CONSTRAINED);

> map(simplify,subs(t=Pi/4,op(rC)));

> map(simplify,subs(t=5*Pi/4,op(rC)));

With the chosen parametrisation the trajectory starts at 5/4 Pi, and ends at Pi/4. However, it passes through the singularity at [0,0,1], and thus we integrate the other way around:

> map(simplify,subs(t=-3*Pi/4,op(rC)));

> t0:=-3*Pi/4: t1:=Pi/4:

> vC:=map(diff,rC,t);

> W2:=int(expand(dotprod(subs(x=rC[1],y=rC[2],z=rC[3],op(F)),vC)),t=t0..t1);

Maple has a hard time doing the integral in closed form. We see that the integrand is non-trivial in the next command. There is no difficulty in integrating numerically if the singularity is avoided.

> arg:=expand(dotprod(subs(x=rC[1],y=rC[2],z=rC[3],op(F)),vC)):

> evalf(Int(arg,t=t0..t1),4);

We are satisfied that the answer agrees with the floating-point evaluation of the previous work calculation for the straight-line trajectory and the negative of the potential difference.

Now that we established the fact that the potential energy allows to calculate the work done on particles by a curl-free force, we illustrate how one uses the contour diagram to read off information about the force. We fixed a plane when graphing the equipotential contours, and thus we choose to illustrate only the force components in that plane. The force is given by the gradient of the potential (up to a minus sign). The direction of the gradient vector originating at any point is orthogonal to an equipotential contour passing through that point. The plots package in Maple contains a command that allows to draw a set of arrows on a discrete mesh for any given potential function in two dimensions. To avoid the singularity we look at the y - z plane defined by x 0=0.5:

> PL1:=gradplot(-U(0.5,y,z),z=-3..3,y=-2..2,axes=boxed,color=red,grid=[16,24],arrows=thick):

> PL2:=contourplot(U(0.5,y,z),z=-3..3,y=-2..2,axes=boxed,scaling=CONSTRAINED,grid=[40,40],contours=[-2,-1.5,-1,-0.5],filled=true,coloring=[white,blue]):

> display(PL1,PL2);

The arrows give a good impression of the variation of the two force components in the y - z plane. No information is conveyed about the z component of the force.

The diagram illustrates the fact that it is not entirely obvious that the work integral is path-independent: while walking different trajectories from one endpoint to another one encounters force vectors pointing in different directions. Think about a path where the electron is first accelerated towards the proton located in the upper half-plane (at [0,0, d ]), and then swings around towards the final location [0,1,1] working against the force field. A very different experience as compared to the straight-line or circular paths calculated before. Thus, the curl-free nature of conservative forces is a very powerful property.

The potential in the present case is given by a superposition of two central potentials with respect to the positions of the nuclei [0,0, d ] and [0,0,- d ] respectively. Central potentials play an important role due to their spherical symmetry (the potential depends only on the distance r , not on the spherical polar angles theta and phi). This symmetry implies angular momentum conservation, both in classical and in quantum physics.

>