Chapter 4: Conservative Forces: Force is minus the gradient of the potential energy

The examples here are supplements to lecture material presented in class.

4.1 Explanation and example

For conservative forces, like gravity or a Hooke's spring, the force is related to derivatives of the potential energy.  For example, for a spring with potential energy 
	U=1/2 k (((x-(Subscript[x, 0 ])))^2), 
the x-component of the force is 
	(Subscript[F, x]) = -((dU )/dx)= -k (x-(Subscript[x, 0])).  
For a 2-dimensional system, the potential energy is 
	U=(1/2) k  ((((x-(Subscript[x, 0])))^2) + (((y-(Subscript[y, 0])))^2)).
This system can be viewed as a spring (of zero length) with one end attached to a peg about which it can freely rotate in a 2-dimensional plane, with a mass attached to the other end.  Given the expression for U as a function of the x and y coordinates of the system, we can obtain the x and y components of the force.
	(Subscript[F, x]) = -((∂U )/(∂x))= -k (x-(Subscript[x, 0])), and
	
	(Subscript[F, y]) = -((∂U )/(∂y))= -k (y-(Subscript[y, 0])).
	
The equations of motion are just
	(Subscript[F, x]) = m (Subscript[a, x])
	-k (x-(Subscript[x, 0])) = m x''	and
	
	(Subscript[F, y]) = m (Subscript[a, y])
	-k (y-(Subscript[y, 0])) = m y''
In this case, we have two equations representing independent harmonic osciallators.  In general, once we have written down the potential energy of the system in terms of the coordinates of the particles, we differentiate to find the forces and then set up the equations we must solve.

More generally, (Overscript[F, ->])= -(Overscript[∇, ->]) U, where  (Overscript[∇, ->])= (Overscript[x, ⋀]) ((∂ )/(∂x)) +  (Overscript[y, ⋀]) ((∂ )/(∂y)) + (Overscript[z, ⋀]) ((∂ )/(∂z)) is the gradient operator in Cartesian coordinates.  In other coordinate systems (polar, spherical coordinates) the gradient operator is more complicated, but the idea is the same.

Remove["Global`*"];

Define the potential energy of the system

u[x_,y_]=((((x-1))^2)+(((y-2))^2))*0.4;

p1=Plot3D[u[x,y],{x,-10,10},{y,-10,10},PlotRange->{{-1,3},{0,4},{0,2}}]

Define the forces

Fx[x_]=-D[u[x,y],x]

Fy[y_]=-D[u[x,y],y]

Having found the components of the force, we can solve Newton's 2nd Law for each component.  Thus, we have two ODE's, one for x'' and one for y'', and both equations depend on x and y, i.e. the equations are coupled.  To fully specify the problem, we also need to provide the components of the initial position and initial velocity of our particle.  Let us assume that the mass of the particle is m=1.

m=1;
sol=NDSolve[{m x''[t]==Fx[x[t]],m y''[t]==Fy[y[t]],x[0]==2,y[0]==3,x'[0]==0,y'[0]==0.2},{x,y},{t,0,200}]

Here, we have solved for 200 units of time.

As we have seen before, NDSolve gives the solution to our coupled ODE's, which should be the two function of time x(t) and y(t).  To get our solutions in a workable form, we do the following.

X[t_]=x[t]/.sol[[1,1]]

Y[t_]=y[t]/.sol[[1,2]]

So now we have the particle's x and y coordinates as functions of time ( X[t] and Y[t] ) under the influence of the potential energy defined by u(x,y).

Let's plot our solution...

p2=ParametricPlot3D[{X[t],Y[t],u[X[t],Y[t]]+0.1},{t,0,200},PlotRange->{{-1,3},{0,4},{0,1}},PlotStyle->{Red}]

Show[p1,p2]

4.2 Example 2

Remove["Global`*"];

We define arbitrarily a potential energy for our system.  We want the potential energy to be a function of the x and y coordinates of the particle position.

h1[x_,y_]=Exp[-((x^2)+(y^2))/3]

Let's see what this potential energy function looks like.  Let's interpret it as the gravitational potential energy of a ball on a hilly landscape.

Plot3D[h1[x,y],{x,-10,10},{y,-10,10},PlotRange->{{-10,10},{-10,10},{0,1}}]

The shape of the function tells us that our particle will feel a repulsive force near the origin, and that this repulsion weakens as the particle moves away from the origin.  There is nothing to stop the particle from going out to infinity.  You can think of this function as the shape of a hill, and the potential energy is given by the height.

Let's add another piece to our potential energy.

h2[x_,y_]=((x^2)+(y^2))/100

It looks like ...

Plot3D[h2[x,y],{x,-10,10},{y,-10,10},PlotRange->{{-10,10},{-10,10},{0,1}}]

The shape of this function tells us that as the particle moves away from the origin, it feels a repulsion pushing it back towards the origin. h2 is just a Hooke's spring kind of potenial.

We make the potential energy of the system the sum of these contributions.

h[x_,y_]=h1[x,y]+h2[x,y]

Our potential energy looks like... (note the change of plotting scale).

p1=Plot3D[h[x,y],{x,-10,10},{y,-10,10},PlotRange->{{-5,5},{-5,5},{0,2}}]

Now we can find the x and y components of the force by taking derivatives of the potential energy with respect to the particle coordinates x and y.

Fx[x_,y_]=-D[h[x,y],x]

Note that in this case the x-component of the force is a function of both x and y.

Fy[x_,y_]=-D[h[x,y],y]

Having found the components of the force, we can solve Newton's 2nd Law for each component.  Thus, we have two ODE's, one for x and one for y, and both equations depend on x and y, i.e. the equations are coupled.  To fully specify the problem, we also need to provide the components of the initial position and initial velocity of our particle.  Let us assume that the mass of the particle is 1.

m=1;
sol=NDSolve[{m x''[t]==Fx[x[t],y[t]],m y''[t]==Fy[x[t],y[t]],x[0]==3,y[0]==0,x'[0]==0,y'[0]==0.2},{x,y},{t,0,200}]

Here, we have solved for 200 units of time.

As we have seen before, NDSolve gives the solution to our coupled ODE's, which should be the two function of time x(t) and y(t).  To get our solutions in a workable form, we do the following.

X[t_]=x[t]/.sol[[1,1]]

Y[t_]=y[t]/.sol[[1,2]]

So now we have the particle's x and y coordinates as functions of time ( X[t] and Y[t] ) under the influence of the potential energy defined by h(x,y).

Let's plot our solution...

p2=ParametricPlot3D[{X[t],Y[t],h[X[t],Y[t]]+0.1},{t,0,200},PlotRange->{{-5,5},{-5,5},{0,1}}];

Show[p1,p2]

Manipulate[Show[p1,ParametricPlot3D[{X[t],Y[t],h[X[t],Y[t]]+0.1},{t,0,s},PlotStyle->{Black,Thick}],PlotRange->{{-5,5},{-5,5},{0,2}}],{s,1,199, 1}]

If you like, go back and give different initial conditions to the particle, or change the potential energy function and see how things change.

4.3 Gravity with 2 particles

Remove["Global`*"]

u[r_]:=-((G m1 m2)/r)

r:=Sqrt[(x2[t]-x1[t])^2+(y2[t]-y1[t])^2]

F1x=-D[u[r],x1[t]]
F1y=-D[u[r],y1[t]]
F2x=-D[u[r],x2[t]]
F2y=-D[u[r],y2[t]]

eq1x=m1 x1''[t]==F1x
eq1y=m1 y1''[t]==F1y
eq2x=m2 x2''[t]==F2x
eq2y=m2 y2''[t]==F2y

G=1;
m1=100;
m2=1;
T=20;

sol=NDSolve[{eq1x,eq1y,eq2x,eq2y,x1[0]==0,y1[0]==0,x2[0]==5,y2[0]==0,x1'[0]==0,y1'[0]==0,x2'[0]==-1,y2'[0]==2},{x1[t],y1[t],x2[t],y2[t]},{t,0,T}]

{X1[t_],Y1[t_],X2[t_],Y2[t_]}={x1[t],y1[t],x2[t],y2[t]}/.sol[[1]]

ParametricPlot[{{X1[t],Y1[t]},{X2[t],Y2[t]}},{t,0,T}]

4.4 Two Lennard-Jones particles starting from rest

We have seen the Lennard-Jones potential in class, which gives the potential energy between two particles as a function of their separation.  Here we will simulate two Lennard-Jones particles starting from rest.  As a result, all their motion will be confined on a line, which will take to be our x-axis.

Remove["Global`*"]

Taking r to be the distance between the two particles, the Lennard-Jones potential is 
V[r_]=4 ϵ(((((σ/r)))^12)-((((σ/r)))^(-6)))
Here, we set ϵ=1 and σ=1.

V[r_]=4 ((r^(-12))-(r^(-6)))

Plot[V[r],{r,0,4},PlotRange->{{0,5},{-1.5,3}}]

The potential energy of the system is simply that of the Lennard-Jones interaction between the two particles
	U = (Subscript[V, 12]).	
Here, (Subscript[V, 12]) is V[(Subscript[r, 12])], the Lennard-Jones potential evaluated at r corresponding to the distance between particles 1 and 2.

Since we are confing our particles to the x-axis, and the repulsion at short range prevents the particles from passing, the distance between the particles is simply

r[x1_,x2_]:=x2-x1

where x1 and x2 are the positions of particles 1 and 2, with particle 2 always being to the right.

As we wrote before, the total potential energy is
	U = (Subscript[V, 12]) = V[r[x2,x1]]

Energy=V[r[x1[t],x2[t]]]

Now that we have the potential energy, we can derive our forces.  For example, the x-component of the force on particle 1 is

-D[Energy,x1[t]]

As you know, we need to give initial conditions before NDSolve can solve our equations of motion.

x1in=0
x2in=1.2
x1vin=0
x2vin=0

Let's also assume the mass of each particle is 1.

m=1;
seteqn={m x1''[t]==-D[Energy,x1[t]],m x2''[t]==-D[Energy,x2[t]],x1'[0]==x1vin,x2'[0]==x2vin,x1[0]==x1in,x2[0]==x2in}

Another list that NDSolve needs is a list of variables...

setvar={x1,x2}

Finally, after much work, we can invoke NDSolve!  Having all the lists set up already really simplifies the look of NDSolve.

tf=5;
sol=Flatten[NDSolve[seteqn,setvar,{t,0,tf}]]

The following command does the substitutions required to extract our functions.

X1[t_]=x1[t]/.sol
X2[t_]=x2[t]/.sol

And now, to plot our results.

Plot[{X1[t],X2[t]},{t,0,tf}]

ListPlot[{{X1[tf],0},{X2[tf],0}},PlotRange->{{-1,2},{-2,2}},PlotStyle->{Red,PointSize[0.1]}]

Animate[ListPlot[{{X1[t],0},{X2[t],0}},PlotRange->{{-1,2},{-2,2}},PlotStyle->{Red,PointSize[0.1]}],{t,0,tf,0.1}]

If all has gone well, we should have a little movie of the two particles vibrating.

Let's plot the potential energy as a function of time.  What do the minima and maxima mean?

En[t_]=V[r[X1[t],X2[t]]];

Plot[En[t],{t,0,tf},PlotRange->{All,All}]

4.5 Mini molecular dynamics simulation

Now let's go through a more complicated example.  We have seen the Lennard-Jones potential in class, which gives the potential energy between two particles as a function of their separation.

Remove["Global`*"]

V[r_]=4 ((r^(-12))-(r^(-6)))

Plot[V[r],{r,0,5},PlotRange->{{0,5},{-1.5,4}}]

Let's have three particles.

n=3

The potential energy of the system is the sum of the pairwise interactions.  There is a potential energy contribution from the interaction between particles 1 and 2, another contribution from the interaction between particles 1 and 3, and finally another contribution from particles 2 and 3.

	U = (Subscript[V, 12]) + (Subscript[V, 13]) + (Subscript[V, 23]).
	
Here, (Subscript[V, 12]) is V[(Subscript[r, 12])], the Lennard-Jones potential evaluated at r corresponding to the distance between particles 1 and 2.

To keep things a little organized, let's define a simpe function that returns the distance given diplacements in x and y.

r[x_,y_]=Sqrt[(x^2)+(y^2)]

So (Subscript[r, 12 ]) would be just r[(Subscript[x, 2])-(Subscript[x, 1]),(Subscript[y, 2])-(Subscript[y, 1])], Where (Subscript[x, 1])is the x-coordinate of particle 1, (Subscript[y, 1]) is the y-coordinate of particle 1 and so on.

r[(Subscript[x, 2])-(Subscript[x, 1]),(Subscript[y, 2])-(Subscript[y, 1])]

As we wrote before, the total potential energy is
	U = (Subscript[V, 12]) + (Subscript[V, 13]) + (Subscript[V, 23])
or	U = V[(Subscript[r, 12])] + V[(Subscript[r, 13])] + V[(Subscript[r, 23])]
or	U = V[r[(Subscript[x, 2])-(Subscript[x, 1]),(Subscript[y, 2])-(Subscript[y, 1])]] + V[r[(Subscript[x, 3])-(Subscript[x, 1]),(Subscript[y, 3])-(Subscript[y, 1])]] + V[r[(Subscript[x, 3])-(Subscript[x, 2]),(Subscript[y, 3])-(Subscript[y, 2])]]

We can enter this by hand.

Energy=V[r[(Subscript[x, 2])[t]-(Subscript[x, 1])[t],(Subscript[y, 2])[t]-(Subscript[y, 1])[t]]]+V[r[(Subscript[x, 3])[t]-(Subscript[x, 1])[t],(Subscript[y, 3])[t]-(Subscript[y, 1])[t]]]+V[r[(Subscript[x, 3])[t]-(Subscript[x, 2])[t],(Subscript[y, 3])[t]-(Subscript[y, 2])[t]]]

Working with 3 particles is still managable.  But what if we had 1000 particles or more?  Writing the sum out by hand gets tedious, and it's easy to introduce errors.  Another way to write down the potential energy is to use summation notation.

Energy=(∑_(i=1)^(n-1))(∑_(j=i+1)^n)V[r[(Subscript[x, j])[t]-(Subscript[x, i])[t],(Subscript[y, j])[t]-(Subscript[y, i])[t]]]

Using the double summation, we have summed over every pair of interaction.  We have been careful not to include distances between a particle and itself, as well we have avoided counting an interaction twice.  It's handy that Mathematica can treat subscipts as counters.

Now that we have the potential energy, we can derive our forces.  For example, the x-component of the force on particle 1 is

(Subscript[Fx, 1])=-D[Energy,(Subscript[x, 1])[t]]

Using a general expression for the x- and y-component of the force on pariticle i, we get

(Subscript[Fx, i_]):=-D[Energy,(Subscript[x, i])[t]]
(Subscript[Fy, i_]):=-D[Energy,(Subscript[y, i])[t]]

Note how we use the underscore on the subscript in a similar fashion to the way we have defined functions.

As you know, we need to give initial conditions before NDSolve can solve our equations of motion.

xin={0,0,3}
yin={-1,1.5,0}
xvin={0,0,0}
yvin={0,0,0}

With this notation
e.g.
	(Subscript[x, 1])[0]=xin[[1]]
and
	(Subscript[x, 1])'[0]=xvin[[1]]
	
i.e. xin is a list (or vector) that contains the initial x-coordinates of the particles.  xin[[1]] is the initial x-position of particle 1, while yin[[2]] is the initial y-position of particle 2 and so on... similarly for the initial velocities.

We could manually type in the six differential equations (3 particles, each with 2 components of position) and twelve initial conditions.  But since all the equations are similar, we can automate the process.  In the end, we want a list of equations, and this list should look like what we normally put into NDSolve.

The following command puts in a list Newton's 2nd Law for x and y along with inital condition equations for x, x', y, y' for each of our particles.  We assume that each particle has mass 1.

m=1;
seteqn=Flatten[Table[{m (Subscript[x, i])''[t]==-D[Energy,(Subscript[x, i])[t]],m (Subscript[y, i])''[t]==-D[Energy,(Subscript[y, i])[t]],(Subscript[x, i])'[0]==xvin[[i]],(Subscript[y, i])'[0]==yvin[[i]],(Subscript[x, i])[0]==xin[[i]],(Subscript[y, i])[0]==yin[[i]]},{i,1,n}]]

In this case, the Flatten command is used to put all the entries into a single list (just like NDSolve needs).  Without Flatten, the set of equations for each particle is put in a minilist.  Try removing Flatten if you like.

By using the Table command, we have provided a way of easily writing down the set of equations for n particles.

Another list that NDSolve needs is a list of variables...

setvar=Flatten[Table[{(Subscript[x, i]),(Subscript[y, i])},{i,1,n}]]

Finally, after much work, we can invoke NDSolve!  Having all the lists set up already really simplifies the look of NDSolve.

sol=Flatten[NDSolve[seteqn,setvar,{t,0,50}]]

The following command does the substitutions required to extract our functions.

Do[(Subscript[X, i])[t_]=(Subscript[x, i])[t]/.sol; (Subscript[Y, i])[t_]=(Subscript[y, i])[t]/.sol,{i,1,n}]

And now, to plot our results.

Table[(Subscript[gr, i])=ParametricPlot[{(Subscript[X, i])[t],(Subscript[Y, i])[t]},{t,0,50}],{i,1,n}]

Show[{(Subscript[gr, 1]),(Subscript[gr, 2]),(Subscript[gr, 3])},PlotRange->{All,All}]

ListPlot[{{(Subscript[X, 2])[50],(Subscript[Y, 2])[50]}},PlotRange->{{-4,4},{-4,4}},PlotStyle->{Red,PointSize[0.1]}]

Manipulate[ListPlot[{{(Subscript[X, 1])[t],(Subscript[Y, 1])[t]},{(Subscript[X, 2])[t],(Subscript[Y, 2])[t]},{(Subscript[X, 3])[t],(Subscript[Y, 3])[t]}},PlotRange->{{-5,5},{-5,5}},PlotStyle->{Red,PointSize[0.1]}],{t,0,50,0.25}]

If all has gone well, we should have a little movie of the three particle moving around.  In this particular case, a near three-body collision results in the 3rd particle being kicked out of the group.  Is this how liquids evaporate?

Let's plot the potential energy as a function of time.

En[t_]=(∑_(i=1)^(n-1))(∑_(j=i+1)^n)V[r[(Subscript[X, j])[t]-(Subscript[X, i])[t],(Subscript[Y, j])[t]-(Subscript[Y, i])[t]]];

Plot[En[t],{t,0,50},PlotRange->{All,All}]

You can see where the 3rd particle collides with the other two.  You can also see the other two particles oscillating.

The example is an illustration of Molecular Dynamics simulation.  By simulating many particles, researchers use this tool to calculate (measure?) the physical and chemical properties of materials.

Here is another set of initial conditions... one in which the particles are nearly collinear.  This is an unstable configuration, and so the particles "melt".

xin={-1.12,0,1.12}
yin={-0.1,-0.2,0.1}
xvin={0,0,0}
yvin={0,0,0}