# Particle in  multistable potentials


## Phase portrait for a one-dimensional system

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.


## Example: motion in the $U(x) = x^3-x^2$ potential



```python
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$. 


```python
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. 


```python
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. 


## Harmonic  oscillations  limit  for  one-dimensional systems

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$. 


```python
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)
```


```python
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$


```python
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.

## Period of oscillation around potential minimum

We take a particle with  energy by $dE$ larger than minimum of the potential:


```python
dE = 0.01
```


```python
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](https://en.wikipedia.org/wiki/Machine_epsilon).


```python
_, 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:



```python
period = 2*sqrt(m/2.)*\
  integral_numerical(1/sqrt(E0-U(x)) , x1,x2, algorithm='qags')[0]
period
```

Let is plot this situation:


```python
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:


```python
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
```


```python
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$).


```python
omega = sqrt(U(x).diff(x,2).subs(x==xmin.n()))
period_harm = 2*pi.n()/omega
period_harm
```

#### How does period depent on energy?

We can easily compute plot $T(E_0)$:


```python
TonE = [(E_,T(E_)) for E_ in srange(U(xmin)+1e-6,-1e-5,0.001)]
```


```python
line(TonE, figsize=(6,2),gridlines=[None,[period_harm]])
```

#### How does the tracjectory look like? 

Let us integrate numerically equation of motion:



```python
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)$:


```python
line(sol[::10,:],marker='o',figsize=4)
```

Position oscillates with time:


```python
line(zip(t_lst,sol[:,0]),figsize=(6,2))
```

## Time to reach the hill

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:


```python
var('E')
assume(E>0)
integrate(-1/sqrt(E+x^2),x,1,0)
```

This result is divergent for $E\to0$:


```python
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:


```python
U1(x) = -abs(x)
E0 = 0
```

we can plot velocity and potential:


```python
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:


```python
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. 


```python
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
```


```python
t_hill(9.1)
```


```python
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))
```


```python
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))
```

#### Experiment with Sage!
Examine in a similar way the system corresponding to the movement in the $U(x) = -\cos(x)$ potential - this is a physical pendulum.

\newpage

d8648c0faf3175a2b2d31304aa7d56a2a271123b
ipynb
Jupyter Notebook
013-1d-multistable.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
013-1d-multistable.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
2022-01-30T16:45:58.000Z
2022-01-30T16:45:58.000Z
013-1d-multistable.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
2020-11-15T08:26:14.000Z
2022-02-12T13:07:16.000Z
Qwen/Qwen-72B
1. YES
2. YES