The Newton equation for a point particle in a one-dimensional potential $U(x)$ can be written as a set of two first-order differential equations:
$$ \begin{cases} \dot x = v \\ \dot v = - \frac{1}{m}U'(x) \end{cases} $$We can draw a phase portrait, i.e. parametric solutions $(x(t),v(t))$ and a vector field defined by the right hand sides of equations in the $(x,v)$-phase space.
var('v')
m = 1
U(x) = x^3-x^2
xmax,xmin = sorted([s_.rhs() for s_ in solve(U.diff(x)==0,x)])
Emin = U(xmin)
Etot = 1/2*m*v^2 + U(x)
plot(U(x),(x,-0.4,1.1),figsize=4) +\
point([xmin,U(xmin)],color='red',size=40)+\
point([xmax,U(xmax)],color='red',size=40)
The potential $U(x) =x^3-x^2$.
pkt = point((xmin,0),size=20,color='black')
plt = sum([ implicit_plot(Etot==E0,(x,-1/2,1.2),(v,-1,1),color='blue')\
for E0 in srange(Emin,0.0,0.02)])
plt +=implicit_plot(Etot==0,(x,-1/2,1.2),(v,-1,1),color='black') +pkt
plt +=sum([ implicit_plot(Etot==E0,(x,-1/2,1.2),(v,-1,1),color='green')\
for E0 in srange(0.02,-2*Emin,0.04)])
plt.show()
The phase curves $(x, v)$ depends on initial conditions $(x_0, v_0)$. There are two types of phase curves: closed (periodic motion of the particle in a bounded interval) and open (the motion is unbounded: the particle can escape to -infinity or can return from -infinity.
vector_field = vector([v,-U.diff(x)])
plt + plot_vector_field(vector_field.normalized(),(x,-1/2,1.2),(v,-1,1))
The vector field shows direction of motion of the particle on the $(x, v)$-plane.
Consider a conservative one-dimensional system. In this case, the force $f(x)$ can always be represented as gradient of the potential $U(x)$, namely, $$f(x) = -\frac{\partial U(x)}{dx}.$$ Consider a certain potential that has a minium at some point $x_0$. The necessary condition for minimum of the function is its zero first derivative at this point. Let's expand the potential in the Taylor series around the minimum. We obtain: $$ U(x) = U(x_0) + \underbrace{U'(x_0)( x-x_0)}_{=0}+\frac{1}{2} U''(x_0)(x-x_0)^2+...$$ For small deviation from the minimum this series can be approximated by the function $$U(x) = \frac{1}{2} k (x-x_0)^2,$$
The Newton equation for such motion is as follows:
$$m \ddot x = m a = F = -U'(x) = -k (x-x_0)$$For the new variable $y=x-x_0$ it takes the form
$$m \ddot y =-ky$$This is the already known equation for the harmonic oscillator with the shifted equilibrium point $x_0$.
Now, let us return to the system with the potential $U(x) = x^3-x^2$.
var('x v')
Etot = 1/2*v^2 + U(x)
Elin = 1/2*v^2 + U(xmin)+1/2*U.diff(x,2).subs(x==xmin)*(x-xmin)^2
show(Etot)
show(Elin)
Emin = Etot(x=xmin,v=0)
Emin
Let's have a look at the trajectories for the exact system with $U(x) = x^3-x^2$ and the linearized system with $U(x)=(1/2) x^2$. The blue line below is a separatrix - i.e. a solution with $E=0$
plt = sum([ implicit_plot(Etot==E0,(x,.4,.91),(v,-.3,.3),color='red') \
for E0 in srange(Emin+1e-3,-0.1,0.005)])
plt += implicit_plot(Etot==0.00,(x,0,1.1), (v,-.6,.6),color='blue')
plt_lin =sum([ implicit_plot(Elin==E0+1e-3,(x,.4,.91),(v,-.3,.3),color='gray') \
for E0 in srange(Emin,-0.1,0.005)])
plt+plt_lin
For larger ones, there is a growing discrepancy:
- for the nonlinear system, above certain energy, there are open trajectories - motion in an linearized system is always an ellipse. The period does not depend on the amplitude.
We take a particle with energy by $dE$ larger than minimum of the potential:
dE = 0.01
E0 = U(xmin)+dE
E0
In order to analyze eqaution of motion, we need to solve: $$U(x)=E_0$$ and obtain the extreme values of position to the left and right from potentiall minimum. In our case it con be done analytically.
However, note that it requires to neglect imaginary part, which is approximated small value of the order of machine epsilon.
_, x1,x2 = sorted( [s_.rhs().n().real() for s_ in solve(U(x)==E0,x)])
x1,x2
Now, having $x_{1,2}$ we can integrate numerically equation for period (\ref{eq:1d_TE}) and obtain:
period = 2*sqrt(m/2.)*\
integral_numerical(1/sqrt(E0-U(x)) , x1,x2, algorithm='qags')[0]
period
Let is plot this situation:
U(x) = x^3-x^2
xmax,xmin = sorted([s_.rhs() for s_ in solve(U.diff(x)==0,x)])
plot(U(x),(x,-0.4,1.1),figsize=4,gridlines=[None,[U(xmin),E0]])+\
point([xmin,U(xmin)],color='red',size=40)+\
point([xmax,U(xmax)],color='red',size=40)+\
point([x1,U(x1)],color='green',size=40)+\
point([x2,U(x2)],color='green',size=40)
We can write a function which computes the period based on previous steps:
def T(E0):
m = 1
_, x1,x2 = sorted( [s_.rhs().n().real() \
for s_ in solve(U(x)==E0,x)])
integral, error = \
integral_numerical(1/sqrt(E0-U(x)), x1,x2, algorithm='qags')
period = 2*sqrt(m/2.) * integral
return period
period_num = T(U(xmin)+dE)
period_num
We known exactly the period of oscillation for small neiborhood of potential minimum. It is simply period of the harmonic oscillator with frequecy $\omega=\sqrt{U''(x_{min})}$ (for $m=1$).
omega = sqrt(U(x).diff(x,2).subs(x==xmin.n()))
period_harm = 2*pi.n()/omega
period_harm
We can easily compute plot $T(E_0)$:
TonE = [(E_,T(E_)) for E_ in srange(U(xmin)+1e-6,-1e-5,0.001)]
line(TonE, figsize=(6,2),gridlines=[None,[period_harm]])
Let us integrate numerically equation of motion:
t_lst = srange(0, period_num, 0.01, include_endpoint=True)
sol = desolve_odeint([v,-U.diff(x)],[x2,.0],t_lst,[x,v])
We have in phase space $(x,\dot x)$:
line(sol[::10,:],marker='o',figsize=4)
Position oscillates with time:
line(zip(t_lst,sol[:,0]),figsize=(6,2))
The top of the potential hill is at the beginning of the coordinate system $(x,E)$. Then we examine the limit $E\to0$ boundary
Near zero, we can approximate the potential by the reverse parabola. Then the time to reach the hill from a certain point (for example $x=1$) reads:
var('E')
assume(E>0)
integrate(-1/sqrt(E+x^2),x,1,0)
This result is divergent for $E\to0$:
limit( arcsinh(1/x),x=0)
It means that time needed to climb a hill with just enough kinetic energy is infinite. It is valid only for potential hills which have zero derivative at the top. On the other hand for potential barriers which do not have this property, for example, $U(x)=-|x|$, the particle can reach the top with just enought energy in finite time.
Let analyze it:
U1(x) = -abs(x)
E0 = 0
we can plot velocity and potential:
plot([U1(x), sqrt(2*(E0 - U1(x)))],(x,-1,1),figsize=4)
the time of travel from $x=-1$ to $x=0$ is given by:
\begin{equation} \label{eq:1d_TE} t=\sqrt{\frac{m}{2}} \; \int_{-x1}^{x1}{\frac{dx}{(\sqrt(E-U(x)}} \end{equation}which in this case is:
sqrt(m/2.)*integrate(1/sqrt((E0- U1(x))),x,-1,0).n()
In the case of potentials which behave like $|x|^\alpha$, for $\alpha>1$ we can calculate time of travel if we particle total energy is by $dE$ larger than potential barrier.
def t_hill(E0):
m = 1
x2, = [s_.rhs().n().real() for s_ in solve(U(x)==E0,x) if s_.rhs().n().imag().abs()<1e-6]
integral, error = \
integral_numerical(1/sqrt(E0-U(x)), 0,x2, algorithm='qags')
m = 1
period = 2*sqrt(m/2.) * integral
return period
t_hill(9.1)
import numpy as np
t_E = [(E_,t_hill(E_)) for E_ in np.logspace(-6,1,120)]
line(t_E,axes_labels=['$E$','$t$'], figsize=(6,2))
t_E = [(log(E_),t_hill(E_)) for E_ in np.logspace(-6,16,120)]
line(t_E,axes_labels=['$\log_{10} E$','$t$'], figsize=(6,2))
Examine in a similar way the system corresponding to the movement in the $U(x) = -\cos(x)$ potential - this is a physical pendulum.
\newpage