least_action.ipynb

Kernel: Python 3 (ipykernel)

Least Action

Load modules

%pylab inline

Populating the interactive namespace from numpy and matplotlib

import numpy as np
import scipy.optimize
import pandas as pd
global g  
g=9.8 #m/s^2

Geometry interpretation

Following the geometry theory developed here, we will try to define something called the Action for one small segment of the free fall movement in one-dimension.

For that we need the experimental data consisting on the height of an object of mass $m$ in free fall, and the height $x_i$, for each time $t_i$. This data would be fitted by a polynomial of degree two, as displayed in the figure for one of the fitted segments of the plot of $x$ as a function of $t$. We take the origin of the coordinates at ground level. For each segment we can calculate an average kinetic energy, $T$, and an averge potential energy, $V$, in the limit of $\Delta t=t_2-t_1$ small. From the figure

$$\begin{align} T=\frac12 m v^2\approx &\frac12 m\left(\frac{x_2-x_1}{t_2-t_1}\right)^2\,,& V=mgh\approx& m g \frac{x_2+x_1}{2}\,. \end{align}$$

We can then reformulate the problem of the free fall in the following terms. From all the possible curves that can interpolate the points $(t_1,x_1)$ and $(t_2,x_2)$, which is the correct one?.

The answer obtained by Leonhard Euler [1] can be obtained from the definition of the function "Lagrangian" $$L(t)=T(t)-V(t)\,.$$

With this, we can build the "Action", by integrating the Lagrangian function between the points $(t_1,x_1)$ and $(t_2,x_2)$ as $$S=\int_{t_1}^{t_2} L\, \operatorname{d}t\,.$$

which give us a numerical result with units of energy multiplied by time (J$\cdot$s). What is worth noticing is that if we relax the definition of $V$ and allows for any $h$, but keeping the initial and final points fixed, we can calculate many Action values. This is illustrated in the figure for blue dotted ($S_1$), solid ($S_{\text{min}}$) and dashed ($S_2$)lines. But only the height, $(x_1+x_2)/2$ , associated with the real physical path, has a minimum value for the Action, $S_{\text{min}}$!

In fact, for one segment of the action between $(t_1,x_1)$, and $(t_2,x_2)$, with $\Delta t$ sufficiently small such that $L$ can be considered constant, we have $$\begin{darray}{rcl} S_{\text{min}}&=&\int_{t_1}^{t_2} L dt \\ &\approx& \left[\frac12 m v^2-m g h \right]\Delta t\\ &\approx& \left[\frac12 m\left(\frac{x_2-x_1}{t_2-t_1}\right)^2-m g \frac{x_2+x_1}{2} \right](t_2-t_1) \end{darray}$$ that corresponds to Eq. (11) of Am. J. Phys, Vol. 72(2004)478: http://www.eftaylor.com/pub/Symmetries&ConsLaws.pdf

Least Action method: The least action method consist in the following steps illustrated in the figure

1. Fix the initial and final point of the movement. Example the initial time and the height from which a body is launched upwards, $(t_{\text{ini}},x_{\text{ini}})$, and the final time and height $(t_{\text{end}},x_{\text{end}})$.

2. Divide the problem in small segments of $\Delta t=t_{i+1}-t_i$.

3. Build many paths with the initial and final point fixed and calculates the Action in each segment, $S_i$, for all the paths

4. Choose the minimal Action for each segment, $S_{\text{min}}^i$, and rebuild the full path which minimizes the Action in each segment. This is the physical trajectory!

Code implementation

The Action

We define the Action $S$ such of an object of mass $m$ throw vertically upwards from $x_{\hbox{ini}}$, such that $t_{\hbox{end}}$ seconds later the object return to a height $x_{\hbox{end}}$, as $$\begin{align} S=&\int_{t_{\hbox{ini}}}^{t_{\hbox{end}}} L\, {\rm d}t \\ =& \sum_i S_i\\ =&\sum_i L_i \Delta t\,. \end{align}$$

Example

Consider an object of mass $m=0.2~$Kg throw vertically upwards from an intial height of zero, and a flight time of $3~$s

1. Calculates the action for an intermediate set of points with $\Delta t=0.5$ s and an arbitrary path.

Activity: Change the path to any other value and report with the same initial and final values of zero and report the Action calculations by the chat

tini=0
tend=3
xini=0
xend=0

x=np.array([xini,10,15.,22.,18.,5.,xend])#.size

x

array([ 0., 10., 15., 22., 18.,  5.,  0.])

$$\begin{darray}{rcl} S_{i} &\approx& \left[\frac12 m\left(\frac{x_{i+1}-x_i}{t_{i+1}-t_i}\right)^2-m g \frac{x_{i+1}+x_i}{2} \right](t_{i+1}-t_i)\\ &=& \left[\frac12 m\left(\frac{\Delta x_i}{\Delta t}\right)^2-m g \frac{x_{i+1}+x_i}{2} \right]\Delta t \end{darray}$$$$S=\sum S_i$$

#Initial points: x_i
x[:-1]

array([ 0., 10., 15., 22., 18.,  5.])

#final points: x_{i+1}
x[1:]

array([10., 15., 22., 18.,  5.,  0.])

Δx=x[1:]-x[:-1]
Δx

array([ 10.,   5.,   7.,  -4., -13.,  -5.])

x[:-1].size

6

Δt=tend/x[:-1].size
Δt

0.5

def S(x,tend=3.,m=0.2,xini=0.,xend=0.):
    """
    Calculate the Action of an object of of mass 'm' throw vertically upward from 
       'xini', such that 'tend' seconds later the object return to a height 'xend'.
       Delta t must be constant.
       
    The defaults units for S are J.s   
    """
    x=np.asarray(x)
    Δt=tend/x[:-1].size
    #Fix initial and final point
    x[0]=xini
    x[-1]=xend
    return ( (0.5*m*(x[1:]-x[:-1])**2/Δt**2-0.5*m*g*(x[1:]+x[:-1]) )*Δt).sum()

$$\begin{darray}{rcl} S_i&\approx& \left[\frac12 m\left(\frac{x_2-x_1}{t_2-t_1}\right)^2-m g \frac{x_2+x_1}{2} \right](t_2-t_1) \end{darray}$$

x

array([ 0., 10., 15., 22., 18.,  5.,  0.])

S(x)

8.199999999999996

x=[ 0., 10., 15., 18., 15.,  10,  0.]
S(x)

-13.040000000000008

-0.9600000000000097

t=[i*Δt for i in range(x.size) ]

t

[0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]

df=pd.DataFrame({'t':t,'x':x})

df

plt.plot(df.t,df.x,'ro')
plt.plot(df.t,df.x)
plt.xlabel('$t$')
plt.ylabel('$x(t)$')

Text(0, 0.5, '$x(t)$')

S(df.x) # J.s

8.199999999999996

2. Fit the points with a polynomial of degree 2 and calculates the Action for 20 points along the fitted curve but keeping the same initial and final points!

coeffs=np.polyfit(df.t,df.x,2)
P=np.poly1d(coeffs)
print(P)

        2
-8.905 x + 26.21 x - 0.381

t=np.linspace(0,3,20)

x=P(t)
x[0]=xini
x[-1]=xend
x

array([ 0.        ,  3.53614299,  7.00923361, 10.03831948, 12.62340061,
       14.76447698, 16.46154861, 17.71461549, 18.52367762, 18.888735  ,
       18.80978763, 18.28683551, 17.31987864, 15.90891703, 14.05395067,
       11.75497955,  9.01200369,  5.82502308,  2.19403773,  0.        ])

plt.plot(t,x,'ko')
plt.plot(t,x)

[<matplotlib.lines.Line2D at 0x7f37165ab630>]

S(x)

-10.04616480633495

Exercise

Try to find a large action (larger than 10 J$\cdot$s) and the minimal possible action

Solution to a free fall problem

$$\begin{align} x(t)=&-\frac{1}{2} g t^2+v_0 t \,, & v=&-g t +v_0 \end{align}$$

At maximum height $t=t_{\text{end}}/2$ and $$\begin{align} x_{\text{max}}=&-\frac{1}{2} g (t_{\text{end}}/2)^2+v_0 \frac{t_{\text{end}}}{2} \,, & 0=&-g \frac{t_{\text{end}}}{2} +v_0 \end{align}$$ From the second equation we get $$v_0=\frac{g t_{\text{end}}}{2}=14.7 \text{m/s}\,,$$ and $$\begin{align} x_{\text{max}}=&-\frac{1}{2} g \left(\frac{t_{\text{end}}}{2}\right)^2+ g \frac{t_{\text{end}}}{2} \frac{t_{\text{end}}}{2} \nonumber\\ =&g\frac{ t_{\text{end}}^2 }{8}=11.025\ \text{m}\,. \end{align}$$

Activity: Find the same solution by fitting the polynomial of degree 2 that go trough the points:

• $(t_0,x_0)=(0,0)$,

• $(t_{\rm end}/2,x_{\rm max})=(1.5,11.025)$,

• $(t_{\rm end},x_{\rm end})=(3,0)$, in s and m respectively.

df=pd.DataFrame( {'t':[0,1.5,3],
                  'x':[0,11.025,0]
     })
df

#df.loc[1,'x']=11.025 
coeff=np.polyfit(df.t,df.x,deg=2)
x=np.poly1d(coeff,variable='t')
print(x)

      2
-4.9 t + 14.7 t + 8.205e-15

t=np.linspace(df.t.min(),df.t.max(),20)
plt.plot(df.t,df.x,'ro')
plt.plot(t,x(t))
plt.grid()

t

array([0.        , 0.15789474, 0.31578947, 0.47368421, 0.63157895,
       0.78947368, 0.94736842, 1.10526316, 1.26315789, 1.42105263,
       1.57894737, 1.73684211, 1.89473684, 2.05263158, 2.21052632,
       2.36842105, 2.52631579, 2.68421053, 2.84210526, 3.        ])

x(t)

array([8.20464080e-15, 2.19889197e+00, 4.15346260e+00, 5.86371191e+00,
       7.32963989e+00, 8.55124654e+00, 9.52853186e+00, 1.02614958e+01,
       1.07501385e+01, 1.09944598e+01, 1.09944598e+01, 1.07501385e+01,
       1.02614958e+01, 9.52853186e+00, 8.55124654e+00, 7.32963989e+00,
       5.86371191e+00, 4.15346260e+00, 2.19889197e+00, 1.88627818e-14])

S(x(t))

-21.549141274238227

Activity:

1. Calculates $$S=\int_{t_{\rm min}}^{t_{\rm max}} L(t)\operatorname{d}t$$ for the physical trajectory of a particle in free fall of $m=0.2\$kg, where $$L(t)=\frac{1}{2}m v^2(t)-mgx(t)\,.\qquad\qquad (1)$$ Hint: $$F(t)=\int L(t)\operatorname{d}t$$ $$S=F(t_{\rm max})-F(t_{\rm min})$$

2. Calculates also the instataneous energy $$E(t)=\frac{1}{2}m v^2(t)+mgx(t)$$ for  several times.

print(x)
v=x.deriv()
print(v)

      2
-4.9 t + 14.7 t + 8.205e-15
 
-9.8 x + 14.7

print('x(t) [m]:\n',x)
print("="*45)
m=0.2
g=9.8

#Use derivative of numpy polynomial
v=x.deriv()

# Implement eq. (1) and (2)
L=0.5*m*v**2-m*g*x
F=L.integ() # Analytical integration
S=F( df.t.max() ) - F( df.t.min() ) 
print('S={:.1f} J.s'.format(S))
print("="*45)

#Implement eq. (3)
E=0.5*m*v**2+m*g*x
Deltat=0.2
t=np.arange(df.t.min() , df.t.max()+Deltat,Deltat)
print('instantaneous energy [J]:\n {}'.format(E(t)))

x(t) [m]:
       2
-4.9 t + 14.7 t + 8.205e-15
=============================================
S=-21.6 J.s
=============================================
instantaneous energy [J]:
 [21.609 21.609 21.609 21.609 21.609 21.609 21.609 21.609 21.609 21.609
 21.609 21.609 21.609 21.609 21.609 21.609]

Activity: Changes the maximum height for the previous activity and check what happens with the instantaneous energy.

Comparision with non-physical trajectories

---

I good approximation to the integral is: $$S=\int_{t_{\rm min}}^{t_{\rm max}} L(t)\operatorname{d}t\approx \sum_i \left[\frac{L(t_{i+1})+L(t_i)}{2}\right]\Delta t$$

print('Δt={}'.format( t[1:]-t[:-1] )  )

Δt=[0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2]

print( 'S={:.1f} J.s'.format( 0.5*( L(t)[1:]+L(t)[:-1] ).sum()*Deltat ) )

S=-21.2 J.s

Or even more general

def Saprx(L,t):
    return   0.5*( L(t)[1:]+L(t)[:-1] )*(  t[1:]-t[:-1]  )  

print( 'S={:.1f} J.s'.format(  Saprx(L,t).sum() ) )

S=-21.2 J.s

Ss=[]
Es=[]
m=0.2
g=9.8
Deltat=0.2
#t=np.linspace(df.t.min() , df.t.max(),10)
t=np.arange(df.t.min() , df.t.max()+Deltat,Deltat)
xmax=[6,11.025,15]; x0=0;xend=0
ls=['r--','k-','b-.']
for i in range(len(xmax) ):
    df=pd.DataFrame({'t':[t.min(),t.mean(),t.max()],'x':[x0,xmax[i],xend]})
    coeffs=np.polyfit(df.t,df.x,2)
    x=np.poly1d(coeffs,variable='t')
    v=x.deriv()
    L=0.5*m*v**2-m*g*x
    F=L.integ() # Analytical integration 
    E=0.5*m*v**2+m*g*x
    S=F( df.t.max() ) - F( df.t.min() )
    print('S={:.1f}  J.s'.format( S ) )
    Ss.append( Saprx(L,t)  )
    Es.append( E(t) )
    #Plot
    if ls:
        plt.plot(t,x(t),ls[i],label='S={:.1f} J$\cdot$s'.format( S  )   )

if ls:
    plt.legend(loc='best',fontsize=12)    
    plt.xlabel('$t$ [s]',size=20)
    plt.ylabel('$x$ [m]',size=20)

S=-17.1  J.s
S=-21.6  J.s
S=-18.8  J.s

Only the physical trajectory conserves energy

Es

[array([ 6.4       ,  7.73404444,  8.87751111,  9.8304    , 10.59271111,
        11.16444444, 11.5456    , 11.73617778, 11.73617778, 11.5456    ,
        11.16444444, 10.59271111,  9.8304    ,  8.87751111,  7.73404444,
         6.4       ]),
 array([21.609, 21.609, 21.609, 21.609, 21.609, 21.609, 21.609, 21.609,
        21.609, 21.609, 21.609, 21.609, 21.609, 21.609, 21.609, 21.609]),
 array([40.        , 37.36177778, 35.10044444, 33.216     , 31.70844444,
        30.57777778, 29.824     , 29.44711111, 29.44711111, 29.824     ,
        30.57777778, 31.70844444, 33.216     , 35.10044444, 37.36177778,
        40.        ])]

Approximated result

Ss[0].sum(),Ss[1].sum(),Ss[2].sum()

(-16.958577777777823, -21.224840000000114, -18.1831111111111)

Activity: Check that the Action of a trajectory with a polynomial of order $n$ which fits $n+1$ points with the same extreme values, gives an Action that is larger than the previous minimun

df

ndf=df.append({'t':0.75,'x':5},ignore_index=True)
ndf=ndf.append({'t':2.25,'x':2.7},ignore_index=True)
ndf=ndf.sort_values('t')
ndf

Ss=[]
Es=[]
m=0.2
g=9.8
Deltat=0.2
#t=np.linspace(df.t.min() , df.t.max(),10)
t=np.arange(df.t.min() , df.t.max()+Deltat,Deltat)
xmax=[11.025]; x0=0;xend=0
ls=['r--','k-','b-.']
for i in range(len(xmax) ):
    df=pd.DataFrame({'t':[t.min(),t.mean(),t.max()],'x':[x0,xmax[i],xend]})
    coeffs=np.polyfit(df.t,df.x,2)
    x=np.poly1d(coeffs,variable='t')
    v=x.deriv()
    L=0.5*m*v**2-m*g*x
    F=L.integ() # Analytical integration 
    E=0.5*m*v**2+m*g*x
    S=F( df.t.max() ) - F( df.t.min() )
    print('S={:.1f}  J.s'.format( S ) )
    Ss.append( Saprx(L,t)  )
    Es.append( E(t) )
    #Plot
    if ls:
        plt.plot(t,x(t),ls[i],label='S={:.1f} J$\cdot$s'.format( S  )   )

if ls:
    plt.legend(loc='best',fontsize=12)    
    plt.xlabel('$t$ [s]',size=20)
    plt.ylabel('$x$ [m]',size=20)
    

S=-21.6  J.s

Solution with minimization

To be seen in Minimization

References

[1] L. Euler, M´emoires de l'Acad´emie des Sciences de Berlin 4, 1898 (1748), republished in Ref.2, p. 38-63; L. Euler, M´emoires de l'Acad´emie des Sciences de Berlin 7, 169 (1751), republished in Ref. 2, p. 152. For a recent historical review see: Dias, Penha Maria Cardozo. (2017). Leonhard Euler's “principle of mechanics” (an essay on the foundations of the equations of motion). Revista Brasileira de Ensino de Física, 39(4), e4601. Epub May 22, 2017.https://doi.org/10.1590/1806-9126-rbef-2017-0057

[2] Leonhardi Euleri Opera Omnia, serie secunda, v. V, edited by J.O. Fleckenstein (Societatis Scientiarum Naturalium Helveticæ, Geneva,1957)