# Particle in one-dimensional potential well

## Period of oscillations  in potential well 

Dynamics of a particle of mass $m$ moving in one dimension $OX$ is described the Newton equation 

$$m\ddot x =m\dot v = F(x) = -U'(x),$$ 
 
where $F(x)$ is a force acting on tha particle and $U(x)$ is potential energy of the particle. We assume that there are no dissipative forces. Therefore, as we know from the previous section, the total energy $E$ of the particle is conserved, i.e.,  

\begin{equation}
\label{eq:cons_E1d}
\frac{mv^2}{2} + U(x) = E = const.
\end{equation}

If we treat the velocity $v=\frac{dx}{dt}$ as an  independent variable, it is implicit equation of the orbit (the phase curve) in **phase space**  $(x,v)$ of the  particle with  energy $E$. This equation can be rewritten in the form 

$$\frac{m}{2}\left(\frac{dx}{dt}\right)^2+U(x)=E$$

It is a first-order differential equation and can be solved by the method of variable separation. The result reads 

\begin{equation}
\label{eq:t_cons_E1d}
 t=\sqrt{\frac{m}{2}} \; \int_{a}^{b}{\frac{dx}{\sqrt{E-U(x)}}}
\end{equation}


Here, $a$ is the initial position $x(0)=a$ and $b$ is the final position $x(t)=b$ of the particle. The time $t$ is time for moving the particle from the point $a$ to the point $b$ under the condition that $E \ge U(x)$ for all $x \in (a, b)$. 

This is a useful formula which allows to calculate period of oscillations of the particle in a potential well.  

In this section we will consider a class of potentials in the form  

\begin{equation}
\label{eq:Uxn}
U(x)=A |x|^n,
\end{equation}

where $n$ is a positive real number.

These potentials are similar: they are bounded from below,   have only one minimum at $x=0$ and tends to infinity when $x\to \pm \infty$. In such potential the particle motion is bounded and the particle oscillates between two symmetrical positions $x_0$ and $-x_0$ which in general depend on the total energy $E$ and are determined by the equation  

$$U(\pm x_0) = E$$

Because of symmetry, the period $T$ of oscillations can be determined by the equation 

\begin{equation}
\label{eq:T}
 T=4 \; \sqrt{\frac{m}{2}} \; \int_{0}^{x_0}{\frac{dx}{\sqrt{E-U(x)}}}
\end{equation}

This class of potentials is simple but nevertheless analysis of $T$ in dependence of the index $n$ and the total energy $E$ is very interesting and instructive. 

We will use computer algebra and numerical methods to investigate properties of motion is such potential wells.  



```python
load('cas_utils.sage')
```


```python
t = var('t')
m = var('m')
A = var('A')
assume(A > 0)
assume(m > 0)
y = function('y')(t)
de = m*diff(y,t,2) + 2*A*y == 0
showmath( desolve(de, y,ivar=t) )
```

It is an analytical solution of the Newton equation in the case when $n=2$ (a harmonic potential). 

## Particle in potential $x^2$

For $n=2$ the system is  a harmonic oscillator:

$$U(x)=A x^2.$$


```python
#reset()
var('m A x E')
forget()
assume(A > 0)
assume(E > 0)
assume(E,'real')
```

To obtain the integration limit $x_0$  in the formula for the period of oscillations, we must solve the equation:

$$U(x)=E$$

So for the $Ax^2$ potential, we have:


```python
U(A,x) = A*x^2
xextr = solve (U(A=A,x=x)==E,x)
showmath(xextr)
```

These formulas describe the values of the oscillator's extremal positions for a given energy. Let's put them into the formula for $T$:
:


```python
period = 2*sqrt(m/2)*integrate( 1/sqrt(E-U(A,x)),(x,x.subs(xextr[0]),x.subs(xextr[1])))
period = period.canonicalize_radical()
showmath(period)
```

We see that the period $T$ does not depend on energy of the oscillator. It means that it does not depend on the initial conditions because the total energy of the particle depends on them. In turn, it means that it does not depend on the distance  between the points $-x_0$ and $x_0$. It seems to be unusual behavior: time to travel  from $-1$ to $1$ and back is the same as time to travel from $-10000$ to $10000$ and back. In the second case the distance is much,  much longer but time is the same. This is an exceptional property valid only for the harmonic potential! 

## Particle in $|x|^n$ potential

If $n\neq2$, the general formula for the period can be written as:

$$T=4 \sqrt{\frac{m}{2}} \; \int_0^{x_0}\frac{dx}{\sqrt{E- A x^n}}$$

or in the equivalent form: 

$$T=4 \sqrt{\frac{m}{2}}\frac{1}{\sqrt{E}}\int_0^{x_0}\frac{dx}{\sqrt{1-Ax^n/E}}$$

This  integral can be transformed  to a dimensionless form  by substitution 

$$\frac{A}{E}x^n=y^n.$$

It  is in fact a  linear relationship between $x$ and $y$:

$$\left(\frac{A}{E}\right)^{\frac{1}{n}}x=y.$$

Therefore, we can change the integration variable to $y$.  To do this,  we use SAGE to transform the expression under integral  in the following way:


```python
var('dx dy A E x y')
var('n',domain='integer')
assume(n>=0)
assume(A>0)
assume(E>0)
ex1 = dx/sqrt(1-A/E*x^n)
showmath(ex1)
```

and we substitute:


```python
ex2 = ex1.subs({x:(E/A)^(1/n)*y,dx:dy*(E/A)^(1/n)})
showmath( ex2.canonicalize_radical().full_simplify() )
```

Let's take out the expression that depends on the parameters $A$ and $E$:


```python
expr2 = (ex2/dy*sqrt(-y^n + 1)).full_simplify()
showmath( expr2.canonicalize_radical() )
```


```python
prefactor = expr2*sqrt(m/2)*4*1/sqrt(E)
showmath( prefactor.canonicalize_radical() )
```

Finally, we obtain:

$$T=4 \sqrt{\frac{m}{2}}\frac{1}{A^{1/n}} E^{\frac{1}{n}-\frac{1}{2}}\int_0^{y_0}\frac{dy}{\sqrt{1- y^n}}$$


For $n=2$, dependence on $E$ disappears, as we already have seen in the previous case.

We still need to calculate the upper limit  $y_0$ of integration.  In the integral, the upper limit is the position in which the total energy is the potential energy:

$$U(x)=E$$

In this case $$Ax^n=E.$$

By changing the variables we get:


```python
solve( (A*x^n == E).subs({x:(E/A)^(1/n)*y}), y)
```

That is, the integration limit is $y_0=1.$

Therefore the period of oscillations is given by the relation:

$$T=4 \sqrt{\frac{m}{2}}\frac{1}{A^{1/n}} E^{\frac{1}{n}-\frac{1}{2}}\int_0^{1}\frac{dy}{\sqrt{1- y^n}}$$

We note that only for the case $n=2$, the period $T$ does not depend on E (i.e. on initial conditions, i.e. on the distance). In other cases it depends on the total energy $E$, i.e. it depends on initial conditions, i.e. it depends on the distance  between the points $-x_0$ and $x_0$. 

The above equation shows how much time the particle needs to travel the distance for one oscillation in dependence on $E$ and in consequence on the distance:   If energy $E$ is higher then the distance $4x_0$ is longer. 

The scaled integral can be expressed by the beta function of Euler http://en.wikipedia.org/wiki/Beta_function.
We can calculate it:


```python
var('a')
assume(a,'integer')
assume(a>0)
print( assumptions() )
```


```python
integrate(1/sqrt(1-x^(a)),(x,0,1) )
```

We get a formula containing the beta function. It can be evaluated numerically for any values ​​of the $a$ parameter.


```python
(beta(1/2,1/a)/a).subs({a:2}).n()
```

Let's examine this formula numerically. You can use the $\beta$ function, or numerically estimate the integral. This second approach allows you to explore any potential, not just $U(x)=ax^n$.


```python
def beta2(a,b):
    return gamma(a)*gamma(b)/gamma(a+b)

a_list = srange(0.02,5,0.1)
a_list2 = [1/4,1/3,1/2,1,2,3,4,5]

integr_list = [ integral_numerical( 1/sqrt(1-x^a_) ,0,1, algorithm='qng',rule=2)[0] \
               for a_ in a_list ]
integr_list_analytical = [ beta2(1/2, 1/a_)/a_ for a_ in a_list2 ]
```

we obtain some analytically simple formulas:


```python
showmath( integr_list_analytical )
```

Not  we can compare those analytical numbers with numerical results, for example on the plot:


```python
plt_num = list_plot(zip( a_list,integr_list), plotjoined=True )
plt_anal = list_plot(zip( a_list2,integr_list_analytical),color='red')
(plt_num + plt_anal).show(ymin=0,figsize=(6,2))
```

Having an analytical solution, you can examine the asymptotics for large $n$:


```python
var('x')
asympt = limit( beta2(1/2, 1/x)/x,x=oo )
asympt
```


```python
plt_asympt = plot(asympt,(x,0,5),linestyle='dashed',color='gray')
```

Let's add a few points for which the integral takes exact values


```python
l = zip(a_list2[:5],integr_list_analytical[:5])
showmath(l)
```


```python
def plot_point_labels(l):
    p=[]
    for x,y in l:
        p.append(  text( "$("+latex(x)+", "+latex(y)+")$" ,(x+0.1,y+0.2) , fontsize=14,horizontal_alignment='left',color='gray') )
        p.append( point ( (x,y),size=75,color='red' ) )
    return sum(p)
```


```python
some_points = plot_point_labels(l)
```


```python
plt_all = plt_num+plt_anal+plt_asympt+some_points
plt_all.show(figsize=(6,3),ymin=0,ymax=7)
```

## Numerical convergence

The integral  
$$\int_0^1  \frac{dx}{\sqrt{1-x^n}}$$ seems to be divergent for $n:$


```python
showmath( numerical_integral( 1/sqrt(1-x^(0.25)) , 0, 1)  )
```

However, choosing the right algorithm gives the correct result:


```python
a_ = 1/4. # exponent in integral
integral_numerical( 1/sqrt(1-abs(x)^a_) , 0, 1, algorithm='qags')
```

lets check it out with an exact formula:


```python
(beta(1/2,1/a)/a).subs({a:a_}).n()
```

Indeed, we  see that carefull numerical integration gives finite result.

## The dependence of  period on  energy for different $n$.


```python
var('E x n')
def draw_E(n,figsize=(6,2.5)): 
    p = []
    p2 = []
    p.append( plot(abs(x)^n,(x,-2,2),\
                   ymin=0,ymax=4,legend_label=r"$U(x)=|x|^{%s}$" % n )   )
    p.append( plot( (x)^2,(x,-2,2),\
                   color='gray',legend_label=r"$U(x)=x^{2}$",\
                   axes_labels=[r"$x$",r"$U(x)$"]  ))
    
    p2.append(  plot(  4/sqrt(2)*(beta(1/2, 1/n)/n)* E^(1/n-1/2),\
                     (E,0.00,3),ymin=0,ymax=7,axes_labels=[r"$E$",r"$T$"]  ) )
    p2.append(  plot(  4/sqrt(2)*(beta(1/2, 1/2)/2),\
                     (E,0.00,3) ,color='gray' ) )
    
    show( sum(p), figsize=figsize )
    show( sum(p2), figsize=figsize )
    
```


```python
@interact
def _(n=slider([1/4,1/2,2/3,3/4,1,3/2,2,3])):
    draw_E(n)
```


```python
import os
if 'PDF' in os.environ.keys():
    draw_E(1/2)
    draw_E(1)
    draw_E(4)
```

We can plot the dependence of period $T$ on energy (i.e. amplitude) $T(E)$ for different values of $n$. In figure belowe we see that if $n>2$ then oscillations are faster as energy grows. On the other hand if $n<1$, oscillations are getting slower with growing energy.

Another interesting observation is that for potentials with $n>>1$ and $n<<1$ oscillations will become arbitrarily slow and fast, respectively when $E\to0$. For $n>1$ the potential well is *shallow* at the bottom and *steep* far away from the minimum and for $n<1$ the opposite is true.


```python
n_s  = [1/3,2/3, 1, 2, 3] 
plot( [4/sqrt(2)*(beta(1/2, 1/n_)/n_)* E^(1/n_-1/2) \
                 for n_ in  n_s],\
    (E,0.00,2),axes_labels=[r"$E$",r"$T$"],\
     legend_label=['$n='+str(latex(n_))+'$' for n_ in n_s],\
     ymax=7, figsize=(6,3) )
            
```

## Numerical integration of equations of motion 

Here we will investigate numerically period  $T$  and compare with the analytical formula.
First, let's have a closer look how the potential and force behave for different index $n$:



```python
def plot_Uf(n_):
    U(x) = abs(x)^(n_)
    plt = plot( [U(x),-diff( U(x),x)],(x,-1,1),\
               detect_poles='show',ymin=-3,ymax=3,
              legend_label=[r"$U(x)$",'$f(x)$'])
    plt.axes_labels([r"$x$",r"$U(x)=|x|^%s$"%latex(n_)])
    show(plt,figsize=(6,3))
```


```python
@interact
def _(n_=slider(0.1,2,0.1)):
    plot_Uf(n_)
```


```python
if 'PDF' in os.environ.keys():
    plot_Uf(1)
    plot_Uf(2)
    plot_Uf(1/2)
```

We can see that for $n \ge 1$ the force and potential are continuous. If $n=1$ then force has finite jump (discontinuity). Both those cases should not be a problem for numerical integration. 

However for $n<1$ we have infinite force at $x=0$.

#### Problems with numerical integration

There is a possibility that  if the ODE integrator comes too close, it will blow out!

We can fix this problem by softening the  potential singularity by adding small number: 
$$ |x| \to |x| + \epsilon. $$


```python
var('x',domain='real')
var('v t')
eps = 1e-6
U(x) = (abs(x)+eps)^(1/2)
showmath( U.diff(x).expand().simplify() )
```

to make sure that Sage will not leave $x/|x|$ unsimplified we can do:


```python
w0 = SR.wild(0)
w1 = SR.wild(1)
f = -U.diff(x).subs({w0*w1/abs(w1):w0*sign(w1)})
```


```python
showmath( f(x) )
```


```python
ode_pot = [v,f(x)]

t_lst = srange(0,10,0.01)
sol = desolve_odeint(ode_pot,[1,.0],t_lst,[x,v])
```


```python
p = line(zip(t_lst, sol[:,0])) + line(zip(t_lst, sol[:,1]), color='red')
p.axes_labels(['$t$','$x(t),v(t)$'])
p + plot(1,(x,0,10),linestyle='dashed',color='gray')
```

We can evaluate the period $T$ from the trajectory obtained via numerical solutions. For this purpose one might need an interpolation of numerical table returned by `desolve_odeint`:


```python
import numpy as np
def find_period(x,t):
    zero_list=[]
    x0 = x[0]
    for i in range(1,len(x)):
        if x[i]*x[i-1] < 0:
            zero_list.append( - (t[i-1]*x[i] - t[i]*x[i-1])/(x[i-1] - x[i]) )
    lnp = np.array(zero_list)
    return 2*(  (lnp-np.roll(lnp,1))[1:] ).mean()
```


```python
var('x1 x2 t1 t2 a b ')
showmath(  (-b/a).subs(  solve([a*t1+b==x1,a*t2+b==x2],[a,b], solution_dict=True)[0] ) )
```

We find numerically a period of trajectory:


```python
T = find_period( sol[:,0],t_lst)
T
```


```python
assert abs(T-7.54250)<1e-4
```

Exact results  for comparison:


```python
# for n=2 2*pi/sqrt(2)==(2*pi/sqrt(2)).n()
table( [["n","T"]]+[ [n_,((4/sqrt(2)*(beta(1/2, 1/n_)/n_)* E^(1/n_-1/2)).subs({E:1})).n()] 
                    for n_ in [1/4,1/3,1/2,2/3,1,2,3,4,5] ] )

```

## Using the formula for the period to reproduce the trajectory of movement

We take $m=1$ and $A=1$ (then $x=E$), then we can reproduce the trajectory reversing the formula for $T(E)$. 


```python
var('x')
U(A,x) = A*x^2
A = 1/2
E = 1
m = 1.
x1=0.1
showmath( solve(E-U(A,x), x) )
```


```python
t_lst  = [ (sqrt(m/2.)*integrate( 1/sqrt(E-U(A,x)),(x,0,x1)).n(),x1)  \
          for x1 in srange(0,sqrt(2.)+1e-10,1e-2)]
```


```python
point(t_lst ,color='red')+\
 plot(sqrt(2)*sin(x),(x,0,pi),figsize=(6,2))
```

Interestingly, if we known the dependence of $T(E)$ then we can calculate exactly the potential! 


\newpage

7f7f755c9b262e1286dd6f7e81e491b365331074
ipynb
Jupyter Notebook
012-1d_potential_well.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
012-1d_potential_well.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
2022-01-30T16:45:58.000Z
2022-01-30T16:45:58.000Z
012-1d_potential_well.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