Classical mechanics is a topic which has been taught intensively over several centuries. It is, with its many variants and ways of presenting the educational material, normally the first real physics course many of us meet and it lays the foundation for further physics studies. Many of the equations and ways of reasoning about the underlying laws of motion and pertinent forces, shape our approaches and understanding of the scientific method and discourse, as well as the way we develop our insights and deeper understanding about physical systems.
There is a wealth of well-tested (from both a physics point of view and a pedagogical standpoint) exercises and problems which can be solved analytically. However, many of these problems represent idealized and less realistic situations. The large majority of these problems are solved by paper and pencil and are traditionally aimed at what we normally refer to as continuous models from which we may find an analytical solution. As a consequence, when teaching mechanics, it implies that we can seldomly venture beyond an idealized case in order to develop our understandings and insights about the underlying forces and laws of motion.
On the other hand, numerical algorithms call for approximate discrete models and much of the development of methods for continuous models are nowadays being replaced by methods for discrete models in science and industry, simply because much larger classes of problems can be addressed with discrete models, often by simpler and more generic methodologies.
As we will see below, when properly scaling the equations at hand, discrete models open up for more advanced abstractions and the possibility to study real life systems, with the added bonus that we can explore and deepen our basic understanding of various physical systems
Analytical solutions are as important as before. In addition, such solutions provide us with invaluable benchmarks and tests for our discrete models. Such benchmarks, as we will see below, allow us to discuss possible sources of errors and their behaviors. And finally, since most of our models are based on various algorithms from numerical mathematics, we have a unique oppotunity to gain a deeper understanding of the mathematical approaches we are using.
With computing and data science as important elements in essentially all aspects of a modern society, we could then try to define Computing as solving scientific problems using all possible tools, including symbolic computing, computers and numerical algorithms, and analytical paper and pencil solutions. Computing provides us with the tools to develope our own understanding of the scientific method by enhancing algorithmic thinking.
The way we will teach this course reflects this definition of computing. The course contains both classical paper and pencil exercises as well as computational projects and exercises. The hope is that this will allow you to explore the physics of systems governed by the degrees of freedom of classical mechanics at a deeper level, and that these insights about the scientific method will help you to develop a better understanding of how the underlying forces and equations of motion and how they impact a given system. Furthermore, by introducing various numerical methods via computational projects and exercises, we aim at developing your competences and skills about these topics.
These competences will enable you to
The power of the scientific method lies in identifying a given problem as a special case of an abstract class of problems, identifying general solution methods for this class of problems, and applying a general method to the specific problem (applying means, in the case of computing, calculations by pen and paper, symbolic computing, or numerical computing by ready-made and/or self-written software). This generic view on problems and methods is particularly important for understanding how to apply available, generic software to solve a particular problem.
However, verification of algorithms and understanding their limitations requires much of the classical knowledge about continuous models.
Before we venture into a reminder on Python and mechanics relevant applications, let us briefly outline some of the abovementioned topics using an example many of you may have seen before in for example CMSE201. A simple algorithm for integration is the Trapezoidal rule. Integration of a function \( f(x) \) by the Trapezoidal Rule is given by following algorithm for an interval \( x \in [a,b] \) $$ \int_a^b(f(x) dx = \frac{1}{2}\left [f(a)+2f(a+h)+\dots+2f(b-h)+f(b)\right] +O(h^2), $$ where \( h \) is the so-called stepsize defined by the number of integration points \( N \) as \( h=(b-a)/(n) \). Python offers an extremely versatile programming environment, allowing for the inclusion of analytical studies in a numerical program. Here we show an example code with the trapezoidal rule. We use also SymPy to evaluate the exact value of the integral and compute the absolute error with respect to the numerically evaluated one of the integral \( \int_0^1 dx x^2 = 1/3 \). The following code for the trapezoidal rule allows you to plot the relative error by comparing with the exact result. By increasing to \( 10^8 \) points one arrives at a region where numerical errors start to accumulate.
from math import log10
import numpy as np
from sympy import Symbol, integrate
import matplotlib.pyplot as plt
# function for the trapezoidal rule
def Trapez(a,b,f,n):
h = (b-a)/float(n)
s = 0
x = a
for i in range(1,n,1):
x = x+h
s = s+ f(x)
s = 0.5*(f(a)+f(b)) +s
return h*s
# function to compute pi
def function(x):
return x*x
# define integration limits
a = 0.0; b = 1.0;
# find result from sympy
# define x as a symbol to be used by sympy
x = Symbol('x')
exact = integrate(function(x), (x, a, b))
# set up the arrays for plotting the relative error
n = np.zeros(9); y = np.zeros(9);
# find the relative error as function of integration points
for i in range(1, 8, 1):
npts = 10**i
result = Trapez(a,b,function,npts)
RelativeError = abs((exact-result)/exact)
n[i] = log10(npts); y[i] = log10(RelativeError);
plt.plot(n,y, 'ro')
plt.xlabel('n')
plt.ylabel('Relative error')
plt.show()
This example shows the potential of combining numerical algorithms with symbolic calculations, allowing us to
Our studies will start with the motion of different types of objects such as a falling ball, a runner, a bicycle etc etc. It means that an object's position in space varies with time. In order to study such systems we need to define
Throughout this course we will use the standardized SI units. The standard unit for length is thus one meter 1m, for mass one kilogram 1kg, for time one second 1s, for force one Newton 1kgm/s$^2$ and for energy 1 Joule 1kgm$^2$s$^{-2}$.
We will use the following notations for various variables (vectors are always boldfaced in these lecture notes):
It is also important to keep track of dimensionalities. Don't mix this up with a chosen unit for a given variable. We mark the dimensionality in these lectures as \( [a] \), where \( a \) is the quantity we are interested in. Thus
Note: This section is under revision
In these lectures we will use boldfaced lower-case letters to label a vector. A vector \( \boldsymbol{a} \) in three dimensions is thus defined as $$ \boldsymbol{a} =(a_x,a_y, a_z), $$ and using the unit vectors in a cartesian system we have $$ \boldsymbol{a} = a_x\boldsymbol{e}_x+a_y\boldsymbol{e}_y+a_z\boldsymbol{e}_z, $$ where the unit vectors have magnitude \( \vert\boldsymbol{e}_i\vert = 1 \) with \( i=x,y,z \).
Using the fact that multiplication of reals is distributive we can show that $$ \boldsymbol{a}(\boldsymbol{b}+\boldsymbol{c})=\boldsymbol{a}\boldsymbol{b}+\boldsymbol{a}\boldsymbol{c}, $$ Similarly we can also show that (using product rule for differentiating reals) $$ \frac{d}{dt}(\boldsymbol{a}\boldsymbol{b})=\boldsymbol{a}\frac{d\boldsymbol{b}}{dt}+\boldsymbol{b}\frac{d\boldsymbol{a}}{dt}. $$
We can repeat these operations for the cross products and show that they are distribuitive $$ \boldsymbol{a}\times(\boldsymbol{b}+\boldsymbol{c})=\boldsymbol{a}\times\boldsymbol{b}+\boldsymbol{a}\times\boldsymbol{c}. $$ We have also that $$ \frac{d}{dt}(\boldsymbol{a}\times\boldsymbol{b})=\boldsymbol{a}\times\frac{d\boldsymbol{b}}{dt}+\boldsymbol{b}\times\frac{d\boldsymbol{a}}{dt}. $$ The rotation of a three-dimensional vector \( \boldsymbol{a}=(a_x,a_y,a_z) \) in the \( xy \) plane around an angle \( \phi \) results in a new vector \( \boldsymbol{b}=(b_x,b_y,b_z) \). This operation can be expressed in terms of linear algebra as a matrix (the rotation matrix) multiplied with a vector. We can write this as $$ \begin{bmatrix} b_x \\ b_y \\ b_z \end{bmatrix} = \begin{bmatrix} \cos{\phi} & \sin{\phi} & 0 \\ -\sin{\phi} & \cos{\phi} & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix} a_x \\ a_y \\ a_z \end{bmatrix}. $$ We can write this in a more compact form as \( \boldsymbol{b} = \boldsymbol{R}\boldsymbol{a} \), where the rotation matrix is defined as $$ \boldsymbol{R} = \begin{bmatrix} \cos{\phi} & \sin{\phi} & 0 \\ -\sin{\phi} & \cos{\phi} & 0 \\ 0 & 0 & 1\end{bmatrix}. $$
We anticipate the mathematical model to come and assume that we have a model for the motion of a falling baseball without air resistance. Our system (the baseball) is at an initial height \( y_0 \) (which we will specify in the program below) at the initial time \( t_0=0 \). In our program example here we will plot the position in steps of \( \Delta t \) up to a final time \( t_f \). The mathematical formula for the position \( y(t) \) as function of time \( t \) is $$ y(t) = y_0-\frac{1}{2}gt^2, $$ where $g=9.80665=0.980655\times 10^1$m/s${}^2$ is a constant representing the standard acceleration due to gravity. We have here adopted the conventional standard value. This does not take into account other effects, such as buoyancy or drag. Furthermore, we stop when the ball hits the ground, which takes place at $$ y(t) = 0= y_0-\frac{1}{2}gt^2, $$ which gives us a final time \( t_f=\sqrt{2y_0/g} \).
As of now we simply assume that we know the formula for the falling object. Afterwards, we will derive it.
We start with preparing folders for storing our calculations, figures and if needed, specific data files we use as input or output files.
# Common imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
#in case we have an input file we wish to read in
#infile = open(data_path("MassEval2016.dat"),'r')
You could also define a function for making our plots. You can obviously avoid this and simply set up various matplotlib commands every time you need them. You may however find it convenient to collect all such commands in one function and simply call this function.
from pylab import plt, mpl
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
def MakePlot(x,y, styles, labels, axlabels):
plt.figure(figsize=(10,6))
for i in range(len(x)):
plt.plot(x[i], y[i], styles[i], label = labels[i])
plt.xlabel(axlabels[0])
plt.ylabel(axlabels[1])
plt.legend(loc=0)
Thereafter we start setting up the code for the falling object.
%matplotlib inline
import matplotlib.patches as mpatches
g = 9.80655 #m/s^2
y_0 = 10.0 # initial position in meters
DeltaT = 0.1 # time step
# final time when y = 0, t = sqrt(2*10/g)
tfinal = np.sqrt(2.0*y_0/g)
#set up arrays
t = np.arange(0,tfinal,DeltaT)
y =y_0 -g*.5*t**2
# Then make a nice printout in table form using Pandas
import pandas as pd
from IPython.display import display
data = {'t[s]': t,
'y[m]': y
}
RawData = pd.DataFrame(data)
display(RawData)
plt.style.use('ggplot')
plt.figure(figsize=(8,8))
plt.scatter(t, y, color = 'b')
blue_patch = mpatches.Patch(color = 'b', label = 'Height y as function of time t')
plt.legend(handles=[blue_patch])
plt.xlabel("t[s]")
plt.ylabel("y[m]")
save_fig("FallingBaseball")
plt.show()
Here we used pandas (see below) to systemize the output of the position as function of time.
Since we have discretized the variables, we introduce the counter \( i \) and let \( y(t)\rightarrow y(t_i)=y_i \) and \( t\rightarrow t_i \) with \( i=0,1,\dots, n \). This gives us the following shorthand notations that we will use for the rest of this course. We define $$ y_i = y(t_i),\hspace{0.2cm} i=0,1,2,\dots,n. $$ This applies to other variables which depend on say time. Examples are the velocities, accelerations, momenta etc. Furthermore we use the shorthand $$ y_{i\pm 1} = y(t_i\pm \Delta t),\hspace{0.12cm} i=0,1,2,\dots,n. $$
A note on notations: When writing for example the velocity as \( v(t) \) we are then referring to the continuous and instantaneous value. A subscript like \( v_i \) refers always to the discretized values.
We can rewrite the instantaneous acceleration as $$ a(t) = \frac{dv}{dt}=\frac{d}{dt}\frac{dy}{dt}=\frac{d^2y}{dt^2}. $$ This forms the starting point for our definition of forces later. It is a famous second-order differential equation. If the acceleration is constant we can now recover the formula for the falling ball we started with. The acceleration can depend on the position and the velocity. To be more formal we should then write the above differential equation as $$ \frac{d^2y}{dt^2}=a(t,y(t),\frac{dy}{dt}). $$
With given initial conditions for \( y(t_0) \) and \( v(t_0) \) we can then integrate the above equation and find the velocities and positions at a given time \( t \).
If we multiply with mass, we have one of the famous expressions for Newton's second law, $$ F(y,v,t)=m\frac{d^2y}{dt^2}=ma(t,y(t),\frac{dy}{dt}), $$ where \( F \) is the force acting on an object with mass \( m \). We see that it also has the right dimension, mass times length divided by time squared. We will come back to this soon.
Formally we can then, starting with the acceleration (suppose we have measured it, how could we do that?) compute say the height of a building. To see this we perform the following integrations from an initial time \( t_0 \) to a given time \( t \) $$ \int_{t_0}^t dt a(t) = \int_{t_0}^t dt \frac{dv}{dt} = v(t)-v(t_0), $$ or as $$ v(t)=v(t_0)+\int_{t_0}^t dt a(t). $$
When we know the velocity as function of time, we can find the position as function of time starting from the defintion of velocity as the derivative with respect to time, that is we have $$ \int_{t_0}^t dt v(t) = \int_{t_0}^t dt \frac{dy}{dt} = y(t)-y(t_0), $$ or as $$ y(t)=y(t_0)+\int_{t_0}^t dt v(t). $$
These equations define what is called the integration method for finding the position and the velocity as functions of time. There is no loss of generality if we extend these equations to more than one spatial dimension.
Let us compute the velocity using the constant value for the acceleration given by \( -g \). We have $$ v(t)=v(t_0)+\int_{t_0}^t dt a(t)=v(t_0)+\int_{t_0}^t dt (-g). $$ Using our initial time as $t_0=0$s and setting the initial velocity $v(t_0)=v_0=0$m/s we get when integrating $$ v(t)=-gt. $$ The more general case is $$ v(t)=v_0-g(t-t_0). $$ We can then integrate the velocity and obtain the final formula for the position as function of time through $$ y(t)=y(t_0)+\int_{t_0}^t dt v(t)=y_0+\int_{t_0}^t dt v(t)=y_0+\int_{t_0}^t dt (-gt), $$ With $y_0=10$m and $t_0=0$s, we obtain the equation we started with $$ y(t)=10-\frac{1}{2}gt^2. $$
After this mathematical background we are now ready to compute the mean velocity using our data.
# Now we can compute the mean velocity using our data
# We define first an array Vaverage
n = np.size(t)
Vaverage = np.zeros(n)
for i in range(1,n-1):
Vaverage[i] = (y[i+1]-y[i])/DeltaT
# Now we can compute the mean accelearatio using our data
# We define first an array Aaverage
n = np.size(t)
Aaverage = np.zeros(n)
Aaverage[0] = -g
for i in range(1,n-1):
Aaverage[i] = (Vaverage[i+1]-Vaverage[i])/DeltaT
data = {'t[s]': t,
'y[m]': y,
'v[m/s]': Vaverage,
'a[m/s^2]': Aaverage
}
NewData = pd.DataFrame(data)
display(NewData[0:n-2])
Note that we don't print the last values!
In our discussions till now of the falling baseball, we have ignored air resistance and simply assumed that our system is only influenced by the gravitational force. We will postpone the derivation of air resistance till later, after our discussion of Newton's laws and forces.
For our discussions here it suffices to state that the accelerations is now modified to $$ \boldsymbol{a}(t) = -g +D\boldsymbol{v}(t)\vert v(t)\vert, $$ where \( \vert v(t)\vert \) is the absolute value of the velocity and \( D \) is a constant which pertains to the specific object we are studying. Since we are dealing with motion in one dimension, we can simplify the above to $$ a(t) = -g +Dv^2(t). $$ We can rewrite this as a differential equation $$ a(t) = \frac{dv}{dt}=\frac{d^2y}{dt^2}= -g +Dv^2(t). $$
Using the integral equations discussed above we can integrate twice and obtain first the velocity as function of time and thereafter the position as function of time.
For this particular case, we can actually obtain an analytical solution for the velocity and for the position. Here we will first compute the solutions analytically, thereafter we will derive Euler's method for solving these differential equations numerically.
For simplicity let us just write \( v(t) \) as \( v \). We have $$ \frac{dv}{dt}= -g +Dv^2(t). $$
We can solve this using the technique of separation of variables. We isolate on the left all terms that involve \( v \) and on the right all terms that involve time. We get then $$ \frac{dv}{g -Dv^2(t) }= -dt, $$
We scale now the equation to the left by introducing a constant \( v_T=\sqrt{g/D} \). This constant has dimension length/time. Can you show this?
Next we integrate the left-hand side (lhs) from \( v_0=0 \) m/s to \( v \) and the right-hand side (rhs) from \( t_0=0 \) to \( t \) and obtain $$ \int_{0}^v\frac{dv}{g -Dv^2(t) }= \frac{v_T}{g}\mathrm{arctanh}(\frac{v}{v_T}) =-\int_0^tdt = -t. $$ We can reorganize these equations as $$ v_T\mathrm{arctanh}(\frac{v}{v_T}) =-gt, $$ which gives us \( v \) as function of time $$ v(t)=v_T\tanh{-(\frac{gt}{v_T})}. $$
With the velocity we can then find the height \( y(t) \) by integrating yet another time, that is $$ y(t)=y(t_0)+\int_{t_0}^t dt v(t)=\int_{0}^t dt[v_T\tanh{-(\frac{gt}{v_T})}]. $$ This integral is a little bit trickier but we can look it up in a table over known integrals and we get $$ y(t)=y(t_0)-\frac{v_T^2}{g}\log{[\cosh{(\frac{gt}{v_T})}]}. $$ Alternatively we could have used the symbolic Python package Sympy (example will be inserted later).
In most cases however, we need to revert to numerical solutions.
Here we will try the simplest possible approach to solving the second-order differential equation $$ a(t) =\frac{d^2y}{dt^2}= -g +Dv^2(t). $$ We rewrite it as two coupled first-order equations (this is a standard approach) $$ \frac{dy}{dt} = v(t), $$ with initial condition \( y(t_0)=y_0 \) and $$ a(t) =\frac{dv}{dt}= -g +Dv^2(t), $$ with initial condition \( v(t_0)=v_0 \).
Many of the algorithms for solving differential equations start with simple Taylor equations. If we now Taylor expand \( y \) and \( v \) around a value \( t+\Delta t \) we have $$ y(t+\Delta t) = y(t)+\Delta t \frac{dy}{dt}+\frac{\Delta t^2}{2!} \frac{d^2y}{dt^2}+O(\Delta t^3), $$ and $$ v(t+\Delta t) = v(t)+\Delta t \frac{dv}{dt}+\frac{\Delta t^2}{2!} \frac{d^2v}{dt^2}+O(\Delta t^3). $$ Using the fact that \( dy/dt = v \) and \( dv/dt=a \) and keeping only terms up to \( \Delta t \) we have $$ y(t+\Delta t) = y(t)+\Delta t v(t)+O(\Delta t^2), $$ and $$ v(t+\Delta t) = v(t)+\Delta t a(t)+O(\Delta t^2). $$
Using our discretized versions of the equations with for example \( y_{i}=y(t_i) \) and \( y_{i\pm 1}=y(t_i+\Delta t) \), we can rewrite the above equations as (and truncating at \( \Delta t \)) $$ y_{i+1} = y_i+\Delta t v_i, $$ and $$ v_{i+1} = v_i+\Delta t a_i. $$ These are the famous Euler equations (forward Euler).
To solve these equations numerically we start at a time \( t_0 \) and simply integrate up these equations to a final time \( t_f \), The step size \( \Delta t \) is an input parameter in our code. You can define it directly in the code below as
DeltaT = 0.1
With a given final time tfinal we can then find the number of integration points via the ceil function included in the math package of Python as
#define final time, assuming that initial time is zero
from math import ceil
tfinal = 0.5
n = ceil(tfinal/DeltaT)
print(n)
The ceil function returns the smallest integer not less than the input in say
x = 21.15
print(ceil(x))
which in the case here is 22.
x = 21.75
print(ceil(x))
which also yields 22. The floor function in the math package is used to return the closest integer value which is less than or equal to the specified expression or value. Compare the previous result to the usage of floor
from math import floor
x = 21.75
print(floor(x))
Alternatively, we can define ourselves the number of integration(mesh) points. In this case we could have
n = 10
tinitial = 0.0
tfinal = 0.5
DeltaT = (tfinal-tinitial)/(n)
print(DeltaT)
Since we will set up one-dimensional arrays that contain the values of various variables like time, position, velocity, acceleration etc, we need to know the value of \( n \), the number of data points (or integration or mesh points). With \( n \) we can initialize a given array by setting all elelements to zero, as done here
# define array a
a = np.zeros(n)
print(a)
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
g = 9.80655 #m/s^2
D = 0.00245 #m/s
DeltaT = 0.1
#set up arrays
tfinal = 0.5
n = ceil(tfinal/DeltaT)
# define scaling constant vT
vT = sqrt(g/D)
# set up arrays for t, a, v, and y and we can compare our results with analytical ones
t = np.zeros(n)
a = np.zeros(n)
v = np.zeros(n)
y = np.zeros(n)
yanalytic = np.zeros(n)
# Initial conditions
v[0] = 0.0 #m/s
y[0] = 10.0 #m
yanalytic[0] = y[0]
# Start integrating using Euler's method
for i in range(n-1):
# expression for acceleration
a[i] = -g + D*v[i]*v[i]
# update velocity and position
y[i+1] = y[i] + DeltaT*v[i]
v[i+1] = v[i] + DeltaT*a[i]
# update time to next time step and compute analytical answer
t[i+1] = t[i] + DeltaT
yanalytic[i+1] = y[0]-(vT*vT/g)*log(cosh(g*t[i+1]/vT))
if ( y[i+1] < 0.0):
break
a[n-1] = -g + D*v[n-1]*v[n-1]
data = {'t[s]': t,
'y[m]': y-yanalytic,
'v[m/s]': v,
'a[m/s^2]': a
}
NewData = pd.DataFrame(data)
display(NewData)
#finally we plot the data
fig, axs = plt.subplots(3, 1)
axs[0].plot(t, y, t, yanalytic)
axs[0].set_xlim(0, tfinal)
axs[0].set_ylabel('y and exact')
axs[1].plot(t, v)
axs[1].set_ylabel('v[m/s]')
axs[2].plot(t, a)
axs[2].set_xlabel('time[s]')
axs[2].set_ylabel('a[m/s^2]')
fig.tight_layout()
save_fig("EulerIntegration")
plt.show()
Try different values for \( \Delta t \) and study the difference between the exact solution and the numerical solution.
The Euler-Cromer method is a simple variant of the standard Euler method. We use the newly updated velocity \( v_{i+1} \) as an input to the new position, that is, instead of $$ y_{i+1} = y_i+\Delta t v_i, $$ and $$ v_{i+1} = v_i+\Delta t a_i, $$ we use now the newly calculate for \( v_{i+1} \) as input to \( y_{i+1} \), that is we compute first $$ v_{i+1} = v_i+\Delta t a_i, $$ and then $$ y_{i+1} = y_i+\Delta t v_{i+1}, $$
Implementing the Euler-Cromer method yields a simple change to the previous code. We only need to change the following line in the loop over time steps
for i in range(n-1):
# more codes in between here
v[i+1] = v[i] + DeltaT*a[i]
y[i+1] = y[i] + DeltaT*v[i+1]
# more code
We will make extensive use of Python as programming language and its myriad of available libraries. You will find Jupyter notebooks invaluable in your work.
If you have Python installed (we strongly recommend Python3) and you feel pretty familiar with installing different packages, we recommend that you install the following Python packages via pip as
For OSX users we recommend, after having installed Xcode, to install brew. Brew allows for a seamless installation of additional software via for example
If you don't want to perform these operations separately and venture into the hassle of exploring how to set up dependencies and paths, we recommend two widely used distrubutions which set up all relevant dependencies for Python, namely
which is an open source distribution of the Python and R programming languages for large-scale data processing, predictive analytics, and scientific computing, that aims to simplify package management and deployment. Package versions are managed by the package management system conda. is a Python distribution for scientific and analytic computing distribution and analysis environment, available for free and under a commercial license.Furthermore, Google's Colab is a free Jupyter notebook environment that requires no setup and runs entirely in the cloud. Try it out!
pycod jupyter nbconvert filename.ipynb --to latex
And to add more versatility, the Python package SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) and is entirely written in Python.
There are several central software libraries for linear algebra and eigenvalue problems. Several of the more popular ones have been wrapped into ofter software packages like those from the widely used text Numerical Recipes. The original source codes in many of the available packages are often taken from the widely used software package LAPACK, which follows two other popular packages developed in the 1970s, namely EISPACK and LINPACK. We describe them shortly here.
$$ \mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix}\qquad \mathbf{I} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} $$
The inverse of a matrix is defined by $$ \mathbf{A}^{-1} \cdot \mathbf{A} = I $$
| Relations | Name | matrix elements |
|---|---|---|
| \( A = A^{T} \) | symmetric | \( a_{ij} = a_{ji} \) |
| \( A = \left (A^{T} \right )^{-1} \) | real orthogonal | \( \sum_k a_{ik} a_{jk} = \sum_k a_{ki} a_{kj} = \delta_{ij} \) |
| \( A = A^{ * } \) | real matrix | \( a_{ij} = a_{ij}^{ * } \) |
| \( A = A^{\dagger} \) | hermitian | \( a_{ij} = a_{ji}^{ * } \) |
| \( A = \left (A^{\dagger} \right )^{-1} \) | unitary | \( \sum_k a_{ik} a_{jk}^{ * } = \sum_k a_{ki}^{ * } a_{kj} = \delta_{ij} \) |
For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all equivalent
import numpy as np
Here follows a simple example where we set up an array of ten elements, all determined by random numbers drawn according to the normal distribution,
n = 10
x = np.random.normal(size=n)
print(x)
We defined a vector \( x \) with \( n=10 \) elements with its values given by the Normal distribution \( N(0,1) \). Another alternative is to declare a vector as follows
import numpy as np
x = np.array([1, 2, 3])
print(x)
Here we have defined a vector with three elements, with \( x_0=1 \), \( x_1=2 \) and \( x_2=3 \). Note that both Python and C++ start numbering array elements from \( 0 \) and on. This means that a vector with \( n \) elements has a sequence of entities \( x_0, x_1, x_2, \dots, x_{n-1} \). We could also let (recommended) Numpy to compute the logarithms of a specific array as
import numpy as np
x = np.log(np.array([4, 7, 8]))
print(x)
In the last example we used Numpy's unary function \( np.log \). This function is highly tuned to compute array elements since the code is vectorized and does not require looping. We normaly recommend that you use the Numpy intrinsic functions instead of the corresponding log function from Python's math module. The looping is done explicitely by the np.log function. The alternative, and slower way to compute the logarithms of a vector would be to write
import numpy as np
from math import log
x = np.array([4, 7, 8])
for i in range(0, len(x)):
x[i] = log(x[i])
print(x)
We note that our code is much longer already and we need to import the log function from the math module. The attentive reader will also notice that the output is \( [1, 1, 2] \). Python interprets automagically our numbers as integers (like the automatic keyword in C++). To change this we could define our array elements to be double precision numbers as
import numpy as np
x = np.log(np.array([4, 7, 8], dtype = np.float64))
print(x)
or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is
import numpy as np
x = np.log(np.array([4.0, 7.0, 8.0])
print(x)
To check the number of bytes (remember that one byte contains eight bits for double precision variables), you can use simple use the itemsize functionality (the array \( x \) is actually an object which inherits the functionalities defined in Numpy) as
import numpy as np
x = np.log(np.array([4.0, 7.0, 8.0])
print(x.itemsize)
Having defined vectors, we are now ready to try out matrices. We can define a \( 3 \times 3 \) real matrix \( \hat{A} \) as (recall that we user lowercase letters for vectors and uppercase letters for matrices)
import numpy as np
A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
print(A)
If we use the shape function we would get \( (3, 3) \) as output, that is verifying that our matrix is a \( 3\times 3 \) matrix. We can slice the matrix and print for example the first column (Python organized matrix elements in a row-major order, see below) as
import numpy as np
A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
# print the first column, row-major order and elements start with 0
print(A[:,0])
We can continue this was by printing out other columns or rows. The example here prints out the second column
import numpy as np
A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
# print the first column, row-major order and elements start with 0
print(A[1,:])
Numpy contains many other functionalities that allow us to slice, subdivide etc etc arrays. We strongly recommend that you look up the Numpy website for more details. Useful functions when defining a matrix are the np.zeros function which declares a matrix of a given dimension and sets all elements to zero
import numpy as np
n = 10
# define a matrix of dimension 10 x 10 and set all elements to zero
A = np.zeros( (n, n) )
print(A)
or initializing all elements to
import numpy as np
n = 10
# define a matrix of dimension 10 x 10 and set all elements to one
A = np.ones( (n, n) )
print(A)
or as unitarily distributed random numbers (see the material on random number generators in the statistics part)
import numpy as np
n = 10
# define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]
A = np.random.rand(n, n)
print(A)
Another useful Python package is pandas, which is an open source library providing high-performance, easy-to-use data structures and data analysis tools for Python. pandas stands for panel data, a term borrowed from econometrics and is an efficient library for data analysis with an emphasis on tabular data. pandas has two major classes, the DataFrame class with two-dimensional data objects and tabular data organized in columns and the class Series with a focus on one-dimensional data objects. Both classes allow you to index data easily as we will see in the examples below. pandas allows you also to perform mathematical operations on the data, spanning from simple reshapings of vectors and matrices to statistical operations.
The following simple example shows how we can, in an easy way make tables of our data. Here we define a data set which includes names, place of birth and date of birth, and displays the data in an easy to read way. We will see repeated use of pandas, in particular in connection with classification of data.
import pandas as pd
from IPython.display import display
data = {'First Name': ["Frodo", "Bilbo", "Aragorn II", "Samwise"],
'Last Name': ["Baggins", "Baggins","Elessar","Gamgee"],
'Place of birth': ["Shire", "Shire", "Eriador", "Shire"],
'Date of Birth T.A.': [2968, 2890, 2931, 2980]
}
data_pandas = pd.DataFrame(data)
display(data_pandas)
In the above we have imported pandas with the shorthand pd, the latter has become the standard way we import pandas. We make then a list of various variables and reorganize the above lists into a DataFrame and then print out a neat table with specific column labels as Name, place of birth and date of birth. Displaying these results, we see that the indices are given by the default numbers from zero to three. pandas is extremely flexible and we can easily change the above indices by defining a new type of indexing as
data_pandas = pd.DataFrame(data,index=['Frodo','Bilbo','Aragorn','Sam'])
display(data_pandas)
Thereafter we display the content of the row which begins with the index Aragorn
display(data_pandas.loc['Aragorn'])
We can easily append data to this, for example
new_hobbit = {'First Name': ["Peregrin"],
'Last Name': ["Took"],
'Place of birth': ["Shire"],
'Date of Birth T.A.': [2990]
}
data_pandas=data_pandas.append(pd.DataFrame(new_hobbit, index=['Pippin']))
display(data_pandas)
Here are other examples where we use the DataFrame functionality to handle arrays, now with more interesting features for us, namely numbers. We set up a matrix of dimensionality \( 10\times 5 \) and compute the mean value and standard deviation of each column. Similarly, we can perform mathematial operations like squaring the matrix elements and many other operations.
import numpy as np
import pandas as pd
from IPython.display import display
np.random.seed(100)
# setting up a 10 x 5 matrix
rows = 10
cols = 5
a = np.random.randn(rows,cols)
df = pd.DataFrame(a)
display(df)
print(df.mean())
print(df.std())
display(df**2)
Thereafter we can select specific columns only and plot final results
df.columns = ['First', 'Second', 'Third', 'Fourth', 'Fifth']
df.index = np.arange(10)
display(df)
print(df['Second'].mean() )
print(df.info())
print(df.describe())
from pylab import plt, mpl
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
df.cumsum().plot(lw=2.0, figsize=(10,6))
plt.show()
df.plot.bar(figsize=(10,6), rot=15)
plt.show()
We can produce a \( 4\times 4 \) matrix
b = np.arange(16).reshape((4,4))
print(b)
df1 = pd.DataFrame(b)
print(df1)
and many other operations.
The Series class is another important class included in pandas. You can view it as a specialization of DataFrame but where we have just a single column of data. It shares many of the same features as _DataFrame. As with DataFrame, most operations are vectorized, achieving thereby a high performance when dealing with computations of arrays, in particular labeled arrays. As we will see below it leads also to a very concice code close to the mathematical operations we may be interested in. For multidimensional arrays, we recommend strongly xarray. xarray has much of the same flexibility as pandas, but allows for the extension to higher dimensions than two.
An overarching aim in this course is to give you a deeper understanding of the scientific method. The problems we study will all involve cases where we can apply classical mechanics. In our previous material we already assumed that we had a model for the motion of an object. Alternatively we could have data from experiment (like Usain Bolt's 100m world record run in 2008). Or we could have performed ourselves an experiment and we want to understand which forces are at play and whether these forces can be understood in terms of fundamental forces.
Our first step consists in identifying the problem. What we sketch here may include a mix of experiment and theoretical simulations, or just experiment or only theory.
Here we can ask questions like
With our eventual data and observations we would now like to develop a model for the system. In the end we want obviously to be able to understand which forces are at play and how they influence our specific system. That is, can we extract some deeper insights about a system?
We need then to
With the model at hand, we can then solve the equations. In classical mechanics we normally end up with solving sets of coupled ordinary differential equations or partial differential equations.
Now it is time to ask the big questions. What do our results mean? Can we give a simple interpretation in terms of fundamental laws? What do our results mean? Are they correct? Thus, typical questions we may ask are
All these steps are essential elements in a scientific enquiry. Hopefully, through a mix of numerical simulations, analytical calculations and experiments we may gain a deeper insight about the physics of a specific system.
Let us now remind ourselves of Newton's laws, since these are the laws of motion we will study in this course.
When analyzing a physical system we normally start with distinguishing between the object we are studying (we will label this in more general terms as our system) and how this system interacts with the environment (which often means everything else!)
In our investigations we will thus analyze a specific physics problem in terms of the system and the environment. In doing so we need to identify the forces that act on the system and assume that the forces acting on the system must have a source, an identifiable cause in the environment.
A force acting on for example a falling object must be related to an interaction with something in the environment. This also means that we do not consider internal forces. The latter are forces between one part of the object and another part. In this course we will mainly focus on external forces.
Forces are either contact forces or long-range forces.
Contact forces, as evident from the name, are forces that occur at the contact between the system and the environment. Well-known long-range forces are the gravitional force and the electromagnetic force.
In order to set up the forces which act on an object, the following steps may be useful
Newton’s second law of motion: The force \( \boldsymbol{F} \) on an object of inertial mass \( m \) is related to the acceleration a of the object through $$ \boldsymbol{F} = m\boldsymbol{a}, $$ where \( \boldsymbol{a} \) is the acceleration.
Newton’s laws of motion are laws of nature that have been found by experimental investigations and have been shown to hold up to continued experimental investigations. Newton’s laws are valid over a wide range of length- and time-scales. We use Newton’s laws of motion to describe everything from the motion of atoms to the motion of galaxies.
The second law is a vector equation with the acceleration having the same direction as the force. The acceleration is proportional to the force via the mass \( m \) of the system under study.
Newton’s second law introduces a new property of an object, the so-called inertial mass \( m \). We determine the inertial mass of an object by measuring the acceleration for a given applied force.
What happens if the net external force on a body is zero? Applying Newton’s second law, we find: $$ \boldsymbol{F} = 0 = m\boldsymbol{a}, $$ which gives using the definition of the acceleration $$ \boldsymbol{a} = \frac{d\boldsymbol{v}}{dt}=0. $$ The acceleration is zero, which means that the velocity of the object is constant. This is often referred to as Newton’s first law. An object in a state of uniform motion tends to remain in that state unless an external force changes its state of motion. Why do we need a separate law for this? Is it not simply a special case of Newton’s second law? Yes, Newton’s first law can be deduced from the second law as we have illustrated. However, the first law is often used for a different purpose: Newton’s First Law tells us about the limit of applicability of Newton’s Second law. Newton’s Second law can only be used in reference systems where the First law is obeyed. But is not the First law always valid? No! The First law is only valid in reference systems that are not accelerated. If you observe the motion of a ball from an accelerating car, the ball will appear to accelerate even if there are no forces acting on it. We call systems that are not accelerating inertial systems, and Newton’s first law is often called the law of inertia. Newton’s first and second laws of motion are only valid in inertial systems.
A system is an inertial system if it is not accelerated. It means that the reference system must not be accelerating linearly or rotating. Unfortunately, this means that most systems we know are not really inertial systems. For example, the surface of the Earth is clearly not an inertial system, because the Earth is rotating. The Earth is also not an inertial system, because it ismoving in a curved path around the Sun. However, even if the surface of the Earth is not strictly an inertial system, it may be considered to be approximately an inertial system for many laboratory-size experiments.
If there is a force from object A on object B, there is also a force from object B on object A. This fundamental principle of interactions is called Newton’s third law. We do not know of any force that do not obey this law: All forces appear in pairs. Newton’s third law is usually formulated as: For every action there is an equal and opposite reaction.
Here we consider the motion of a single particle moving under the influence of some set of forces. We will consider some problems where the force does not depend on the position. In that case Newton's law \( m\dot{\boldsymbol{v}}=\boldsymbol{F}(\boldsymbol{v}) \) is a first-order differential equation and one solves for \( \boldsymbol{v}(t) \), then moves on to integrate \( \boldsymbol{v} \) to get the position. In essentially all of these cases we cna find an analytical solution.
Air resistance tends to scale as the square of the velocity. This is in contrast to many problems chosen for textbooks, where it is linear in the velocity. The choice of a linear dependence is motivated by mathematical simplicity (it keeps the differential equation linear) rather than by physics. One can see that the force should be quadratic in velocity by considering the momentum imparted on the air molecules. If an object sweeps through a volume \( dV \) of air in time \( dt \), the momentum imparted on the air is $$ \begin{equation} dP=\rho_m dV v, \label{_auto1} \end{equation} $$
where \( v \) is the velocity of the object and \( \rho_m \) is the mass density of the air. If the molecules bounce back as opposed to stop you would double the size of the term. The opposite value of the momentum is imparted onto the object itself. Geometrically, the differential volume is $$ \begin{equation} dV=Avdt, \label{_auto2} \end{equation} $$
where \( A \) is the cross-sectional area and \( vdt \) is the distance the object moved in time \( dt \).
This is the force felt by the particle, and is opposite to its direction of motion. Now, because air doesn't stop when it hits an object, but flows around the best it can, the actual force is reduced by a dimensionless factor \( c_W \), called the drag coefficient. $$ \begin{equation} F_{\rm drag}=-c_W\rho_m Av^2, \label{_auto4} \end{equation} $$
and the acceleration is $$ \begin{eqnarray} \frac{dv}{dt}=-\frac{c_W\rho_mA}{m}v^2. \end{eqnarray} $$
For a particle with initial velocity \( v_0 \), one can separate the \( dt \) to one side of the equation, and move everything with $v$s to the other side. We did this in our discussion of simple motion and will not repeat it here.
On more general terms, for many systems, e.g. an automobile, there are multiple sources of resistance. In addition to wind resistance, where the force is proportional to \( v^2 \), there are dissipative effects of the tires on the pavement, and in the axel and drive train. These other forces can have components that scale proportional to \( v \), and components that are independent of \( v \). Those independent of \( v \), e.g. the usual \( f=\mu_K N \) frictional force you consider in your first Physics courses, only set in once the object is actually moving. As speeds become higher, the \( v^2 \) components begin to dominate relative to the others. For automobiles at freeway speeds, the \( v^2 \) terms are largely responsible for the loss of efficiency. To travel a distance \( L \) at fixed speed \( v \), the energy/work required to overcome the dissipative forces are \( fL \), which for a force of the form \( f=\alpha v^n \) becomes $$ \begin{eqnarray} W=\int dx~f=\alpha v^n L. \end{eqnarray} $$ For \( n=0 \) the work is independent of speed, but for the wind resistance, where \( n=2 \), slowing down is essential if one wishes to reduce fuel consumption. It is also important to consider that engines are designed to be most efficient at a chosen range of power output. Thus, some cars will get better mileage at higher speeds (They perform better at 50 mph than at 5 mph) despite the considerations mentioned above.
As an example of Newton's Laws we consider projectile motion (or a falling raindrop or a ball we throw up in the air) with a drag force. Even though air resistance is largely proportional to the square of the velocity, we will consider the drag force to be linear to the velocity, \( \boldsymbol{F}=-m\gamma\boldsymbol{v} \), for the purposes of this exercise. The acceleration for a projectile moving upwards, \( \boldsymbol{a}=\boldsymbol{F}/m \), becomes $$ \begin{eqnarray} \frac{dv_x}{dt}=-\gamma v_x,\\ \nonumber \frac{dv_y}{dt}=-\gamma v_y-g, \end{eqnarray} $$ and \( \gamma \) has dimensions of inverse time.
If you on the other hand have a falling raindrop, how do these equations change? See for example Figure 2.1 in Taylor. Let us stay with a ball which is thrown up in the air at \( t=0 \).
We will go over two different ways to solve this equation. The first by direct integration, and the second as a differential equation. To do this by direct integration, one simply multiplies both sides of the equations above by \( dt \), then divide by the appropriate factors so that the $v$s are all on one side of the equation and the \( dt \) is on the other. For the \( x \) motion one finds an easily integrable equation, $$ \begin{eqnarray} \frac{dv_x}{v_x}&=&-\gamma dt,\\ \nonumber \int_{v_{0x}}^{v_{x}}\frac{dv_x}{v_x}&=&-\gamma\int_0^{t}dt,\\ \nonumber \ln\left(\frac{v_{x}}{v_{0x}}\right)&=&-\gamma t,\\ \nonumber v_{x}(t)&=&v_{0x}e^{-\gamma t}. \end{eqnarray} $$
This is very much the result you would have written down by inspection. For the \( y \)-component of the velocity, $$ \begin{eqnarray} \frac{dv_y}{v_y+g/\gamma}&=&-\gamma dt\\ \nonumber \ln\left(\frac{v_{y}+g/\gamma}{v_{0y}-g/\gamma}\right)&=&-\gamma t_f,\\ \nonumber v_{fy}&=&-\frac{g}{\gamma}+\left(v_{0y}+\frac{g}{\gamma}\right)e^{-\gamma t}. \end{eqnarray} $$
Whereas \( v_x \) starts at some value and decays exponentially to zero, \( v_y \) decays exponentially to the terminal velocity, \( v_t=-g/\gamma \).
Although this direct integration is simpler than the method we invoke below, the method below will come in useful for some slightly more difficult differential equations in the future. The differential equation for \( v_x \) is straight-forward to solve. Because it is first order there is one arbitrary constant, \( A \), and by inspection the solution is $$ \begin{equation} v_x=Ae^{-\gamma t}. \label{_auto5} \end{equation} $$
The arbitrary constants for equations of motion are usually determined by the initial conditions, or more generally boundary conditions. By inspection \( A=v_{0x} \), the initial \( x \) component of the velocity.
The differential equation for \( v_y \) is a bit more complicated due to the presence of \( g \). Differential equations where all the terms are linearly proportional to a function, in this case \( v_y \), or to derivatives of the function, e.g., \( v_y \), \( dv_y/dt \), \( d^2v_y/dt^2\cdots \), are called linear differential equations. If there are terms proportional to \( v^2 \), as would happen if the drag force were proportional to the square of the velocity, the differential equation is not longer linear. Because this expression has only one derivative in \( v \) it is a first-order linear differential equation. If a term were added proportional to \( d^2v/dt^2 \) it would be a second-order differential equation. In this case we have a term completely independent of \( v \), the gravitational acceleration \( g \), and the usual strategy is to first rewrite the equation with all the linear terms on one side of the equal sign, $$ \begin{equation} \frac{dv_y}{dt}+\gamma v_y=-g. \label{_auto6} \end{equation} $$
Now, the solution to the equation can be broken into two parts. Because this is a first-order differential equation we know that there will be one arbitrary constant. Physically, the arbitrary constant will be determined by setting the initial velocity, though it could be determined by setting the velocity at any given time. Like most differential equations, solutions are not "solved". Instead, one guesses at a form, then shows the guess is correct. For these types of equations, one first tries to find a single solution, i.e. one with no arbitrary constants. This is called the {\it particular} solution, \( y_p(t) \), though it should really be called "a" particular solution because there are an infinite number of such solutions. One then finds a solution to the {\it homogenous} equation, which is the equation with zero on the right-hand side, $$ \begin{equation} \frac{dv_{y,h}}{dt}+\gamma v_{y,h}=0. \label{_auto7} \end{equation} $$
Homogenous solutions will have arbitrary constants.
The particular solution will solve the same equation as the original general equation $$ \begin{equation} \frac{dv_{y,p}}{dt}+\gamma v_{y,p}=-g. \label{_auto8} \end{equation} $$
However, we don't need find one with arbitrary constants. Hence, it is called a particular solution.
The sum of the two, $$ \begin{equation} v_y=v_{y,p}+v_{y,h}, \label{_auto9} \end{equation} $$
is a solution of the total equation because of the linear nature of the differential equation. One has now found a general solution encompassing all solutions, because it both satisfies the general equation (like the particular solution), and has an arbitrary constant that can be adjusted to fit any initial condition (like the homogneous solution). If the equation were not linear, e.g if there were a term such as \( v_y^2 \) or \( v_y\dot{v}_y \), this technique would not work.
Returning to the example above, the homogenous solution is the same as that for \( v_x \), because there was no gravitational acceleration in that case, $$ \begin{equation} v_{y,h}=Be^{-\gamma t}. \label{_auto10} \end{equation} $$ In this case a particular solution is one with constant velocity, $$ \begin{equation} v_{y,p}=-g/\gamma. \label{_auto11} \end{equation} $$ Note that this is the terminal velocity of a particle falling from a great height. The general solution is thus, $$ \begin{equation} v_y=Be^{-\gamma t}-g/\gamma, \label{_auto12} \end{equation} $$
and one can find \( B \) from the initial velocity, $$ \begin{equation} v_{0y}=B-g/\gamma,~~~B=v_{0y}+g/\gamma. \label{_auto13} \end{equation} $$ Plugging in the expression for \( B \) gives the \( y \) motion given the initial velocity, $$ \begin{equation} v_y=(v_{0y}+g/\gamma)e^{-\gamma t}-g/\gamma. \label{_auto14} \end{equation} $$ It is easy to see that this solution has \( v_y=v_{0y} \) when \( t=0 \) and \( v_y=-g/\gamma \) when \( t\rightarrow\infty \).
One can also integrate the two equations to find the coordinates \( x \) and \( y \) as functions of \( t \), $$ \begin{eqnarray} x&=&\int_0^t dt'~v_{0x}(t')=\frac{v_{0x}}{\gamma}\left(1-e^{-\gamma t}\right),\\ \nonumber y&=&\int_0^t dt'~v_{0y}(t')=-\frac{gt}{\gamma}+\frac{v_{0y}+g/\gamma}{\gamma}\left(1-e^{-\gamma t}\right). \end{eqnarray} $$
If the question was to find the position at a time \( t \), we would be finished. However, the more common goal in a projectile equation problem is to find the range, i.e. the distance \( x \) at which \( y \) returns to zero. For the case without a drag force this was much simpler. The solution for the \( y \) coordinate would have been \( y=v_{0y}t-gt^2/2 \). One would solve for \( t \) to make \( y=0 \), which would be \( t=2v_{0y}/g \), then plug that value for \( t \) into \( x=v_{0x}t \) to find \( x=2v_{0x}v_{0y}/g=v_0\sin(2\theta_0)/g \). One follows the same steps here, except that the expression for \( y(t) \) is more complicated. Searching for the time where \( y=0 \), and we get $$ \begin{equation} 0=-\frac{gt}{\gamma}+\frac{v_{0y}+g/\gamma}{\gamma}\left(1-e^{-\gamma t}\right). \label{_auto15} \end{equation} $$
This cannot be inverted into a simple expression \( t=\cdots \). Such expressions are known as "transcendental equations", and are not the rare instance, but are the norm. In the days before computers, one might plot the right-hand side of the above graphically as a function of time, then find the point where it crosses zero.
Now, the most common way to solve for an equation of the above type would be to apply Newton's method numerically. This involves the following algorithm for finding solutions of some equation \( F(t)=0 \).
Another example of a velocity-dependent force is magnetism, $$ \begin{eqnarray} \boldsymbol{F}&=&q\boldsymbol{v}\times\boldsymbol{B},\\ \nonumber F_i&=&q\sum_{jk}\epsilon_{ijk}v_jB_k. \end{eqnarray} $$
For a uniform field in the \( z \) direction \( \boldsymbol{B}=B\hat{z} \), the force can only have \( x \) and \( y \) components, $$ \begin{eqnarray} F_x&=&qBv_y\\ \nonumber F_y&=&-qBv_x. \end{eqnarray} $$
The differential equations are $$ \begin{eqnarray} \dot{v}_x&=&\omega_c v_y,\omega_c= qB/m\\ \nonumber \dot{v}_y&=&-\omega_c v_x. \end{eqnarray} $$ One can solve the equations by taking time derivatives of either equation, then substituting into the other equation, $$ \begin{eqnarray} \ddot{v}_x=\omega_c\dot{v_y}=-\omega_c^2v_x,\\ \nonumber \ddot{v}_y&=&-\omega_c\dot{v}_x=-\omega_cv_y. \end{eqnarray} $$
The solution to these equations can be seen by inspection, $$ \begin{eqnarray} v_x&=&A\sin(\omega_ct+\phi),\\ \nonumber v_y&=&A\cos(\omega_ct+\phi). \end{eqnarray} $$
One can integrate the equations to find the positions as a function of time, $$ \begin{eqnarray} x-x_0&=&\int_{x_0}^x dx=\int_0^t dt v(t)\\ \nonumber &=&\frac{-A}{\omega_c}\cos(\omega_ct+\phi),\\ \nonumber y-y_0&=&\frac{A}{\omega_c}\sin(\omega_ct+\phi). \end{eqnarray} $$ The trajectory is a circle centered at \( x_0,y_0 \) with amplitude \( A \) rotating in the clockwise direction.
The equations of motion for the \( z \) motion are $$ \begin{equation} \dot{v_z}=0, \label{_auto16} \end{equation} $$ which leads to $$ \begin{equation} z-z_0=V_zt. \label{_auto17} \end{equation} $$ Added onto the circle, the motion is helical.
Note that the kinetic energy, $$ \begin{equation} T=\frac{1}{2}m(v_x^2+v_y^2+v_z^2)=\frac{1}{2}m(\omega_c^2A^2+V_z^2), \label{_auto18} \end{equation} $$ is constant. This is because the force is perpendicular to the velocity, so that in any differential time element \( dt \) the work done on the particle \( \boldsymbol{F}\cdot{dr}=dt\boldsymbol{F}\cdot{v}=0 \).
One should think about the implications of a velocity dependent force. Suppose one had a constant magnetic field in deep space. If a particle came through with velocity \( v_0 \), it would undergo cyclotron motion with radius \( R=v_0/\omega_c \). However, if it were still its motion would remain fixed. Now, suppose an observer looked at the particle in one reference frame where the particle was moving, then changed their velocity so that the particle's velocity appeared to be zero. The motion would change from circular to fixed. Is this possible?
The solution to the puzzle above relies on understanding relativity. Imagine that the first observer believes \( \boldsymbol{B}\ne 0 \) and that the electric field \( \boldsymbol{E}=0 \). If the observer then changes reference frames by accelerating to a velocity \( \boldsymbol{v} \), in the new frame \( \boldsymbol{B} \) and \( \boldsymbol{E} \) both change. If the observer moved to the frame where the charge, originally moving with a small velocity \( v \), is now at rest, the new electric field is indeed \( \boldsymbol{v}\times\boldsymbol{B} \), which then leads to the same acceleration as one had before. If the velocity is not small compared to the speed of light, additional \( \gamma \) factors come into play, \( \gamma=1/\sqrt{1-(v/c)^2} \). Relativistic motion will not be considered in this course.
Another classical case is that of simple harmonic oscillations, here represented by a block sliding on a horizontal frictionless surface. The block is tied to a wall with a spring. If the spring is not compressed or stretched too far, the force on the block at a given position \( x \) is $$ F=-kx. $$ The negative sign means that the force acts to restore the object to an equilibrium position. Newton's equation of motion for this idealized system is then $$ m\frac{d^2x}{dt^2}=-kx, $$ or we could rephrase it as $$ \frac{d^2x}{dt^2}=-\frac{k}{m}x=-\omega_0^2x, \label{eq:newton1} $$ with the angular frequency \( \omega_0^2=k/m \).
The above differential equation has the advantage that it can be solved analytically with solutions on the form $$ x(t)=Acos(\omega_0t+\nu), $$ where \( A \) is the amplitude and \( \nu \) the phase constant. This provides in turn an important test for the numerical solution and the development of a program for more complicated cases which cannot be solved analytically.
With the position \( x(t) \) and the velocity \( v(t)=dx/dt \) we can reformulate Newton's equation in the following way $$ \frac{dx(t)}{dt}=v(t), $$ and $$ \frac{dv(t)}{dt}=-\omega_0^2x(t). $$ We are now going to solve these equations using first the standard forward Euler method. Later we will try to improve upon this.
Before proceeding however, it is important to note that in addition to the exact solution, we have at least two further tests which can be used to check our solution.
Since functions like \( cos \) are periodic with a period \( 2\pi \), then the solution \( x(t) \) has also to be periodic. This means that $$ x(t+T)=x(t), $$ with \( T \) the period defined as $$ T=\frac{2\pi}{\omega_0}=\frac{2\pi}{\sqrt{k/m}}. $$ Observe that \( T \) depends only on \( k/m \) and not on the amplitude of the solution.
In addition to the periodicity test, the total energy has also to be conserved.
Suppose we choose the initial conditions $$ x(t=0)=1\hspace{0.1cm} \mathrm{m}\hspace{1cm} v(t=0)=0\hspace{0.1cm}\mathrm{m/s}, $$ meaning that block is at rest at \( t=0 \) but with a potential energy $$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$ The total energy at any time \( t \) has however to be conserved, meaning that our solution has to fulfil the condition $$ E_0=\frac{1}{2}kx(t)^2+\frac{1}{2}mv(t)^2. $$ We will derive this equation in our discussion on energy conservation.
An algorithm which implements these equations is included below.
#
# This program solves Newtons equation for a block sliding on
# an horizontal frictionless surface.
# The block is tied to the wall with a spring, so N's eq takes the form:
#
# m d^2x/dt^2 = - kx
#
# In order to make the solution dimless, we set k/m = 1.
# This results in two coupled diff. eq's that may be written as:
#
# dx/dt = v
# dv/dt = -x
#
# The user has to specify the initial velocity and position,
# and the number of steps. The time interval is fixed to
# t \in [0, 4\pi) (two periods)
#
The angular equation of motion of the pendulum is given by Newton's equation and with no external force it reads $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+mgsin(\theta)=0, \label{_auto19} \end{equation} $$ with an angular velocity and acceleration given by $$ \begin{equation} v=l\frac{d\theta}{dt}, \label{_auto20} \end{equation} $$ and $$ \begin{equation} a=l\frac{d^2\theta}{dt^2}. \label{_auto21} \end{equation} $$
We do however expect that the motion will gradually come to an end due a viscous drag torque acting on the pendulum. In the presence of the drag, the above equation becomes $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=0, \label{eq:pend1} \end{equation} $$ where \( \nu \) is now a positive constant parameterizing the viscosity of the medium in question. In order to maintain the motion against viscosity, it is necessary to add some external driving force. We choose here a periodic driving force. The last equation becomes then $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Asin(\omega t), \label{eq:pend2} \end{equation} $$ with \( A \) and \( \omega \) two constants representing the amplitude and the angular frequency respectively. The latter is called the driving frequency.
We define $$ \omega_0=\sqrt{g/l}, $$ the so-called natural frequency and the new dimensionless quantities $$ \hat{t}=\omega_0t, $$ with the dimensionless driving frequency $$ \hat{\omega}=\frac{\omega}{\omega_0}, $$ and introducing the quantity \( Q \), called the quality factor, $$ Q=\frac{mg}{\omega_0\nu}, $$ and the dimensionless amplitude $$ \hat{A}=\frac{A}{mg} $$
We have $$ \frac{d^2\theta}{d\hat{t}^2}+\frac{1}{Q}\frac{d\theta}{d\hat{t}} +sin(\theta)=\hat{A}cos(\hat{\omega}\hat{t}). $$ This equation can in turn be recast in terms of two coupled first-order differential equations as follows $$ \frac{d\theta}{d\hat{t}}=\hat{v}, $$ and $$ \frac{d\hat{v}}{d\hat{t}}=-\frac{\hat{v}}{Q}-sin(\theta)+\hat{A}cos(\hat{\omega}\hat{t}). $$ These are the equations to be solved. The factor \( Q \) represents the number of oscillations of the undriven system that must occur before its energy is significantly reduced due to the viscous drag. The amplitude \( \hat{A} \) is measured in units of the maximum possible gravitational torque while \( \hat{\omega} \) is the angular frequency of the external torque measured in units of the pendulum's natural frequency.
Energy conservation is most convenient as a strategy for addressing problems where time does not appear. For example, a particle goes from position \( x_0 \) with speed \( v_0 \), to position \( x_f \); what is its new speed? However, it can also be applied to problems where time does appear, such as in solving for the trajectory \( x(t) \), or equivalently \( t(x) \).
More material to be added here.
The net energy, \( E=V+K \) where \( K \) is the kinetic energy, is then conserved, $$ \begin{eqnarray} \frac{d}{dt}(K+V)&=&\frac{d}{dt}\left(\frac{m}{2}(v_x^2+v_y^2+v_z^2)+V(\boldsymbol{r})\right)\\ \nonumber &=&m\left(v_x\frac{dv_x}{dt}+v_y\frac{dv_y}{dt}+v_z\frac{dv_z}{dt}\right) +\partial_xV\frac{dx}{dt}+\partial_yV\frac{dy}{dt}+\partial_zV\frac{dz}{dt}\\ \nonumber &=&v_xF_x+v_yF_y+v_zF_z-F_xv_x-F_yv_y-F_zv_z=0. \end{eqnarray} $$
The same proof can be written more compactly with vector notation, $$ \begin{eqnarray} \frac{d}{dt}\left(\frac{m}{2}v^2+V(\boldsymbol{r})\right) &=&m\boldsymbol{v}\cdot\dot{\boldsymbol{v}}+\nabla V(\boldsymbol{r})\cdot\dot{\boldsymbol{r}}\\ \nonumber &=&\boldsymbol{v}\cdot\boldsymbol{F}-\boldsymbol{F}\cdot\boldsymbol{v}=0. \end{eqnarray} $$
Inverting the expression for kinetic energy, $$ \begin{equation} v=\sqrt{2K/m}=\sqrt{2(E-V)/m}, \label{_auto23} \end{equation} $$
allows one to solve for the one-dimensional trajectory \( x(t) \), by finding \( t(x) \), $$ \begin{equation} t=\int_{x_0}^x \frac{dx'}{v(x')}=\int_{x_0}^x\frac{dx'}{\sqrt{2(E-V(x'))/m}}. \label{_auto24} \end{equation} $$
Note this would be much more difficult in higher dimensions, because you would have to determine which points, \( x,y,z \), the particles might reach in the trajectory, whereas in one dimension you can typically tell by simply seeing whether the kinetic energy is positive at every point between the old position and the new position.
Consider a simple harmonic oscillator potential, \( V(x)=kx^2/2 \), with a particle emitted from \( x=0 \) with velocity \( v_0 \). Solve for the trajectory \( t(x) \), $$ \begin{eqnarray} t&=&\int_{0}^x \frac{dx'}{\sqrt{2(E-kx^2/2)/m}}\\ \nonumber &=&\sqrt{m/k}\int_0^x~\frac{dx'}{\sqrt{x_{\rm max}^2-x^{\prime 2}}},~~~x_{\rm max}^2=2E/k. \end{eqnarray} $$
Here \( E=mv_0^2/2 \) and \( x_{\rm max} \) is defined as the maximum displacement before the particle turns around. This integral is done by the substitution \( \sin\theta=x/x_{\rm max} \). $$ \begin{eqnarray} (k/m)^{1/2}t&=&\sin^{-1}(x/x_{\rm max}),\\ \nonumber x&=&x_{\rm max}\sin\omega t,~~~\omega=\sqrt{k/m}. \end{eqnarray} $$
Newton's third law which we met earlier states that For every action there is an equal and opposite reaction, is more accurately stated as If two bodies exert forces on each other, these forces are equal in magnitude and opposite in direction.
This means that for two bodies \( i \) and \( j \), if the force on \( i \) due to \( j \) is called \( \boldsymbol{F}_{ij} \), then $$ \begin{equation} \boldsymbol{F}_{ij}=-\boldsymbol{F}_{ji}. \label{_auto25} \end{equation} $$ Newton's second law, \( \boldsymbol{F}=m\boldsymbol{a} \), can be written for a particle \( i \) as $$ \begin{equation} \boldsymbol{F}_i=\sum_{j\ne i} \boldsymbol{F}_{ij}=m_i\boldsymbol{a}_i, \label{_auto26} \end{equation} $$
where \( \boldsymbol{F}_i \) (a single subscript) denotes the net force acting on \( i \). Because the mass of \( i \) is fixed, one can see that $$ \begin{equation} \boldsymbol{F}_i=\frac{d}{dt}m_i\boldsymbol{v}_i=\sum_{j\ne i}\boldsymbol{F}_{ij}. \label{_auto27} \end{equation} $$
Now, one can sum over all the particles and obtain $$ \begin{eqnarray} \frac{d}{dt}\sum_i m_iv_i&=&\sum_{ij, i\ne j}\boldsymbol{F}_{ij}\\ \nonumber &=&0. \end{eqnarray} $$
The last step made use of the fact that for every term \( ij \), there is an equivalent term \( ji \) with opposite force. Because the momentum is defined as \( m\boldsymbol{v} \), for a system of particles, $$ \begin{equation} \frac{d}{dt}\sum_im_i\boldsymbol{v}_i=0,~~{\rm for~isolated~particles}. \label{_auto28} \end{equation} $$
By "isolated" one means that the only force acting on any particle \( i \) are those originating from other particles in the sum, i.e. "no external" forces. Thus, Newton's third law leads to the conservation of total momentum, $$ \begin{eqnarray} \boldsymbol{P}&=&\sum_i m_i\boldsymbol{v}_i,\\ \nonumber \frac{d}{dt}\boldsymbol{P}&=&0. \end{eqnarray} $$
Consider the rocket of mass \( M \) moving with velocity \( v \). After a brief instant, the velocity of the rocket is \( v+\Delta v \) and the mass is \( M-\Delta M \). Momentum conservation gives $$ \begin{eqnarray*} Mv&=&(M-\Delta M)(v+\Delta v)+\Delta M(v-v_e)\\ 0&=&-\Delta Mv+M\Delta v+\Delta M(v-v_e),\\ 0&=&M\Delta v-\Delta Mv_e. \end{eqnarray*} $$
In the second step we ignored the term \( \Delta M\Delta v \) because it is doubly small. The last equation gives $$ \begin{eqnarray} \Delta v&=&\frac{v_e}{M}\Delta M,\\ \nonumber \frac{dv}{dt}&=&\frac{v_e}{M}\frac{dM}{dt}. \end{eqnarray} $$
Integrating the expression with lower limits \( v_0=0 \) and \( M_0 \), one finds $$ \begin{eqnarray*} v&=&v_e\int_{M_0}^M \frac{dM'}{M'}\\ v&=&-v_e\ln(M/M_0)\\ &=&-v_e\ln[(M_0-\alpha t)/M_0]. \end{eqnarray*} $$
Because the total momentum of an isolated system is constant, one can also quickly see that the center of mass of an isolated system is also constant. The center of mass is the average position of a set of masses weighted by the mass, $$ \begin{equation} \bar{x}=\frac{\sum_im_ix_i}{\sum_i m_i}. \label{_auto29} \end{equation} $$
The rate of change of \( \bar{x} \) is $$ \begin{eqnarray} \dot{\bar{x}}&=&\frac{1}{M}\sum_i m_i\dot{x}_i=\frac{1}{M}P_x. \end{eqnarray} $$
Thus if the total momentum is constant the center of mass moves at a constant velocity, and if the total momentum is zero the center of mass is fixed.
Consider a case where the force always points radially, $$ \begin{equation} \boldsymbol{F}(\boldsymbol{r})=F(r)\hat{r}, \label{_auto30} \end{equation} $$
where \( \hat{r} \) is a unit vector pointing outward from the origin. The angular momentum is defined as $$ \begin{equation} \boldsymbol{L}=\boldsymbol{r}\times\boldsymbol{p}=m\boldsymbol{r}\times\boldsymbol{v}. \label{_auto31} \end{equation} $$
The rate of change of the angular momentum is $$ \begin{eqnarray} \frac{d\boldsymbol{L}}{dt}&=&m\boldsymbol{v}\times\boldsymbol{v}+m\boldsymbol{r}\times\dot{\boldsymbol{v}}\\ \nonumber &=&m\boldsymbol{v}\times\boldsymbol{v}+\boldsymbol{r}\times{\boldsymbol{F}}=0. \end{eqnarray} $$
The first term is zero because \( \boldsymbol{v} \) is parallel to itself, and the second term is zero because \( \boldsymbol{F} \) is parallel to \( \boldsymbol{r} \).
As an aside, one can see from the Levi-Civita symbol that the cross product of a vector with itself is zero. Here, we consider a vector $$ \begin{eqnarray} \boldsymbol{V}&=&\boldsymbol{A}\times\boldsymbol{A},\\ \nonumber V_i&=&(\boldsymbol{A}\times\boldsymbol{A})_i=\sum_{jk}\epsilon_{ijk}A_jA_k. \end{eqnarray} $$
For any term \( i \), there are two contributions. For example, for \( i \) denoting the \( x \) direction, either \( j \) denotes the \( y \) direction and \( k \) denotes the \( z \) direction, or vice versa, so $$ \begin{equation} V_1=\epsilon_{123}A_2A_3+\epsilon_{132}A_3A_2. \label{_auto32} \end{equation} $$
This is zero by the antisymmetry of \( \epsilon \) under permutations.
If the force is not radial, \( \boldsymbol{r}\times\boldsymbol{F}\ne 0 \) as above, and angular momentum is no longer conserved, $$ \begin{equation} \frac{d\boldsymbol{L}}{dt}=\boldsymbol{r}\times\boldsymbol{F}\equiv\boldsymbol{\tau}, \label{_auto33} \end{equation} $$
where \( \boldsymbol{\tau} \) is the torque.
For a system of isolated particles, one can write $$ \begin{eqnarray} \frac{d}{dt}\sum_i\boldsymbol{L}_i&=&\sum_{i\ne j}\boldsymbol{r}_i\times \boldsymbol{F}_{ij}\\ \nonumber &=&\frac{1}{2}\sum_{i\ne j} \boldsymbol{r}_i\times \boldsymbol{F}_{ij}+\boldsymbol{r}_j\times\boldsymbol{F}_{ji}\\ \nonumber &=&\frac{1}{2}\sum_{i\ne j} (\boldsymbol{r}_i-\boldsymbol{r}_j)\times\boldsymbol{F}_{ij}=0, \end{eqnarray} $$
where the last step used Newton's third law, \( \boldsymbol{F}_{ij}=-\boldsymbol{F}_{ji} \). If the forces between the particles are radial, i.e. \( \boldsymbol{F}_{ij} ~||~ (\boldsymbol{r}_i-\boldsymbol{r}_j) \), then each term in the sum is zero and the net angular momentum is fixed. Otherwise, you could imagine an isolated system that would start spinning spontaneously.
One can write the torque about a given axis, which we will denote as \( \hat{z} \), in polar coordinates, where $$ \begin{eqnarray} x&=&r\sin\theta\cos\phi,~~y=r\sin\theta\cos\phi,~~z=r\cos\theta, \end{eqnarray} $$ to find the \( z \) component of the torque, $$ \begin{eqnarray} \tau_z&=&xF_y-yF_x\\ \nonumber &=&-r\sin\theta\left\{\cos\phi \partial_y-\sin\phi \partial_x\right\}V(x,y,z). \end{eqnarray} $$
One can use the chain rule to write the partial derivative w.r.t. \( \phi \) (keeping \( r \) and \( \theta \) fixed), $$ \begin{eqnarray} \partial_\phi&=&\frac{\partial x}{\partial\phi}\partial_x+\frac{\partial_y}{\partial\phi}\partial_y +\frac{\partial z}{\partial\phi}\partial_z\\ \nonumber &=&-r\sin\theta\sin\phi\partial_x+\sin\theta\cos\phi\partial_y. \end{eqnarray} $$
Combining the two equations, $$ \begin{eqnarray} \tau_z&=&-\partial_\phi V(r,\theta,\phi). \end{eqnarray} $$
Thus, if the potential is independent of the azimuthal angle \( \phi \), there is no torque about the \( z \) axis and \( L_z \) is conserved.
When we derived the conservation of energy, we assumed that the potential depended only on position, not on time. If it depended explicitly on time, one can quickly see that the energy would have changed at a rate \( \partial_tV(x,y,z,t) \). Note that if there is no explicit dependence on time, i.e. \( V(x,y,z) \), the potential energy can depend on time through the variations of \( x,y,z \) with time. However, that variation does not lead to energy non-conservation. Further, we just saw that if a potential does not depend on the azimuthal angle about some axis, \( \phi \), that the angular momentum about that axis is conserved.
Now, we relate momentum conservation to translational invariance. Considering a system of particles with positions, \( \boldsymbol{r}_i \), if one changed the coordinate system by a translation by a differential distance \( \boldsymbol{\epsilon} \), the net potential would change by $$ \begin{eqnarray} \delta V(\boldsymbol{r}_1,\boldsymbol{r}_2\cdots)&=&\sum_i \boldsymbol{\epsilon}\cdot\nabla_i V(\boldsymbol{r}_1,\boldsymbol{r}_2,\cdots)\\ \nonumber &=&-\sum_i \boldsymbol{\epsilon}\cdot\boldsymbol{F}_i\\ \nonumber &=&-\frac{d}{dt}\sum_i \boldsymbol{\epsilon}\cdot\boldsymbol{p}_i. \end{eqnarray} $$
Thus, if the potential is unchanged by a translation of the coordinate system, the total momentum is conserved. If the potential is translationally invariant in a given direction, defined by a unit vector, \( \hat{\epsilon} \) in the \( \boldsymbol{\epsilon} \) direction, one can see that $$ \begin{eqnarray} \hat{\epsilon}\cdot\nabla_i V(\boldsymbol{r}_i)&=&0. \end{eqnarray} $$
The component of the total momentum along that axis is conserved. This is rather obvious for a single particle. If \( V(\boldsymbol{r}) \) does not depend on some coordinate \( x \), then the force in the \( x \) direction is \( F_x=-\partial_xV=0 \), and momentum along the \( x \) direction is constant.
We showed how the total momentum of an isolated system of particle was conserved, even if the particles feel internal forces in all directions. In that case the potential energy could be written $$ \begin{eqnarray} V=\sum_{i,j\le i}V_{ij}(\boldsymbol{r}_i-\boldsymbol{r}_j). \end{eqnarray} $$
In this case, a translation leads to \( \boldsymbol{r}_i\rightarrow \boldsymbol{r}_i+\boldsymbol{\epsilon} \), with the translation equally affecting the coordinates of each particle. Because the potential depends only on the relative coordinates, \( \delta V \) is manifestly zero. If one were to go through the exercise of calculating \( \delta V \) for small \( \boldsymbol{\epsilon} \), one would find that the term \( \nabla_i V(\boldsymbol{r}_i-\boldsymbol{r}_j) \) would be canceled by the term \( \nabla_jV(\boldsymbol{r}_i-\boldsymbol{r}_j) \).
The relation between symmetries of the potential and conserved quantities (also called constants of motion) is one of the most profound concepts one should gain from this course. It plays a critical role in all fields of physics. This is especially true in quantum mechanics, where a quantity \( A \) is conserved if its operator commutes with the Hamiltonian. For example if the momentum operator \( -i\hbar\partial_x \) commutes with the Hamiltonian, momentum is conserved, and clearly this operator commutes if the Hamiltonian (which represents the total energy, not just the potential) does not depend on \( x \). Also in quantum mechanics the angular momentum operator is \( L_z=-i\hbar\partial_\phi \). In fact, if the potential is unchanged by rotations about some axis, angular momentum about that axis is conserved. We return to this concept, from a more formal perspective, later in the course when Lagrangian mechanics is presented.
We will now venture into a study of a system which is energy conserving. The aim is to see if we (since it is not possible to solve the general equations analytically) we can develop stable numerical algorithms whose results we can trust!
We solve the equations of motion numerically. We will also compute quantities like the energy numerically.
We start with a simpler case first, the Earth-Sun system in two dimensions only. The gravitational force \( F_G \) on the earth from the sun is $$ \boldsymbol{F}_G=-\frac{GM_{\odot}M_E}{r^3}\boldsymbol{r}, $$ where \( G \) is the gravitational constant, $$ M_E=6\times 10^{24}\mathrm{Kg}, $$ the mass of Earth, $$ M_{\odot}=2\times 10^{30}\mathrm{Kg}, $$ the mass of the Sun and $$ r=1.5\times 10^{11}\mathrm{m}, $$ is the distance between Earth and the Sun. The latter defines what we call an astronomical unit AU. From Newton's second law we have then for the \( x \) direction $$ \frac{d^2x}{dt^2}=-\frac{F_{x}}{M_E}, $$ and $$ \frac{d^2y}{dt^2}=-\frac{F_{y}}{M_E}, $$ for the \( y \) direction.
Here we will use that \( x=r\cos{(\theta)} \), \( y=r\sin{(\theta)} \) and $$ r = \sqrt{x^2+y^2}. $$ We can rewrite $$ F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x, $$ and $$ F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y, $$ for the \( y \) direction.
We can rewrite these two equations $$ F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x, $$ and $$ F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y, $$ as four first-order coupled differential equations $$ \frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x, $$ $$ \frac{dx}{dt}=v_x, $$ $$ \frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y, $$ $$ \frac{dy}{dt}=v_y. $$
The four coupled differential equations $$ \frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x, $$ $$ \frac{dx}{dt}=v_x, $$ $$ \frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y, $$ $$ \frac{dy}{dt}=v_y, $$ can be turned into dimensionless equations or we can introduce astronomical units with \( 1 \) AU = \( 1.5\times 10^{11} \).
Using the equations from circular motion (with \( r =1\mathrm{AU} \)) $$ \frac{M_E v^2}{r} = F = \frac{GM_{\odot}M_E}{r^2}, $$ we have $$ GM_{\odot}=v^2r, $$ and using that the velocity of Earth (assuming circular motion) is \( v = 2\pi r/\mathrm{yr}=2\pi\mathrm{AU}/\mathrm{yr} \), we have $$ GM_{\odot}= v^2r = 4\pi^2 \frac{(\mathrm{AU})^3}{\mathrm{yr}^2}. $$
The four coupled differential equations can then be discretized using Euler's method as (with step length \( h \)) $$ v_{x,i+1}=v_{x,i}-h\frac{4\pi^2}{r_i^3}x_i, $$ $$ x_{i+1}=x_i+hv_{x,i}, $$ $$ v_{y,i+1}=v_{y,i}-h\frac{4\pi^2}{r_i^3}y_i, $$ $$ y_{i+1}=y_i+hv_{y,i}, $$
The code here implements Euler's method for the Earth-Sun system using a more compact way of representing the vectors. Alternatively, you could have spelled out all the variables \( v_x \), \( v_y \), \( x \) and \( y \) as one-dimensional arrays.
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
DeltaT = 0.001
#set up arrays
tfinal = 10 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Start integrating using Euler's method
for i in range(n-1):
# Set up the acceleration
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using Euler's forward method
v[i+1] = v[i] + DeltaT*a
r[i+1] = r[i] + DeltaT*v[i]
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
#ax.set_xlim(0, tfinal)
ax.set_ylabel('x[m]')
ax.set_xlabel('y[m]')
ax.plot(r[:,0], r[:,1])
fig.tight_layout()
save_fig("EarthSunEuler")
plt.show()
We notice here that Euler's method doesn't give a stable orbit. It means that we cannot trust Euler's method. In a deeper way, as we will see in homework 5, Euler's method does not conserve energy. It is an example of an integrator which is not symplectic.
Here we present thus two methods, which with simple changes allow us to avoid these pitfalls. The simplest possible extension is the so-called Euler-Cromer method. The changes we need to make to our code are indeed marginal here. We need simply to replace
r[i+1] = r[i] + DeltaT*v[i]
in the above code with the velocity at the new time \( t_{i+1} \)
r[i+1] = r[i] + DeltaT*v[i+1]
By this simple caveat we get stable orbits. Below we derive the Euler-Cromer method as well as one of the most utlized algorithms for sovling the above type of problems, the so-called Velocity-Verlet method.
Let us repeat Euler's method. We have a differential equation $$ \begin{equation} y'(t_i)=f(t_i,y_i) \label{_auto34} \end{equation} $$ and if we truncate at the first derivative, we have from the Taylor expansion $$ \begin{equation} y_{i+1}=y(t_i) + (\Delta t) f(t_i,y_i) + O(\Delta t^2), \label{eq:euler} \end{equation} $$ which when complemented with \( t_{i+1}=t_i+\Delta t \) forms the algorithm for the well-known Euler method. Note that at every step we make an approximation error of the order of \( O(\Delta t^2) \), however the total error is the sum over all steps \( N=(b-a)/(\Delta t) \) for \( t\in [a,b] \), yielding thus a global error which goes like \( NO(\Delta t^2)\approx O(\Delta t) \).
To make Euler's method more precise we can obviously decrease \( \Delta t \) (increase \( N \)), but this can lead to loss of numerical precision. Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like.
Euler's method is asymmetric in time, since it uses information about the derivative at the beginning of the time interval. This means that we evaluate the position at \( y_1 \) using the velocity at \( v_0 \). A simple variation is to determine \( x_{n+1} \) using the velocity at \( v_{n+1} \), that is (in a slightly more generalized form) $$ \begin{equation} y_{n+1}=y_{n}+ v_{n+1}+O(\Delta t^2) \label{_auto35} \end{equation} $$ and $$ \begin{equation} v_{n+1}=v_{n}+(\Delta t) a_{n}+O(\Delta t^2). \label{_auto36} \end{equation} $$ The acceleration \( a_n \) is a function of \( a_n(y_n, v_n, t_n) \) and needs to be evaluated as well. This is the Euler-Cromer method.
Exercise: go back to the above code with Euler's method and add the Euler-Cromer method.
Let us stay with \( x \) (position) and \( v \) (velocity) as the quantities we are interested in.
We have the Taylor expansion for the position given by $$ x_{i+1} = x_i+(\Delta t)v_i+\frac{(\Delta t)^2}{2}a_i+O((\Delta t)^3). $$ The corresponding expansion for the velocity is $$ v_{i+1} = v_i+(\Delta t)a_i+\frac{(\Delta t)^2}{2}v^{(2)}_i+O((\Delta t)^3). $$ Via Newton's second law we have normally an analytical expression for the derivative of the velocity, namely $$ a_i= \frac{d^2 x}{dt^2}\vert_{i}=\frac{d v}{dt}\vert_{i}= \frac{F(x_i,v_i,t_i)}{m}. $$
If we add to this the corresponding expansion for the derivative of the velocity $$ v^{(1)}_{i+1} = a_{i+1}= a_i+(\Delta t)v^{(2)}_i+O((\Delta t)^2)=a_i+(\Delta t)v^{(2)}_i+O((\Delta t)^2), $$ and retain only terms up to the second derivative of the velocity since our error goes as \( O(h^3) \), we have $$ (\Delta t)v^{(2)}_i\approx a_{i+1}-a_i. $$ We can then rewrite the Taylor expansion for the velocity as $$ v_{i+1} = v_i+\frac{(\Delta t)}{2}\left( a_{i+1}+a_{i}\right)+O((\Delta t)^3). $$
Our final equations for the position and the velocity become then $$ x_{i+1} = x_i+(\Delta t)v_i+\frac{(\Delta t)^2}{2}a_{i}+O((\Delta t)^3), $$ and $$ v_{i+1} = v_i+\frac{(\Delta t)}{2}\left(a_{i+1}+a_{i}\right)+O((\Delta t)^3). $$ Note well that the term \( a_{i+1} \) depends on the position at \( x_{i+1} \). This means that you need to calculate the position at the updated time \( t_{i+1} \) before the computing the next velocity. Note also that the derivative of the velocity at the time \( t_i \) used in the updating of the position can be reused in the calculation of the velocity update as well.
We can now easily add the Verlet method to our original code as
DeltaT = 0.01
#set up arrays
tfinal = 10
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up forces, air resistance FD, note now that we need the norm of the vecto
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
rabs = sqrt(sum(r[i+1]*r[i+1]))
anew = -4*(pi**2)*r[i+1]/(rabs**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
ax.set_ylabel('x[m]')
ax.set_xlabel('y[m]')
ax.plot(r[:,0], r[:,1])
fig.tight_layout()
save_fig("EarthSunVV")
plt.show()
You can easily generalize the calculation of the forces by defining a function which takes in as input the various variables. We leave this as a challenge to you.
In order to study the conservation of energy, we will need to perform a numerical integration, unless we can integrate analytically. Here we present the Trapezoidal rule as a the simplest possible approximation.
It is also useful to consider methods to integrate numerically. Let us consider the following case. We have classical electron which moves in the \( x \)-direction along a surface. The force from the surface is $$ \boldsymbol{F}(x)=-F_0\sin{(\frac{2\pi x}{b})}\boldsymbol{e}_x. $$ The constant \( b \) represents the distance between atoms at the surface of the material, \( F_0 \) is a constant and \( x \) is the position of the electron. Using the work-energy theorem we can find the work \( W \) done when moving an electron from a position \( x_0 \) to a final position \( x \) through the integral $$ W=-\int_{x_0}^x \boldsymbol{F}(x')dx' = \int_{x_0}^x F_0\sin{(\frac{2\pi x'}{b})} dx', $$ which results in $$ W=\frac{F_0b}{2\pi}\left[\cos{(\frac{2\pi x}{b})}-\cos{(\frac{2\pi x_0}{b})}\right]. $$
There are several numerical algorithms for finding an integral numerically. The more familiar ones like the rectangular rule or the trapezoidal rule have simple geometric interpretations.
Let us look at the mathematical details of what are called equal-step methods, also known as Newton-Cotes quadrature.
The integral $$ \begin{equation} I=\int_a^bf(x) dx \label{eq:integraldef} \end{equation} $$ has a very simple meaning. The integral is the area enscribed by the function \( f(x) \) starting from \( x=a \) to \( x=b \). It is subdivided in several smaller areas whose evaluation is to be approximated by different techniques. The areas under the curve can for example be approximated by rectangular boxes or trapezoids.
In considering equal step methods, our basic approach is that of approximating a function \( f(x) \) with a polynomial of at most degree \( N-1 \), given \( N \) integration points. If our polynomial is of degree \( 1 \), the function will be approximated with \( f(x)\approx a_0+a_1x \).
The algorithm for these integration methods is rather simple, and the number of approximations perhaps unlimited!
A very simple approach is the so-called midpoint or rectangle method. In this case the integration area is split in a given number of rectangles with length \( h \) and height given by the mid-point value of the function. This gives the following simple rule for approximating an integral $$ \begin{equation} I=\int_a^bf(x) dx \approx h\sum_{i=1}^N f(x_{i-1/2}), \label{eq:rectangle} \end{equation} $$ where \( f(x_{i-1/2}) \) is the midpoint value of \( f \) for a given rectangle. We will discuss its truncation error below. It is easy to implement this algorithm, as shown below
The correct mathematical expression for the local error for the rectangular rule \( R_i(h) \) for element \( i \) is $$ \begin{equation*} \int_{-h}^hf(x)dx - R_i(h)=-\frac{h^3}{24}f^{(2)}(\xi), \end{equation*} $$ and the global error reads $$ \begin{equation*} \int_a^bf(x)dx -R_h(f)=-\frac{b-a}{24}h^2f^{(2)}(\xi), \end{equation*} $$ where \( R_h \) is the result obtained with rectangular rule and \( \xi \in [a,b] \).
We go back to our simple example above and set \( F_0=b=1 \) and choose \( x_0=0 \) and \( x=1/2 \), and have $$ W=\frac{1}{\pi}. $$ The code here computes the integral using the rectangle rule and \( n=100 \) integration points we have a relative error of \( 10^{-5} \).
from math import sin, pi
import numpy as np
from sympy import Symbol, integrate
# function for the Rectangular rule
def Rectangular(a,b,f,n):
h = (b-a)/float(n)
s = 0
for i in range(0,n,1):
x = (i+0.5)*h
s = s+ f(x)
return h*s
# function to integrate
def function(x):
return sin(2*pi*x)
# define integration limits and integration points
a = 0.0; b = 0.5;
n = 100
Exact = 1./pi
print("Relative error= ", abs( (Rectangular(a,b,function,n)-Exact)/Exact))
The other integral gives $$ \begin{equation*} \int_{x_0-h}^{x_0}f(x)dx=\frac{h}{2}\left(f(x_0) + f(x_0-h)\right)+O(h^3), \end{equation*} $$ and adding up we obtain $$ \begin{equation} \int_{x_0-h}^{x_0+h}f(x)dx=\frac{h}{2}\left(f(x_0+h) + 2f(x_0) + f(x_0-h)\right)+O(h^3), \label{eq:trapez} \end{equation} $$ which is the well-known trapezoidal rule. Concerning the error in the approximation made, \( O(h^3)=O((b-a)^3/N^3) \), you should note that this is the local error. Since we are splitting the integral from \( a \) to \( b \) in \( N \) pieces, we will have to perform approximately \( N \) such operations.
This means that the global error goes like \( \approx O(h^2) \). The trapezoidal reads then $$ \begin{equation} I=\int_a^bf(x) dx=h\left(f(a)/2 + f(a+h) +f(a+2h)+ \dots +f(b-h)+ f_{b}/2\right), \label{eq:trapez1} \end{equation} $$ with a global error which goes like \( O(h^2) \).
Hereafter we use the shorthand notations \( f_{-h}=f(x_0-h) \), \( f_{0}=f(x_0) \) and \( f_{h}=f(x_0+h) \).
The correct mathematical expression for the local error for the trapezoidal rule is $$ \begin{equation*} \int_a^bf(x)dx -\frac{b-a}{2}\left[f(a)+f(b)\right]=-\frac{h^3}{12}f^{(2)}(\xi), \end{equation*} $$ and the global error reads $$ \begin{equation*} \int_a^bf(x)dx -T_h(f)=-\frac{b-a}{12}h^2f^{(2)}(\xi), \end{equation*} $$ where \( T_h \) is the trapezoidal result and \( \xi \in [a,b] \).
The trapezoidal rule is easy to implement numerically through the following simple algorithm
We use the same function and integrate now using the trapoezoidal rule.
import numpy as np
from sympy import Symbol, integrate
# function for the trapezoidal rule
def Trapez(a,b,f,n):
h = (b-a)/float(n)
s = 0
x = a
for i in range(1,n,1):
x = x+h
s = s+ f(x)
s = 0.5*(f(a)+f(b)) +s
return h*s
# function to integrate
def function(x):
return sin(2*pi*x)
# define integration limits and integration points
a = 0.0; b = 0.5;
n = 100
Exact = 1./pi
print("Relative error= ", abs( (Trapez(a,b,function,n)-Exact)/Exact))
Instead of using the above first-order polynomials approximations for \( f \), we attempt at using a second-order polynomials. In this case we need three points in order to define a second-order polynomial approximation $$ \begin{equation*} f(x) \approx P_2(x)=a_0+a_1x+a_2x^2. \end{equation*} $$ Using again Lagrange's interpolation formula we have $$ \begin{equation*} P_2(x)=\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_2+ \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_1+ \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0. \end{equation*} $$ Inserting this formula in the integral of Eq. \eqref{eq:hhint} we obtain $$ \begin{equation*} \int_{-h}^{+h}f(x)dx=\frac{h}{3}\left(f_h + 4f_0 + f_{-h}\right)+O(h^5), \end{equation*} $$ which is Simpson's rule.
Note that the improved accuracy in the evaluation of the derivatives gives a better error approximation, \( O(h^5) \) vs.\ \( O(h^3) \) . But this is again the local error approximation. Using Simpson's rule we can easily compute the integral of Eq. \eqref{eq:integraldef} to be $$ \begin{equation} I=\int_a^bf(x) dx=\frac{h}{3}\left(f(a) + 4f(a+h) +2f(a+2h)+ \dots +4f(b-h)+ f_{b}\right), \label{eq:simpson} \end{equation} $$ with a global error which goes like \( O(h^4) \).
More formal expressions for the local and global errors are for the local error $$ \begin{equation*} \int_a^bf(x)dx -\frac{b-a}{6}\left[f(a)+4f((a+b)/2)+f(b)\right]=-\frac{h^5}{90}f^{(4)}(\xi), \end{equation*} $$ and for the global error $$ \begin{equation*} \int_a^bf(x)dx -S_h(f)=-\frac{b-a}{180}h^4f^{(4)}(\xi). \end{equation*} $$ with \( \xi\in[a,b] \) and \( S_h \) the results obtained with Simpson's method.
from math import sin, pi
import numpy as np
from sympy import Symbol, integrate
# function for the trapezoidal rule
def Simpson(a,b,f,n):
h = (b-a)/float(n)
sum = f(a)/float(2);
for i in range(1,n):
sum = sum + f(a+i*h)*(3+(-1)**(i+1))
sum = sum + f(b)/float(2)
return sum*h/3.0
# function to integrate
def function(x):
return sin(2*pi*x)
# define integration limits and integration points
a = 0.0; b = 0.5;
n = 100
Exact = 1./pi
print("Relative error= ", abs( (Simpson(a,b,function,n)-Exact)/Exact))
We see that Simpson's rule gives a much better estimation of the relative error with the same amount of points as we had for the Rectangle rule and the Trapezoidal rule.
The harmonic oscillator is omnipresent in physics. Although you may think of this as being related to springs, it, or an equivalent mathematical representation, appears in just about any problem where a mode is sitting near its potential energy minimum. At that point, \( \partial_x V(x)=0 \), and the first non-zero term (aside from a constant) in the potential energy is that of a harmonic oscillator. In a solid, sound modes (phonons) are built on a picture of coupled harmonic oscillators, and in relativistic field theory the fundamental interactions are also built on coupled oscillators positioned infinitesimally close to one another in space. The phenomena of a resonance of an oscillator driven at a fixed frequency plays out repeatedly in atomic, nuclear and high-energy physics, when quantum mechanically the evolution of a state oscillates according to \( e^{-iEt} \) and exciting discrete quantum states has very similar mathematics as exciting discrete states of an oscillator.
The potential energy for a single particle as a function of its position \( x \) can be written as a Taylor expansion about some point \( x_0 \) $$ \begin{equation} V(x)=V(x_0)+(x-x_0)\left.\partial_xV(x)\right|_{x_0}+\frac{1}{2}(x-x_0)^2\left.\partial_x^2V(x)\right|_{x_0} +\frac{1}{3!}\left.\partial_x^3V(x)\right|_{x_0}+\cdots \label{_auto37} \end{equation} $$ If the position \( x_0 \) is at the minimum of the resonance, the first two non-zero terms of the potential are $$ \begin{eqnarray} V(x)&\approx& V(x_0)+\frac{1}{2}(x-x_0)^2\left.\partial_x^2V(x)\right|_{x_0},\\ \nonumber &=&V(x_0)+\frac{1}{2}k(x-x_0)^2,~~~~k\equiv \left.\partial_x^2V(x)\right|_{x_0},\\ \nonumber F&=&-\partial_xV(x)=-k(x-x_0). \end{eqnarray} $$ Put into Newton's 2nd law (assuming \( x_0=0 \)), $$ \begin{eqnarray} m\ddot{x}&=&-kx,\\ x&=&A\cos(\omega_0 t-\phi),~~~\omega_0=\sqrt{k/m}. \end{eqnarray} $$
Here \( A \) and \( \phi \) are arbitrary. Equivalently, one could have written this as \( A\cos(\omega_0 t)+B\sin(\omega_0 t) \), or as the real part of \( Ae^{i\omega_0 t} \). In this last case \( A \) could be an arbitrary complex constant. Thus, there are 2 arbitrary constants (either \( A \) and \( B \) or \( A \) and \( \phi \), or the real and imaginary part of one complex constant. This is the expectation for a second order differential equation, and also agrees with the physical expectation that if you know a particle's initial velocity and position you should be able to define its future motion, and that those two arbitrary conditions should translate to two arbitrary constants.
A key feature of harmonic motion is that the system repeats itself after a time \( T=1/f \), where \( f \) is the frequency, and \( \omega=2\pi f \) is the angular frequency. The period of the motion is independent of the amplitude. However, this independence is only exact when one can neglect higher terms of the potential, \( x^3, x^4\cdots \). Once can neglect these terms for sufficiently small amplitudes, and for larger amplitudes the motion is no longer purely sinusoidal, and even though the motion repeats itself, the time for repeating the motion is no longer independent of the amplitude.
One can also calculate the velocity and the kinetic energy as a function of time, $$ \begin{eqnarray} \dot{x}&=&-\omega_0A\sin(\omega_0 t-\phi),\\ \nonumber K&=&\frac{1}{2}m\dot{x}^2=\frac{m\omega_0^2A^2}{2}\sin^2(\omega_0t-\phi),\\ \nonumber &=&\frac{k}{2}A^2\sin^2(\omega_0t-\phi). \end{eqnarray} $$
The total energy is then $$ \begin{equation} E=K+V=\frac{1}{2}m\dot{x}^2+\frac{1}{2}kx^2=\frac{1}{2}kA^2. \label{_auto38} \end{equation} $$
The total energy then goes as the square of the amplitude.
A pendulum is an example of a harmonic oscillator. By expanding the kinetic and potential energies for small angles find the frequency for a pendulum of length \( L \) with all the mass \( m \) centered at the end by writing the eq.s of motion in the form of a harmonic oscillator.
The potential energy and kinetic energies are (for \( x \) being the displacement) $$ \begin{eqnarray*} V&=&mgL(1-\cos\theta)\approx mgL\frac{x^2}{2L^2},\\ K&=&\frac{1}{2}mL^2\dot{\theta}^2\approx \frac{m}{2}\dot{x}^2. \end{eqnarray*} $$
For small \( x \) Newton's 2nd law becomes $$ m\ddot{x}=-\frac{mg}{L}x, $$
and the spring constant would appear to be \( k=mg/L \), which makes the frequency equal to \( \omega_0=\sqrt{g/L} \). Note that the frequency is independent of the mass.
We consider only the case where the damping force is proportional to the velocity. This is counter to dragging friction, where the force is proportional in strength to the normal force and independent of velocity, and is also inconsistent with wind resistance, where the magnitude of the drag force is proportional the square of the velocity. Rolling resistance does seem to be mainly proportional to the velocity. However, the main motivation for considering damping forces proportional to the velocity is that the math is more friendly. This is because the differential equation is linear, i.e. each term is of order \( x \), \( \dot{x} \), \( \ddot{x}\cdots \), or even terms with no mention of \( x \), and there are no terms such as \( x^2 \) or \( x\ddot{x} \). The equations of motion for a spring with damping force \( -b\dot{x} \) are $$ \begin{equation} m\ddot{x}+b\dot{x}+kx=0. \label{_auto39} \end{equation} $$
Just to make the solution a bit less messy, we rewrite this equation as $$ \begin{equation} \label{eq:dampeddiffyq} \ddot{x}+2\beta\dot{x}+\omega_0^2x=0,~~~~\beta\equiv b/2m,~\omega_0\equiv\sqrt{k/m}. \end{equation} $$
Both \( \beta \) and \( \omega \) have dimensions of inverse time. To find solutions (see appendix C in the text) you must make an educated guess at the form of the solution. To do this, first realize that the solution will need an arbitrary normalization \( A \) because the equation is linear. Secondly, realize that if the form is $$ \begin{equation} x=Ae^{rt} \label{_auto40} \end{equation} $$
that each derivative simply brings out an extra power of \( r \). This means that the \( Ae^{rt} \) factors out and one can simply solve for an equation for \( r \). Plugging this form into Eq. \eqref{eq:dampeddiffyq}, $$ \begin{equation} r^2+2\beta r+\omega_0^2=0. \label{_auto41} \end{equation} $$
Because this is a quadratic equation there will be two solutions, $$ \begin{equation} r=-\beta\pm\sqrt{\beta^2-\omega_0^2}. \label{_auto42} \end{equation} $$
We refer to the two solutions as \( r_1 \) and \( r_2 \) corresponding to the \( + \) and \( - \) roots. As expected, there should be two arbitrary constants involved in the solution, $$ \begin{equation} x=A_1e^{r_1t}+A_2e^{r_2t}, \label{_auto43} \end{equation} $$
where the coefficients \( A_1 \) and \( A_2 \) are determined by initial conditions.
The roots listed above, \( \sqrt{\omega_0^2-\beta_0^2} \), will be imaginary if the damping is small and \( \beta < \omega_0 \). In that case, \( r \) is complex and the factor \( e{rt} \) will have some oscillatory behavior. If the roots are real, there will only be exponentially decaying solutions. There are three cases:
Here we have made use of the identity \( e^{i\omega't}=\cos\omega't+i\sin\omega't \). Because the constants are arbitrary, and because the real and imaginary parts are both solutions individually, we can simply consider the real part of the solution alone: $$ \begin{eqnarray} \label{eq:homogsolution} x&=&B_1e^{-\beta t}\cos\omega't+B_2e^{-\beta t}\sin\omega't,\\ \nonumber \omega'&\equiv&\sqrt{\omega_0^2-\beta^2}. \end{eqnarray} $$
In this case the two terms involving \( r_1 \) and \( r_2 \) are identical because \( \omega'=0 \). Because we need to arbitrary constants, there needs to be another solution. This is found by simply guessing, or by taking the limit of \( \omega'\rightarrow 0 \) from the underdamped solution. The solution is then $$ \begin{equation} \label{eq:criticallydamped} x=Ae^{-\beta t}+Bte^{-\beta t}. \end{equation} $$
The critically damped solution is interesting because the solution approaches zero quickly, but does not oscillate. For a problem with zero initial velocity, the solution never crosses zero. This is a good choice for designing shock absorbers or swinging doors.
This solution will also never pass the origin more than once, and then only if the initial velocity is strong and initially toward zero.
Given \( b \), \( m \) and \( \omega_0 \), find \( x(t) \) for a particle whose initial position is \( x=0 \) and has initial velocity \( v_0 \) (assuming an underdamped solution).
The solution is of the form, $$ \begin{eqnarray*} x&=&e^{-\beta t}\left[A_1\cos(\omega' t)+A_2\sin\omega't\right],\\ \dot{x}&=&-\beta x+\omega'e^{-\beta t}\left[-A_1\sin\omega't+A_2\cos\omega't\right].\\ \omega'&\equiv&\sqrt{\omega_0^2-\beta^2},~~~\beta\equiv b/2m. \end{eqnarray*} $$
From the initial conditions, \( A_1=0 \) because \( x(0)=0 \) and \( \omega'A_2=v_0 \). So $$ x=\frac{v_0}{\omega'}e^{-\beta t}\sin\omega't. $$
Here we study first the case without additional friction term and scale our equation in terms of a dimensionless time \( \tau \).
Let us remind ourselves about the differential equation we want to solve (the general case with damping due to friction) $$ m\frac{d^2x}{dt^2} + b\frac{dx}{dt}+kx(t) =0. $$ We divide by \( m \) and introduce \( \omega_0^2=\sqrt{k/m} \) and obtain $$ \frac{d^2x}{dt^2} + \frac{b}{m}\frac{dx}{dt}+\omega_0^2x(t) =0. $$
Thereafter we introduce a dimensionless time \( \tau = t\omega_0 \) (check that the dimensionality is correct) and rewrite our equation as $$ \frac{d^2x}{d\tau^2} + \frac{b}{m\omega_0}\frac{dx}{d\tau}+x(\tau) =0, $$ which gives us $$ \frac{d^2x}{d\tau^2} + \frac{b}{m\omega_0}\frac{dx}{d\tau}+x(\tau) =0. $$
We then define \( \gamma = b/(2m\omega_0) \) and rewrite our equations as $$ \frac{d^2x}{d\tau^2} + 2\gamma\frac{dx}{d\tau}+x(\tau) =0. $$ This is the equation we will code below. The first version employs the Euler-Cromer method.
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
from pylab import plt, mpl
plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
DeltaT = 0.001
#set up arrays
tfinal = 20 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, v, and x
t = np.zeros(n)
v = np.zeros(n)
x = np.zeros(n)
# Initial conditions as simple one-dimensional arrays of time
x0 = 1.0
v0 = 0.0
x[0] = x0
v[0] = v0
gamma = 0.0
# Start integrating using Euler-Cromer's method
for i in range(n-1):
# Set up the acceleration
# Here you could have defined your own function for this
a = -2*gamma*v[i]-x[i]
# update velocity, time and position
v[i+1] = v[i] + DeltaT*a
x[i+1] = x[i] + DeltaT*v[i+1]
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
#ax.set_xlim(0, tfinal)
ax.set_ylabel('x[m]')
ax.set_xlabel('t[s]')
ax.plot(t, x)
fig.tight_layout()
save_fig("BlockEulerCromer")
plt.show()
When setting up the value of \( \gamma \) we see that for \( \gamma=0 \) we get the simple oscillatory motion with no damping. Choosing \( \gamma < 1 \) leads to the classical underdamped case with oscillatory motion, but where the motion comes to an end.
Choosing \( \gamma =1 \) leads to what normally is called critical damping and \( \gamma> 1 \) leads to critical overdamping. Try it out and try also to change the initial position and velocity. Setting \( \gamma=1 \) yields a situation, as discussed above, where the solution approaches quickly zero and does not oscillate. With zero initial velocity it will never cross zero.
Here, we consider the force $$ \begin{equation} F=-kx-b\dot{x}+F_0\cos\omega t, \label{_auto44} \end{equation} $$
which leads to the differential equation $$ \begin{equation} \label{eq:drivenosc} \ddot{x}+2\beta\dot{x}+\omega_0^2x=(F_0/m)\cos\omega t. \end{equation} $$
Consider a single solution with no arbitrary constants, which we will call a {\it particular solution}, \( x_p(t) \). It should be emphasized that this is {\bf A} particular solution, because there exists an infinite number of such solutions because the general solution should have two arbitrary constants. Now consider solutions to the same equation without the driving term, which include two arbitrary constants. These are called either {\it homogenous solutions} or {\it complementary solutions}, and were given in the previous section, e.g. Eq. \eqref{eq:homogsolution} for the underdamped case. The homogenous solution already incorporates the two arbitrary constants, so any sum of a homogenous solution and a particular solution will represent the {\it general solution} of the equation. The general solution incorporates the two arbitrary constants \( A \) and \( B \) to accommodate the two initial conditions. One could have picked a different particular solution, i.e. the original particular solution plus any homogenous solution with the arbitrary constants \( A_p \) and \( B_p \) chosen at will. When one adds in the homogenous solution, which has adjustable constants with arbitrary constants \( A' \) and \( B' \), to the new particular solution, one can get the same general solution by simply adjusting the new constants such that \( A'+A_p=A \) and \( B'+B_p=B \). Thus, the choice of \( A_p \) and \( B_p \) are irrelevant, and when choosing the particular solution it is best to make the simplest choice possible.
To find a particular solution, one first guesses at the form, $$ \begin{equation} \label{eq:partform} x_p(t)=D\cos(\omega t-\delta), \end{equation} $$
and rewrite the differential equation as $$ \begin{equation} D\left\{-\omega^2\cos(\omega t-\delta)-2\beta\omega\sin(\omega t-\delta)+\omega_0^2\cos(\omega t-\delta)\right\}=\frac{F_0}{m}\cos(\omega t). \label{_auto45} \end{equation} $$
One can now use angle addition formulas to get $$ \begin{eqnarray} D\left\{(-\omega^2\cos\delta+2\beta\omega\sin\delta+\omega_0^2\cos\delta)\cos(\omega t)\right.&&\\ \nonumber \left.+(-\omega^2\sin\delta-2\beta\omega\cos\delta+\omega_0^2\sin\delta)\sin(\omega t)\right\} &=&\frac{F_0}{m}\cos(\omega t). \end{eqnarray} $$
Both the \( \cos \) and \( \sin \) terms need to equate if the expression is to hold at all times. Thus, this becomes two equations $$ \begin{eqnarray} D\left\{-\omega^2\cos\delta+2\beta\omega\sin\delta+\omega_0^2\cos\delta\right\}&=&\frac{F_0}{m}\\ \nonumber -\omega^2\sin\delta-2\beta\omega\cos\delta+\omega_0^2\sin\delta&=&0. \end{eqnarray} $$
After dividing by \( \cos\delta \), the lower expression leads to $$ \begin{equation} \tan\delta=\frac{2\beta\omega}{\omega_0^2-\omega^2}. \label{_auto46} \end{equation} $$
Using the identities \( \tan^2+1=\csc^2 \) and \( \sin^2+\cos^2=1 \), one can also express \( \sin\delta \) and \( \cos\delta \), $$ \begin{eqnarray} \sin\delta&=&\frac{2\beta\omega}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}},\\ \nonumber \cos\delta&=&\frac{(\omega_0^2-\omega^2)}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}} \end{eqnarray} $$
Inserting the expressions for \( \cos\delta \) and \( \sin\delta \) into the expression for \( D \), $$ \begin{equation} \label{eq:Ddrive} D=\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}}. \end{equation} $$
For a given initial condition, e.g. initial displacement and velocity, one must add the homogenous solution then solve for the two arbitrary constants. However, because the homogenous solutions decay with time as \( e^{-\beta t} \), the particular solution is all that remains at large times, and is therefore the steady state solution. Because the arbitrary constants are all in the homogenous solution, all memory of the initial conditions are lost at large times, \( t>>1/\beta \).
The amplitude of the motion, \( D \), is linearly proportional to the driving force (\( F_0/m \)), but also depends on the driving frequency \( \omega \). For small \( \beta \) the maximum will occur at \( \omega=\omega_0 \). This is referred to as a resonance. In the limit \( \beta\rightarrow 0 \) the amplitude at resonance approaches infinity.
Here, we derive the same expressions as in Equations \eqref{eq:partform} and \eqref{eq:Ddrive} but express the driving forces as $$ \begin{eqnarray} F(t)&=&F_0e^{i\omega t}, \end{eqnarray} $$ rather than as \( F_0\cos\omega t \). The real part of \( F \) is the same as before. For the differential equation, $$ \begin{eqnarray} \label{eq:compdrive} \ddot{x}+2\beta\dot{x}+\omega_0^2x&=&\frac{F_0}{m}e^{i\omega t}, \end{eqnarray} $$
one can treat \( x(t) \) as an imaginary function. Because the operations \( d^2/dt^2 \) and \( d/dt \) are real and thus do not mix the real and imaginary parts of \( x(t) \), Eq. \eqref{eq:compdrive} is effectively 2 equations. Because \( e^{\omega t}=\cos\omega t+i\sin\omega t \), the real part of the solution for \( x(t) \) gives the solution for a driving force \( F_0\cos\omega t \), and the imaginary part of \( x \) corresponds to the case where the driving force is \( F_0\sin\omega t \). It is rather easy to solve for the complex \( x \) in this case, and by taking the real part of the solution, one finds the answer for the \( \cos\omega t \) driving force.
We assume a simple form for the particular solution $$ \begin{equation} x_p=De^{i\omega t}, \label{_auto47} \end{equation} $$
where \( D \) is a complex constant.
From Eq. \eqref{eq:compdrive} one inserts the form for \( x_p \) above to get $$ \begin{eqnarray} D\left\{-\omega^2+2i\beta\omega+\omega_0^2\right\}e^{i\omega t}=(F_0/m)e^{i\omega t},\\ \nonumber D=\frac{F_0/m}{(\omega_0^2-\omega^2)+2i\beta\omega}. \end{eqnarray} $$
The norm and phase for \( D=|D|e^{-i\delta} \) can be read by inspection, $$ \begin{equation} |D|=\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}},~~~~\tan\delta=\frac{2\beta\omega}{\omega_0^2-\omega^2}. \label{_auto48} \end{equation} $$
This is the same expression for \( \delta \) as before. One then finds \( x_p(t) \), $$ \begin{eqnarray} \label{eq:fastdriven1} x_p(t)&=&\Re\frac{(F_0/m)e^{i\omega t-i\delta}}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}\\ \nonumber &=&\frac{(F_0/m)\cos(\omega t-\delta)}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}. \end{eqnarray} $$ This is the same answer as before. If one wished to solve for the case where \( F(t)= F_0\sin\omega t \), the imaginary part of the solution would work $$ \begin{eqnarray} \label{eq:fastdriven2} x_p(t)&=&\Im\frac{(F_0/m)e^{i\omega t-i\delta}}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}\\ \nonumber &=&\frac{(F_0/m)\sin(\omega t-\delta)}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}. \end{eqnarray} $$
Consider the damped and driven harmonic oscillator worked out above. Given \( F_0, m,\beta \) and \( \omega_0 \), solve for the complete solution \( x(t) \) for the case where \( F=F_0\sin\omega t \) with initial conditions \( x(t=0)=0 \) and \( v(t=0)=0 \). Assume the underdamped case.
The general solution including the arbitrary constants includes both the homogenous and particular solutions, $$ \begin{eqnarray*} x(t)&=&\frac{F_0}{m}\frac{\sin(\omega t-\delta)}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}} +A\cos\omega't e^{-\beta t}+B\sin\omega't e^{-\beta t}. \end{eqnarray*} $$
The quantities \( \delta \) and \( \omega' \) are given earlier in the section, \( \omega'=\sqrt{\omega_0^2-\beta^2}, \delta=\tan^{-1}(2\beta\omega/(\omega_0^2-\omega^2) \). Here, solving the problem means finding the arbitrary constants \( A \) and \( B \). Satisfying the initial conditions for the initial position and velocity: $$ \begin{eqnarray*} x(t=0)=0&=&-\eta\sin\delta+A,\\ v(t=0)=0&=&\omega\eta\cos\delta-\beta A+\omega'B,\\ \eta&\equiv&\frac{F_0}{m}\frac{1}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}. \end{eqnarray*} $$
The problem is now reduced to 2 equations and 2 unknowns, \( A \) and \( B \). The solution is $$ \begin{eqnarray} A&=& \eta\sin\delta ,~~~B=\frac{-\omega\eta\cos\delta+\beta\eta\sin\delta}{\omega'}. \end{eqnarray} $$
From the previous two sections, the particular solution for a driving force, \( F=F_0\cos\omega t \), is $$ \begin{eqnarray} x_p(t)&=&\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}}\cos(\omega_t-\delta),\\ \nonumber \delta&=&\tan^{-1}\left(\frac{2\beta\omega}{\omega_0^2-\omega^2}\right). \end{eqnarray} $$
If one fixes the driving frequency \( \omega \) and adjusts the fundamental frequency \( \omega_0=\sqrt{k/m} \), the maximum amplitude occurs when \( \omega_0=\omega \) because that is when the term from the denominator \( (\omega_0^2-\omega^2)^2+4\omega^2\beta^2 \) is at a minimum. This is akin to dialing into a radio station. However, if one fixes \( \omega_0 \) and adjusts the driving frequency one minimize with respect to \( \omega \), e.g. set $$ \begin{equation} \frac{d}{d\omega}\left[(\omega_0^2-\omega^2)^2+4\omega^2\beta^2\right]=0, \label{_auto49} \end{equation} $$
and one finds that the maximum amplitude occurs when \( \omega=\sqrt{\omega_0^2-2\beta^2} \). If \( \beta \) is small relative to \( \omega_0 \), one can simply state that the maximum amplitude is $$ \begin{equation} x_{\rm max}\approx\frac{F_0}{2m\beta \omega_0}. \label{_auto50} \end{equation} $$ $$ \begin{eqnarray} \frac{4\omega^2\beta^2}{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}=\frac{1}{2}. \end{eqnarray} $$
For small damping this occurs when \( \omega=\omega_0\pm \beta \), so the \( FWHM\approx 2\beta \). For the purposes of tuning to a specific frequency, one wants the width to be as small as possible. The ratio of \( \omega_0 \) to \( FWHM \) is known as the {\it quality} factor, or \( Q \) factor, $$ \begin{equation} Q\equiv \frac{\omega_0}{2\beta}. \label{_auto51} \end{equation} $$
Solving the problem of driven oscillations numerically gives us much more flexibility to study different types of driving forces. We can reuse our earlier code by simply adding a driving force. If we stay in the \( x \)-direction only this can be easily done by adding a term \( F_{\mathrm{ext}}(x,t) \). Note that we have kept it rather general here, allowing for both a spatial and a temporal dependence.
Before we dive into the code, we need to briefly remind ourselves about the equations we started with for the case with damping, namely $$ m\frac{d^2x}{dt^2} + b\frac{dx}{dt}+kx(t) =0, $$
with no external force applied to the system.
Let us now for simplicty assume that our external force is given by $$ F_{\mathrm{ext}}(t) = F_0\cos{(\omega t)}, $$ where \( F_0 \) is a constant (what is its dimension?) and \( \omega \) is the frequency of the applied external driving force. Small question: would you expect energy to be conserved now?
Introducing the external force into our lovely differential equation and dividing by \( m \) and introducing \( \omega_0^2=\sqrt{k/m} \) we have $$ \frac{d^2x}{dt^2} + \frac{b}{m}\frac{dx}{dt}+\omega_0^2x(t) =\frac{F_0}{m}\cos{(\omega t)}, $$
Thereafter we introduce a dimensionless time \( \tau = t\omega_0 \) and a dimensionless frequency \( \tilde{\omega}=\omega/\omega_0 \). We have then $$ \frac{d^2x}{d\tau^2} + \frac{b}{m\omega_0}\frac{dx}{d\tau}+x(\tau) =\frac{F_0}{m\omega_0^2}\cos{(\tilde{\omega}\tau)}, $$ Introducing a new amplitude \( \tilde{F} =F_0/(m\omega_0^2) \) (check dimensionality again) we have $$ \frac{d^2x}{d\tau^2} + \frac{b}{m\omega_0}\frac{dx}{d\tau}+x(\tau) =\tilde{F}\cos{(\tilde{\omega}\tau)}. $$
Our final step, as we did in the case of various types of damping, is to define \( \gamma = b/(2m\omega_0) \) and rewrite our equations as $$ \frac{d^2x}{d\tau^2} + 2\gamma\frac{dx}{d\tau}+x(\tau) =\tilde{F}\cos{(\tilde{\omega}\tau)}. $$ This is the equation we will code below using the Euler-Cromer method.
DeltaT = 0.001
#set up arrays
tfinal = 20 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, v, and x
t = np.zeros(n)
v = np.zeros(n)
x = np.zeros(n)
# Initial conditions as one-dimensional arrays of time
x0 = 1.0
v0 = 0.0
x[0] = x0
v[0] = v0
gamma = 0.2
Omegatilde = 0.5
Ftilde = 1.0
# Start integrating using Euler-Cromer's method
for i in range(n-1):
# Set up the acceleration
# Here you could have defined your own function for this
a = -2*gamma*v[i]-x[i]+Ftilde*cos(t[i]*Omegatilde)
# update velocity, time and position
v[i+1] = v[i] + DeltaT*a
x[i+1] = x[i] + DeltaT*v[i+1]
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
ax.set_ylabel('x[m]')
ax.set_xlabel('t[s]')
ax.plot(t, x)
fig.tight_layout()
save_fig("ForcedBlockEulerCromer")
plt.show()
In the above example we have focused on the Euler-Cromer method. This method has a local truncation error which is proportional to \( \Delta t^2 \) and thereby a global error which is proportional to \( \Delta t \). We can improve this by using the Runge-Kutta family of methods. The widely popular Runge-Kutta to fourth order or just RK4 has indeed a much better truncation error. The RK4 method has a global error which is proportional to \( \Delta t \).
Let us revisit this method and see how we can implement it for the above example.
Runge-Kutta (RK) methods are based on Taylor expansion formulae, but yield in general better algorithms for solutions of an ordinary differential equation. The basic philosophy is that it provides an intermediate step in the computation of \( y_{i+1} \).
To see this, consider first the following definitions $$ \begin{equation} \frac{dy}{dt}=f(t,y), \label{_auto52} \end{equation} $$ and $$ \begin{equation} y(t)=\int f(t,y) dt, \label{_auto53} \end{equation} $$ and $$ \begin{equation} y_{i+1}=y_i+ \int_{t_i}^{t_{i+1}} f(t,y) dt. \label{_auto54} \end{equation} $$
To demonstrate the philosophy behind RK methods, let us consider the second-order RK method, RK2. The first approximation consists in Taylor expanding \( f(t,y) \) around the center of the integration interval \( t_i \) to \( t_{i+1} \), that is, at \( t_i+h/2 \), \( h \) being the step. Using the midpoint formula for an integral, defining \( y(t_i+h/2) = y_{i+1/2} \) and \( t_i+h/2 = t_{i+1/2} \), we obtain $$ \begin{equation} \int_{t_i}^{t_{i+1}} f(t,y) dt \approx hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \label{_auto55} \end{equation} $$ This means in turn that we have $$ \begin{equation} y_{i+1}=y_i + hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \label{_auto56} \end{equation} $$
However, we do not know the value of \( y_{i+1/2} \). Here comes thus the next approximation, namely, we use Euler's method to approximate \( y_{i+1/2} \). We have then $$ \begin{equation} y_{(i+1/2)}=y_i + \frac{h}{2}\frac{dy}{dt}=y(t_i) + \frac{h}{2}f(t_i,y_i). \label{_auto57} \end{equation} $$ This means that we can define the following algorithm for the second-order Runge-Kutta method, RK2. $$ \begin{equation} k_1=hf(t_i,y_i), \label{_auto58} \end{equation} $$ $$ \begin{equation} k_2=hf(t_{i+1/2},y_i+k_1/2), \label{_auto59} \end{equation} $$ with the final value $$ \begin{equation} y_{i+i}\approx y_i + k_2 +O(h^3). \label{_auto60} \end{equation} $$
The difference between the previous one-step methods is that we now need an intermediate step in our evaluation, namely \( t_i+h/2 = t_{(i+1/2)} \) where we evaluate the derivative \( f \). This involves more operations, but the gain is a better stability in the solution.
The fourth-order Runge-Kutta, RK4, has the following algorithm $$ k_1=hf(t_i,y_i) \hspace{0.5cm} k_2=hf(t_i+h/2,y_i+k_1/2) $$ $$ k_3=hf(t_i+h/2,y_i+k_2/2)\hspace{0.5cm} k_4=hf(t_i+h,y_i+k_3) $$ with the final result $$ y_{i+1}=y_i +\frac{1}{6}\left( k_1 +2k_2+2k_3+k_4\right). $$ Thus, the algorithm consists in first calculating \( k_1 \) with \( t_i \), \( y_1 \) and \( f \) as inputs. Thereafter, we increase the step size by \( h/2 \) and calculate \( k_2 \), then \( k_3 \) and finally \( k_4 \). The global error goes as \( O(h^4) \).
However, at this stage, if we keep adding different methods in our main program, the code will quickly become messy and ugly. Before we proceed thus, we will now introduce functions that enbody the various methods for solving differential equations. This means that we can separate out these methods in own functions and files (and later as classes and more generic functions) and simply call them when needed. Similarly, we could easily encapsulate various forces or other quantities of interest in terms of functions. To see this, let us bring up the code we developed above for the simple sliding block, but now only with the simple forward Euler method. We introduce two functions, one for the simple Euler method and one for the force.
Note that here the forward Euler method does not know the specific force function to be called. It receives just an input the name. We can easily change the force by adding another function.
def ForwardEuler(v,x,t,n,Force):
for i in range(n-1):
v[i+1] = v[i] + DeltaT*Force(v[i],x[i],t[i])
x[i+1] = x[i] + DeltaT*v[i]
t[i+1] = t[i] + DeltaT
def SpringForce(v,x,t):
# note here that we have divided by mass and we return the acceleration
return -2*gamma*v-x+Ftilde*cos(t*Omegatilde)
It is easy to add a new method like the Euler-Cromer
def ForwardEulerCromer(v,x,t,n,Force):
for i in range(n-1):
a = Force(v[i],x[i],t[i])
v[i+1] = v[i] + DeltaT*a
x[i+1] = x[i] + DeltaT*v[i+1]
t[i+1] = t[i] + DeltaT
and the Velocity Verlet method (be careful with time-dependence here, it is not an ideal method for non-conservative forces))
def VelocityVerlet(v,x,t,n,Force):
for i in range(n-1):
a = Force(v[i],x[i],t[i])
x[i+1] = x[i] + DeltaT*v[i]+0.5*a
anew = Force(v[i],x[i+1],t[i+1])
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
Finally, we can now add the Runge-Kutta2 method via a new function
def RK2(v,x,t,n,Force):
for i in range(n-1):
# Setting up k1
k1x = DeltaT*v[i]
k1v = DeltaT*Force(v[i],x[i],t[i])
# Setting up k2
vv = v[i]+k1v*0.5
xx = x[i]+k1x*0.5
k2x = DeltaT*vv
k2v = DeltaT*Force(vv,xx,t[i]+DeltaT*0.5)
# Final result
x[i+1] = x[i]+k2x
v[i+1] = v[i]+k2v
t[i+1] = t[i]+DeltaT
Finally, we can now add the Runge-Kutta2 method via a new function
def RK4(v,x,t,n,Force):
for i in range(n-1):
# Setting up k1
k1x = DeltaT*v[i]
k1v = DeltaT*Force(v[i],x[i],t[i])
# Setting up k2
vv = v[i]+k1v*0.5
xx = x[i]+k1x*0.5
k2x = DeltaT*vv
k2v = DeltaT*Force(vv,xx,t[i]+DeltaT*0.5)
# Setting up k3
vv = v[i]+k2v*0.5
xx = x[i]+k2x*0.5
k3x = DeltaT*vv
k3v = DeltaT*Force(vv,xx,t[i]+DeltaT*0.5)
# Setting up k4
vv = v[i]+k3v
xx = x[i]+k3x
k4x = DeltaT*vv
k4v = DeltaT*Force(vv,xx,t[i]+DeltaT)
# Final result
x[i+1] = x[i]+(k1x+2*k2x+2*k3x+k4x)/6.
v[i+1] = v[i]+(k1v+2*k2v+2*k3v+k4v)/6.
t[i+1] = t[i] + DeltaT
The Runge-Kutta family of methods are particularly useful when we have a time-dependent acceleration. If we have forces which depend only the spatial degrees of freedom (no velocity and/or time-dependence), then energy conserving methods like the Velocity Verlet or the Euler-Cromer method are preferred. As soon as we introduce an explicit time-dependence and/or add dissipitave forces like friction or air resistance, then methods like the family of Runge-Kutta methods are well suited for this. The code below uses the Runge-Kutta4 methods.
DeltaT = 0.001
#set up arrays
tfinal = 20 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, v, and x
t = np.zeros(n)
v = np.zeros(n)
x = np.zeros(n)
# Initial conditions (can change to more than one dim)
x0 = 1.0
v0 = 0.0
x[0] = x0
v[0] = v0
gamma = 0.2
Omegatilde = 0.5
Ftilde = 1.0
# Start integrating using Euler's method
# Note that we define the force function as a SpringForce
RK4(v,x,t,n,SpringForce)
# Plot position as function of time
fig, ax = plt.subplots()
ax.set_ylabel('x[m]')
ax.set_xlabel('t[s]')
ax.plot(t, x)
fig.tight_layout()
save_fig("ForcedBlockRK4")
plt.show()
If one has several driving forces, \( F(t)=\sum_n F_n(t) \), one can find the particular solution to each \( F_n \), \( x_{pn}(t) \), and the particular solution for the entire driving force is $$ \begin{equation} x_p(t)=\sum_nx_{pn}(t). \label{_auto61} \end{equation} $$
This is known as the principal of superposition. It only applies when the homogenous equation is linear. If there were an anharmonic term such as \( x^3 \) in the homogenous equation, then when one summed various solutions, \( x=(\sum_n x_n)^2 \), one would get cross terms. Superposition is especially useful when \( F(t) \) can be written as a sum of sinusoidal terms, because the solutions for each sinusoidal (sine or cosine) term is analytic, as we saw above.
Driving forces are often periodic, even when they are not sinusoidal. Periodicity implies that for some time \( \tau \) $$ \begin{eqnarray} F(t+\tau)=F(t). \end{eqnarray} $$
One example of a non-sinusoidal periodic force is a square wave. Many components in electric circuits are non-linear, e.g. diodes, which makes many wave forms non-sinusoidal even when the circuits are being driven by purely sinusoidal sources.
The code here shows a typical example of such a square wave generated using the functionality included in the scipy Python package. We have used a period of \( \tau=0.2 \).
import numpy as np
import math
from scipy import signal
import matplotlib.pyplot as plt
# number of points
n = 500