💻 Open in Colab#

Ashwin Kharat

1
#-- This will reset the runtime and clear all variable.
2
%reset -f
3
#A new Begining
4

5
#-- importing all the necessary Libraries
6
import numpy as np
7
from matplotlib import pyplot as plt
8
from matplotlib import animation, rc
9
import ipywidgets as widgets
10
from IPython.display import display, HTML
11
## --
12

13
#-- Custom functions to make life easy.
14
# g3 represent Group 3
15
# lable class creates a label object which can be used with new plot and multiPlot functions
16
# lable format is
17
#1.Legend of the plot that u want to make
18
#2.X-axis label,  3. Y-axis label
19
#3.style of line plotted
20
#this can all be skipped and a default plot is made with the default label
21
#newPlot makes a fresh new Figure and plot taking the g3_lable
22
#formate is (x values , y values , g3_lable object )
23
# multiPlot can be used repeatedly to plot on the same figure.
24
#formate is (x values , y values , g3_lable object ).
25

26
class label():
27
  def __init__(self, ll = None ,lx = 'x',ly ='y',style = '-'):
28
    self.ll = ll    # the lable for legend
29
    self.lx = lx    # x axis
30
    self.ly = ly    # y axis
31
    self.style = style
32

33
class g3:
34

35
  def newPlot(x, y, L = label() ):
36
    fig = plt.figure()
37
    easy = fig.add_subplot(111)
38
    easy.plot( x, y,  L.style, label = L.ll )
39
    easy.legend()
40
    plt.xlabel( L.lx )
41
    plt.ylabel( L.ly )
42

43
  def multiPlot(x, y, L = label() ):
44
    plt.plot( x, y, L.style, label = L.ll ,alpha=0.5)
45
    plt.legend()
46
    plt.xlabel( L.lx )
47
    plt.ylabel( L.ly )
48
#--

##Problem 1 : Implementation of Euler and Euler-Cromer ODE.

Part A#

Write a function for solving the Newton’s laws of motion $\frac{d^2 x}{dt^2} = a(x, \dot{x}, t),$ using Euler and Euler-Cromer algorithms. As shown in class, this 2nd order ODE can be broken into a set of 2 coupled 1st order ODEs:

$\frac{d}{dt} \begin{pmatrix}x \\ v \end{pmatrix} = \begin{pmatrix} {v} \\ a(x,v,t) \end{pmatrix}$

The $\texttt{EulerODE}(\, )$ function should take following as input:

a. a list or array containing initial condition for $x$ and $v$ i.e. $[x_0, v_0]$ .

b. a list or array with starting time and end time i.e. $[t_0, t_{max}]$ .

c. The function to calculate the right-hand side of the above equation. It can be defined with a Python function as follows:

$\texttt{rhs} \; \Big(x_{k}, v_{k}, t_{k} \Big):$

and return a list containing the velocity and acceleration at time $t_k$ :

$\texttt{return} \qquad \Big[ v[t_k], \, a[t_k] \Big]$

Hence the functional form for the $\texttt{EulerODE}$ should be as follows:

$\texttt{EulerODE} \Big(\texttt{initCondn}, \; \texttt{tRange}, \; rhs\Big)$ where $\texttt{initCondn} \equiv [x_0, v_0]$ is a list containing the initial conditions at $t=t_0$ and $\texttt{tRange} \equiv [t_0, t_f]$ is the list containing the time-range to integrate the equations.

Make a similar function for Euler-Cromer algorithm $\texttt{EulerCromerODE}(\, )$ .

Solution part A#

1
# Just a BluePrint.
2

3
def rhs(xk,vk,tk):
4
#This is the general definition for any 2nd order vectorised ODE.
5
#we are calling another function 'a(xk,vk,tk)' which has the definition for acceleration written.
6
  return [vk , a(xk,vk,tk)]
7

8

9

10
def EulerODE(initCondn, tRange, rhs):
11
  # we could have given a value to dt but wanted user to have a control over it.
12
  dt = tRange[2]
13
  time = np.arange(tRange[0], tRange[1] + dt , dt)
14

15
  position = np.zeros(len(time))
16
  velocity = np.zeros(len(time))
17

18
  position[0] = initCondn[0]
19
  velocity[0] = initCondn[1]
20

21
  for i in range(0, len(time)-1):
22
    position[i+1] = position[i] + rhs(position[i],velocity[i], time[i] )[0] * dt
23
    velocity[i+1] = velocity[i] + rhs(position[i],velocity[i], time[i] )[1] * dt
24

25
  return position, velocity
26

27

28
def EulerCromerODE(initCondn, tRange, rhs):
29

30
  dt = tRange[2]
31
  time = np.arange(tRange[0], tRange[1] + dt , dt)
32
  position = np.zeros(len(time))
33
  velocity = np.zeros(len(time))
34

35
  position[0] = initCondn[0]
36
  velocity[0] = initCondn[1]
37

38
  for i in range(0, len(time)-1):
39
    velocity[i+1] = velocity[i] + rhs(position[i],velocity[i], None )[1] * dt
40
    position[i+1] = position[i] + rhs(position[i+1],velocity[i+1], None )[0] * dt
41

42
  return position, velocity

1
# An example code is written for a spring mass system with the format given above.
2
def a(xk,vk,t):
3
  mass = 1 #kg
4
  c = 0 # damping constant taking zero for simple case / demo purpose
5

6
  k = 1 # spring constant
7
  return -1/mass*(c*vk + k* xk )
8
# -- Initialization
9
time_0 = 0
10
time_max = 50
11
time_step = 1/10 # 10th part of a second
12
position_0 = 5 # in meters difference from rest position
13
velocity_0 = 0        # angular velocity
14
#--
15

16
initCondn = [position_0, velocity_0]
17
tRange = [time_0, time_max, time_step]
18
dt = tRange[2]
19
time = np.arange(tRange[0], tRange[1] + dt , dt)
20

21
# calling euler fn to solve for the given initial condition.
22
position, velocity = EulerODE(initCondn, tRange, rhs)
23

24
lable = label('trajectory in phase space','position (m)','velocity ( m/s )','--')
25
g3.newPlot(position, velocity,lable)
26
plt.plot(position[0], velocity[0],'*r')
27
#plt.legend()
28
plt.title('Phase Space by Euler method')
29
plt.show()
30

31
lable = label('time vs velocity ','time (sec)','velocity ( m/s ) ','--')
32
g3.newPlot(time, velocity,lable)
33

34
lable = label('oscillation amplitude','time (sec)','position (m)','--')
35
g3.newPlot(time, position  ,lable)

png

Simulation of SHO : Euler.#

1
x_data = []
2
y_data = []
3
l = 10 # length of unstretched spring.
4

5
fig, ax = plt.subplots(figsize=(6,6))
6
ax = plt.axes(xlim=(-l, l*3),ylim=(-l*0.8, l*0.8 ))
7
line, = ax.plot([], [])
8
dot, = ax.plot([], [], 'sk', markersize = 20 )
9
plt.plot( np.zeros(10),np.linspace(0,2,10), '-k')
10
plt.plot( np.linspace(-10,30,10) ,np.zeros(10), '-b')
11

12
plt.xlabel('x (meter)')
13
plt.ylabel('y (meter)')
14
plt.title('Simulation using Euler method')
15

16
def animation_frame(i):
17

18
  x_data = np.linspace(0, (l + position[i]), 50)
19
  temp = 0.5 * np.sin(x_data * (l*2*np.pi)/(l + position[i]))
20
  y_data = temp + 0.5 # just putting spring in place.
21

22
  line.set_xdata(x_data)
23
  line.set_ydata(y_data)
24

25
  dot.set_ydata([0.5])
26
  dot.set_xdata([(l + position[i])])
27

28
  return line, dot
29

30
# simple plot is just the initial condition.
31
# Animination is a video.
32

33
animation_vid = animation.FuncAnimation(fig, func = animation_frame, frames = np.arange(0, len(time)-1 , 1) , interval= dt*1000 )
34
HTML(animation_vid.to_html5_video())

png

The above numerical solution is using the euler method which is why we are getting an error. The solution shows that an extra energy is pumped into the system which is not.

1
position, velocity = EulerCromerODE(initCondn, tRange, rhs)
2

3
lable = label('trajectory in phase space','position (m)','velocity ( m/s )','--')
4
g3.newPlot(position, velocity,lable)
5
plt.plot(position[0], velocity[0],'*r')
6
#plt.legend()
7
plt.title('Phase Space by Euler Cromer method')
8
plt.show()
9

10
lable = label('time vs velocity ','time (sec)','velocity ( m/s ) ','--')
11
g3.newPlot(time, velocity,lable)
12

13
lable = label('oscillation amplitude','time (sec)','position (m)','--')
14
g3.newPlot(time, position  ,lable)

png

Here using the Euler cromer method we get the proper phase space trajectory and the energy of the system is conserved. Here the velocity and position shows an oscillatory variation with a phase difference of $π/2$ i.e. the two are in quadrature just as expected for a Simple harmonic oscillator.

Simulation of SHO: Euler Cromer.#

1
x_data = []
2
y_data = []
3
l = 10 # length of unstretched spring.
4

5
fig, ax = plt.subplots(figsize=(6,6))
6
ax = plt.axes(xlim=(-l, l*3),ylim=(-l*0.8, l*0.8 ))
7
line, = ax.plot([], [])
8
dot, = ax.plot([], [], 'sk', markersize = 20 )
9
plt.plot( np.zeros(10),np.linspace(0,2,10), '-k')
10
plt.plot( np.linspace(-10,30,10) ,np.zeros(10), '-b')
11

12
plt.xlabel('x (meter)')
13
plt.ylabel('y (meter)')
14
plt.title('Simulation using Euler cromer method')
15

16
def animation_frame(i):
17

18
  x_data = np.linspace(0, (l + position[i]), 50)
19
  temp = 0.5 * np.sin(x_data * (l*2*np.pi)/(l + position[i]))
20
  y_data = temp + 0.5 # just puting spring in place.
21

22
  line.set_xdata(x_data)
23
  line.set_ydata(y_data)
24

25
  dot.set_ydata([0.5])
26
  dot.set_xdata([l + position[i]])
27

28
  return line, dot
29

30
# simple plot is just the initial condition
31
# Animination is a video.
32

33
animation_vid = animation.FuncAnimation(fig, func = animation_frame, frames = np.arange(0, len(time)-1 , 1) , interval= dt*1000 )
34
HTML(animation_vid.to_html5_video())

png

Part B#

Use these functions to solve for the simple harmonic pendulum: $\frac{d^2 \theta}{dt^2} = − k \, sin \theta,$ with boundary conditions: at $t = 0$ , $\theta = \pi/4, \dot{\theta} = 0$ . Note that the initial amplitude is quite large, so it is not really a simple harmonic oscillator.

Plot $\theta$ Vs $t$ and the phase space trajectory.

Explore the system for small oscillations and explain any qualitative difference from the previous case. (For this choose suitable boundary conditions).

Solution using Euler Method.#

We are solving an ODE of the form $\frac{d\vec y}{dt} = \vec f$

For a simple harmonic pendulum, the equation of motion is second-order ODE:

$\frac{d^2 \theta}{dt^2} = − k \, sin \theta,$

we can vectorize this as

\frac{d}{dt} \begin{bmatrix} \theta \\ \omega \end{bmatrix} = \begin{bmatrix} \omega \\ - k sin(\theta) \end{bmatrix}

The initial condition is given: $\vec y_0 = \vec y(t=0) = \begin{bmatrix} \theta_0 \\ \omega_0 \end{bmatrix} = \begin{bmatrix} \pi/4 \\ 0.0 \end{bmatrix}$

1
# Functions used in Part B solution.
2
# the Functions are created using same principles from Part A.
3
def rhs(theta, omega, time):
4
  g = 9.8
5
  r = 1
6
  k = g/r
7
  return [omega, -k * np.sin(theta)]
8

9
def EulerODE(initCondn, tRange, rhs):
10

11
  dt = tRange[2]
12
  time = np.arange(tRange[0], tRange[1] + dt , dt)
13
  theta = np.zeros(len(time))
14
  omega = np.zeros(len(time))
15

16
  theta[0] = initCondn[0]
17
  omega[0] = initCondn[1]
18

19
  for i in range(0, len(time)-1):
20
    theta[i+1] = theta[i] + rhs(theta[i],omega[i], None )[0] * dt
21
    omega[i+1] = omega[i] + rhs(theta[i],omega[i], None )[1] * dt
22

23
  return theta, omega
24

25

26
def EulerCromerODE(initCondn, tRange, rhs):
27

28
  dt = tRange[2]
29
  time = np.arange(tRange[0], tRange[1] + dt , dt)
30
  theta = np.zeros(len(time))
31
  omega = np.zeros(len(time))
32

33
  theta[0] = initCondn[0]
34
  omega[0] = initCondn[1]
35

36
  for i in range(0, len(time)-1):
37
    omega[i+1] = omega[i] + rhs( theta[i], omega[i] , None )[1] * dt
38
    theta[i+1] = theta[i] + rhs( theta[i+1], omega[i+1] , None )[0] * dt
39
  return theta, omega

1
time_0 = 0
2
time_max = 5
3
time_step = 1/10 # 10th part of a second
4

5
theta_0 = np.pi/4
6
omega_0 = 0        # angular velocity
7

8
initCondn = [theta_0, omega_0]
9
tRange = [time_0, time_max, time_step]
10

11
dt = tRange[2]
12
time = np.arange(tRange[0], tRange[1] + dt , dt)
13

14
theta, omega = EulerODE(initCondn, tRange, rhs)
15

16
lable = label('trajectory in phase space','angle (radian)','angular velocity (radian/s)','--')
17
g3.newPlot(theta, omega,lable)
18

19
plt.plot(theta[0], omega[0],'*r')
20
#plt.legend()
21
plt.title('Phase Space by Euler method')
22
plt.show()
23

24
lable = label(' ','time (sec)','angular velocity (radian / sec)','--')
25
g3.newPlot(time, omega,lable)
26

27
lable = label(' ','time (sec)','theta (radian)','--')
28
g3.newPlot(time, theta  ,lable)
29
plt.show()
30

31
plt.polar(theta, np.ones(len(theta)), 'g.')
32
plt.title("all the angles taken by Pendulum")
33

34
#saving the data in dictionary form to use later for calculating energy
35
euelrSolution = {'theta' : theta, 'omega': omega}

png

We can see the angle taken by Pendulum are higher then the initial angle, which should not be posible. We get this error because of euler method , as it is just an approximation.

1
x_data = []
2
y_data = []
3
r = 1
4

5
fig, ax = plt.subplots(figsize=(6,6))
6
ax = plt.axes(xlim=(-r*1.1, r*1.1 ), ylim=(-r*1.1, r*1.1))
7
line, = ax.plot([], [])
8
dot, = ax.plot([], [],'ro', markersize = 10 )
9
plt.xlabel('x (meter)')
10
plt.ylabel('y (meter)')
11
plt.title("Simulation of Simple Pendulum: Euler Method")
12

13
def animation_frame(i):
14
  x_data = np.linspace(0, r * np.sin(theta[i]), 10)
15
  y_data = np.linspace(0, - r * np.cos(theta[i]) ,10 )
16
  line.set_xdata(x_data)
17
  line.set_ydata(y_data)
18

19
  dot.set_xdata([ r * np.sin(theta[i]) ])
20
  dot.set_ydata([- r * np.cos(theta[i])])
21

22
  return line, dot
23

24
animation_vid = animation.FuncAnimation(fig, func = animation_frame, frames = np.arange(0, len(time)-1 , 1) , interval= dt*1000 )
25
HTML(animation_vid.to_html5_video())

png

####Solution using Euler-Cromer algorithm

1
time_0 = 0
2
time_max = 5
3
time_step = 1/10 # 10th part of a second
4

5
theta_0 = np.pi/4
6
omega_0 = 0        # angular velocity
7

8
initCondn = [theta_0, omega_0]
9
tRange = [time_0, time_max, time_step]
10

11
dt = tRange[2]
12
time = np.arange(tRange[0], tRange[1] + dt , dt)
13

14
theta, omega = EulerCromerODE(initCondn, tRange, rhs)
15

16
lable = label('trajectory in phase space','angle (radian)','angular velocity (radian/s)','--')
17
g3.newPlot(theta, omega,lable)
18

19
plt.plot(theta[0], omega[0],'*r')
20
#plt.legend()
21
plt.title('Phase Space by Euler method')
22
plt.show()
23

24
lable = label(' ','time (sec)','angular velocity (radian / sec)','--')
25
g3.newPlot(time, omega,lable)
26

27
lable = label(' ','time (sec)','theta (radian)','--')
28
g3.newPlot(time, theta  ,lable)
29
plt.show()
30

31
plt.polar(theta, np.ones(len(theta)), 'g.')
32
plt.title("all the angles taken by Pendulum")
33

34
#saving the data in dictionary form to use later for calculating energy
35
euelrCromerSolution = {'theta' : theta, 'omega': omega}

png

1
x_data = []
2
y_data = []
3
r = 1 # meter
4
#length of pendulume rod
5

6
fig, ax = plt.subplots(figsize=(6,6))
7
ax = plt.axes(xlim=(-r*1.1, r*1.1 ), ylim=(-r*1.1, r*1.1))
8
line, = ax.plot(0, 0)
9
dot, = ax.plot(0, 0,'ro', markersize = 10 )
10
plt.xlabel('x (meter)')
11
plt.ylabel('y (meter)')
12
plt.title("Simulation of Simple Pendulum: Euler cromer Method")
13

14
def animation_frame(i):
15
  x_data = np.linspace(0, r * np.sin(theta[i]), 10)
16
  y_data = np.linspace(0, - r * np.cos(theta[i]) ,10 )
17
  line.set_xdata(x_data)
18
  line.set_ydata(y_data)
19

20
  dot.set_xdata([r * np.sin(theta[i])])
21
  dot.set_ydata([- r * np.cos(theta[i])])
22

23
  return line, dot
24

25
animation_vid = animation.FuncAnimation(fig, func=animation_frame, frames = np.arange(0, len(time)-1 , 1) , interval= dt*1000 )
26
HTML(animation_vid.to_html5_video())

png

Part C#

Plot the total energy as a function of time for both the methods. Notice any oscillations?

####Solution Energy

$\text{kinetic energy} =\frac{m}{2}(\dot r ^2 + r^2\dot \theta ^2)$

$\text{potential energy} = m g (h - r * cos(\theta) )$

1
# befor running this make sure that the initial condition is the same for both Methods above.
2

3
mass = 1 # kg
4
height = radius = 1 # meter
5
g = 9.8 #m/s
6

7
plt.figure(figsize=(7,5))
8

9
energy_euler = mass / 2 *( 0 +radius**2 * np.array(euelrSolution.get('omega'))**2 ) + mass * g * (height - radius*np.cos(euelrSolution.get('theta')))
10

11
energy_eulerCromer = mass / 2 *( 0 +radius**2 * np.array(euelrCromerSolution.get('omega'))**2 ) + mass * g * (height - radius* np.cos(euelrCromerSolution.get('theta')))
12

13
plt.plot(time, energy_euler, label='Euler')
14
plt.plot(time, energy_eulerCromer, label='Euler-Cromer')
15
plt.plot(time, np.zeros(len(time)), label='Zero Energy Level')
16

17
plt.legend()
18
plt.xlabel('time (sec)')
19
plt.ylabel('total energy (joule)')

1
Text(0, 0.5, 'total energy (joule)')

png

The total energy of a system executing simple harmonic motion remains constant i.e. remains conserved. The oscillations seen in the graph are due to the $r * cos(\theta)$ term in the potential energy. It is observed that increasing the radius gives more oscillations in the total energy graph.
For small oscillations, $\sin\theta$ can be approximated to $\theta$ and $\cos\theta$ as 1. Hence the equation of motion will be:

$\frac{d^2 \theta}{dt^2} = − k \,\theta$

##Problem 2: Exploring Liouville’s theorem with symplectic integrators

Part A#

The Euler and Euler-Cromer algorithms define a map from the $n^{th} \rightarrow (n+1)^{th}$ time steps. For a Simple Harmonic Oscillator (SHO), calculate the Jacobian of the transformation. For this, explain why Euler algorithm is not suitable for solving the Newton’s equations of motion?

####Solution For a Simple Harmonic Oscillator, Euler algorithm is not suitable for solving Newton’s equations of motion because along with inducing error in the solution due to cancellation of higher-order differential terms, it does not conserve the energy of the system. Euler algorithm is unstable since the trajectory produced by the algorithm spirals continually outwards away from the exact solution. This means that the Euler algorithm does not conserve energy. In the SHO problem we saw that in the Euler method, the object’s maximum displacement and maximum speed are steadily increasing, so its total energy must also be increasing.
According to Liouville’s Theorem, a conservative system (such as simple harmonic oscillator) must preserve the volume (or area) occupied by an ensemble of trajectories in phase space. The system evolves with time and the cluster of points move through phase space but the region containing the points remains constant. In the Euler method, this region in phase space grows in size with the passage of time. For SHO, the transformation in the Euler method is given by the set of below equations:

$x_{k+1} = x_k + v_{k}\Delta t$
$v_{k+1} = v_k -x_{k} \frac{k}{m} \Delta t$

If we treat the above equations as a two-dimensional map, then this map will preserve phase space volume if and only if the Jacobian of the map is equal to unity. The Jacobian of this transformation is given by:

$J = \begin{pmatrix} {\frac{\partial x_{n+1}}{\partial x_{n}}} & {\frac{\partial x_{n+1}}{\partial v_{n}}}\\ {\frac{\partial v_{n+1}}{\partial x_{n}}} & {\frac{\partial v_{n+1}}{\partial v_{n}}} \end{pmatrix}$ $= \begin{pmatrix} 1 & \Delta t \\ -\frac{k}{m} \Delta t & 1\end{pmatrix}$

$= 1+\frac{k}{m} \Delta t^2$

Thus in contradiction to Liouville’s Theorem, the Euler algorithm does not preserve phase space volume and grows in size with time. Hence it is not suitable for solving Newton’s equations of motion.
The transformation equations for SHO in the case of Euler-Cromer method are given by:

$x_{k+1} = x_k + v_{k+1}\Delta t$
$v_{k+1} = v_k -x_{k} \frac{k}{m} \Delta t$

The Jacobian of this transformation is given by:

$J = \begin{pmatrix} {\frac{\partial x_{n+1}}{\partial x_{n}}} & {\frac{\partial x_{n+1}}{\partial v_{n}}}\\ {\frac{\partial v_{n+1}}{\partial x_{n}}} & {\frac{\partial v_{n+1}}{\partial v_{n}}} \end{pmatrix}$ $= \begin{pmatrix} 1-\frac{k}{m} \Delta t^2 & \Delta t \\ -\frac{k}{m} \Delta t & 1\end{pmatrix}$ $= 1$

Hence Euler-Cromer method preserves the phase space volume and conserves energy for SHO system.

Part B#

For an SHO, consider a small square element of phase-space area at $t = 0$ whose corners are labelled A, B, C and D respectively. In this exercise, we would like to follow the area enclosed as the system evolves in time.

One way to do this is to evolve the four different initial conditions in time and plot the resultant polygon in the phase-space after regular intervals of time. For a conserved system, we expect this area conserved in time. For a symplectic solver (like Euler-Cromer with a global accuracy of $\mathcal{O} (\Delta t)$ ), this is true. However, this is not for the Euler algorithm.

Alternately, consider the trajectory of a cloud of points through phase space. Show that as the cloud stretches in one coordinate – $p$ say – it shrinks in the corresponding ${q}$ direction so that the product ${ \Delta p\,\Delta q}$ remains constant.

####Solution

1
def rhs(position, velocity, time):
2

3
  #k is the spring constant
4
  k = 9.8 #N·m−1
5
  m =1.   #mass kg
6
  return [velocity, -k /m * position]
7

8
#just name a class to store all functions
9
# this are the same definations frome above
10
#for using the class/object method to make code simple.
11

12
# as the following code is written with generalized functions
13
# the theta and omega can be and are  used alternatively for position and velocity and vice-versa
14

15
class Algorithms():
16
  def __init__(self, initCondn , tRange, rhs  ):
17
    self.initCondn = initCondn
18
    self.tRange = tRange
19
    self.rhs = rhs
20

21
  def EulerODE(self):
22

23
    dt = self.tRange[2]
24
    time = np.arange(self.tRange[0], self.tRange[1] + dt , dt)
25
    theta = np.zeros(len(time))
26
    omega = np.zeros(len(time))
27

28
    theta[0] = self.initCondn[0]
29
    omega[0] = self.initCondn[1]
30

31
    for i in range(0, len(time)-1):
32
      theta[i+1] = theta[i] + self.rhs(theta[i],omega[i], None )[0] * dt
33
      omega[i+1] = omega[i] + self.rhs(theta[i],omega[i], None )[1] * dt
34

35
    return theta, omega
36

37
  def EulerCromerODE(self):
38

39
    dt = self.tRange[2]
40
    time = np.arange(self.tRange[0], self.tRange[1] + dt , dt)
41
    theta = np.zeros(len(time))
42
    omega = np.zeros(len(time))
43

44
    theta[0] = self.initCondn[0]
45
    omega[0] = self.initCondn[1]
46

47
    for i in range(0, len(time)-1):
48
      omega[i+1] = omega[i] + self.rhs( theta[i], omega[i] , None )[1] * dt
49
      theta[i+1] = theta[i] + self.rhs( theta[i+1], omega[i+1] , None )[0] * dt
50
    return theta, omega

Euler Method#

1
time_0 = 0
2
time_max = 5
3
time_step = 1/10 # 10th part of a second
4

5
position = 0.3 # meters
6
velocity = 0.3 # m/s
7
dp = 0.1
8

9

10
#all the values are sored in dictionary.
11
phaseSpace = {}
12

13
#giving initial values to A,B,C,D
14
position_0 = position
15
velocity_0 = 0
16
AinitCondn = [position_0, velocity_0]
17

18
position_0 = 2*position
19
velocity_0 = 0
20
BinitCondn = [position_0, velocity_0]
21

22
position_0 = 2 * position
23
velocity_0 = 1 * velocity
24
CinitCondn = [position_0, velocity_0]
25

26

27
position_0 = 1*position
28
velocity_0 = 1*velocity
29
DinitCondn = [position_0, velocity_0]
30

31
tRange = [time_0, time_max, time_step]
32
dt = tRange[2]
33
time = np.arange(tRange[0], tRange[1] + dt , dt)
34

35
A = Algorithms(AinitCondn, tRange, rhs)
36
B = Algorithms(BinitCondn, tRange, rhs)
37
C = Algorithms(CinitCondn, tRange, rhs)
38
D = Algorithms(DinitCondn, tRange, rhs)
39

40

41
# all lines below are just to make the plots.
42
fig = plt.figure()
43

44
lable = label('A','position (meters)','velocity (m/s)','o')
45
phaseSpace['A_theta'],phaseSpace['A_omega'] = theta, omega = A.EulerODE()
46
g3.multiPlot(theta, omega,lable)
47

48
lable = label('B','position (meters)','velocity (m/s)','o')
49
phaseSpace['B_theta'],phaseSpace['B_omega'] = theta, omega = B.EulerODE()
50
g3.multiPlot(theta, omega,lable)
51

52
lable = label('C','position (meters)','velocity (m/s)','o')
53
phaseSpace['C_theta'],phaseSpace['C_omega'] = theta, omega = C.EulerODE()
54
g3.multiPlot(theta, omega,lable)
55

56
lable = label('D','position (meters)','velocity (m/s)','o')
57
phaseSpace['D_theta'],phaseSpace['D_omega'] = theta, omega = D.EulerODE()
58
g3.multiPlot(theta, omega,lable)
59

60
skip = 5
61
x = []
62
y = []
63
areaE = list()
64
fig = plt.figure()
65
#plt.axes(xlim=(-5,5), ylim=(-5,5))
66
plt.figsize=(10,10)
67

68
#This follow definations is to find are of any polygon.
69
#Only this folling 2 lines taken from online Q/A source (stackoverflow).
70
def PolyArea(x,y):
71
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
72
#---
73

74
for i in range(0,len(phaseSpace['A_theta']), skip):
75
  x =  [phaseSpace['A_theta'][i],phaseSpace['B_theta'][i],phaseSpace['C_theta'][i],phaseSpace['D_theta'][i],phaseSpace['A_theta'][i]]
76
  y = [phaseSpace['A_omega'][i],phaseSpace['B_omega'][i],phaseSpace['C_omega'][i],phaseSpace['D_omega'][i],phaseSpace['A_omega'][i]]
77
  areaE.append(PolyArea(x[:4],y[:4]))
78
  plt.plot(x,y,'--r')
79
  plt.plot(x,y,'*b')
80

81
plt.xlabel('position (meters)')
82
plt.ylabel('velocity (m/s)')
83
plt.title('area elements')
84
plt.show()
85

86
skip = 5
87
for i in range(0,len(phaseSpace['A_theta']), skip):
88
  x =  [phaseSpace['A_theta'][i],phaseSpace['B_theta'][i],phaseSpace['C_theta'][i],phaseSpace['D_theta'][i],phaseSpace['A_theta'][i]]
89
  y = [phaseSpace['A_omega'][i],phaseSpace['B_omega'][i],phaseSpace['C_omega'][i],phaseSpace['D_omega'][i],phaseSpace['A_omega'][i]]
90
  plt.fill(x[:4],y[:4])
91

92
lable = label('A','position (meters)','velocity (m/s)','o')
93
phaseSpace['A_theta'],phaseSpace['A_omega'] = theta, omega = A.EulerODE()
94
g3.multiPlot(theta, omega,lable)
95

96
lable = label('B','position (meters)','velocity (m/s)','o')
97
phaseSpace['B_theta'],phaseSpace['B_omega'] = theta, omega = B.EulerODE()
98
g3.multiPlot(theta, omega,lable)
99

100
lable = label('C','position (meters)','velocity (m/s)','o')
101
phaseSpace['C_theta'],phaseSpace['C_omega'] = theta, omega = C.EulerODE()
102
g3.multiPlot(theta, omega,lable)
103

104
lable = label('D','position (meters)','velocity (m/s)','o')
105
phaseSpace['D_theta'],phaseSpace['D_omega'] = theta, omega = D.EulerODE()
106
g3.multiPlot(theta, omega,lable)
107

108
plt.title('Phase Space with area element: euler method')
109
plt.show()

png

Euler Cromer Method#

1
lable = label('A','position (meters)','velocity (m/s)','o')
2
phaseSpace['A_theta'],phaseSpace['A_omega'] = theta, omega = A.EulerCromerODE()
3
g3.multiPlot(theta, omega,lable)
4

5
lable = label('B','position (meters)','velocity (m/s)','o')
6
phaseSpace['B_theta'],phaseSpace['B_omega'] = theta, omega = B.EulerCromerODE()
7
g3.multiPlot(theta, omega,lable)
8

9
lable = label('C','position (meters)','velocity (m/s)','o')
10
phaseSpace['C_theta'],phaseSpace['C_omega'] = theta, omega = C.EulerCromerODE()
11
g3.multiPlot(theta, omega,lable)
12

13
lable = label('D','position (meters)','velocity (m/s)','o')
14
phaseSpace['D_theta'],phaseSpace['D_omega'] = theta, omega = D.EulerCromerODE()
15
g3.multiPlot(theta, omega,lable)
16

17
skip = 5
18
x = []
19
y = []
20
areaEC = list()
21
fig = plt.figure()
22
plt.axes(xlim=(-2,2), ylim=(-2,2))
23
plt.figsize=(10,10)
24

25
def PolyArea(x,y):
26
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
27

28
for i in range(0,len(phaseSpace['A_theta']), skip):
29
  x =  [phaseSpace['A_theta'][i],phaseSpace['B_theta'][i],phaseSpace['C_theta'][i],phaseSpace['D_theta'][i],phaseSpace['A_theta'][i]]
30
  y = [phaseSpace['A_omega'][i],phaseSpace['B_omega'][i],phaseSpace['C_omega'][i],phaseSpace['D_omega'][i],phaseSpace['A_omega'][i]]
31
  areaEC.append(PolyArea(x[:4],y[:4]))
32
  plt.plot(x,y,'--r')
33
  plt.plot(x,y,'*b')
34
plt.xlabel('position (meters)')
35
plt.ylabel('velocity (m/s)')
36
plt.title('area elements')
37
plt.show()
38
for i in range(0,len(phaseSpace['A_theta']), skip):
39
  x =  [phaseSpace['A_theta'][i],phaseSpace['B_theta'][i],phaseSpace['C_theta'][i],phaseSpace['D_theta'][i],phaseSpace['A_theta'][i]]
40
  y = [phaseSpace['A_omega'][i],phaseSpace['B_omega'][i],phaseSpace['C_omega'][i],phaseSpace['D_omega'][i],phaseSpace['A_omega'][i]]
41
  plt.fill(x[:4],y[:4])
42

43
lable = label('A','position (meters)','velocity (m/s)','o')
44
phaseSpace['A_theta'],phaseSpace['A_omega'] = theta, omega = A.EulerCromerODE()
45
g3.multiPlot(theta, omega,lable)
46

47
lable = label('B','position (meters)','velocity (m/s)','o')
48
phaseSpace['B_theta'],phaseSpace['B_omega'] = theta, omega = B.EulerCromerODE()
49
g3.multiPlot(theta, omega,lable)
50

51
lable = label('C','position (meters)','velocity (m/s)','o')
52
phaseSpace['C_theta'],phaseSpace['C_omega'] = theta, omega = C.EulerCromerODE()
53
g3.multiPlot(theta, omega,lable)
54

55
lable = label('D','position (meters)','velocity (m/s)','o')
56
phaseSpace['D_theta'],phaseSpace['D_omega'] = theta, omega = D.EulerCromerODE()
57
g3.multiPlot(theta, omega,lable)
58

59
plt.title('Phase Space with area element : Eueler Cromer method')
60
plt.show()

png

1
#The area is not calculated for all elements just for the coloured ones above.
2
plt.plot(range(len(areaEC)),areaEC,'-g',label ="Euler Cromer , max = {0} m^2".format(round(max(areaEC),2)))
3
plt.plot(range(len(areaE)),areaE,'--r',label ="Euler, max = {0} m^2".format(round(max(areaE),2)))
4
plt.legend()
5
plt.xlabel('area element taken as time passes')
6
plt.ylabel('area (m^2)')
7
plt.show()

png

The area elements generated by choosing four initial conditions and evolving (numerically solving) in time results in an increasing area when the Euler method is used. and when we use the Euler cromer method the area is conserved.

##Problem 3: Linearly-damped harmonic oscillator

In the real world, all classical harmonic oscillators have some form of damping that converts their mechanical energy into heat and eventually brings the system to rest. The simplest way to represent such an effect is to add a linear term $−by$ to the simple harmonic oscillator corresponding to a frictional force that is proportional to the velocity, \begin{eqnarray} \dot{x} &=& v, \ \dot{v} &=& -x. \end{eqnarray}

The resulting system is called a damped harmonic oscillator, and b (if positive) is the damping constant. It describes the exponential rate at which orbits spiral into the origin of the phase-space at $(0,0)$ and is related to the $Q$ (quality factor) of the oscillator by $b = 1/Q$ . The $Q$ of an oscillator is the number of radians of oscillation required for the energy to decay to $1/e$ of its original value. A system with $Q = 1/2$ (or $b = 2$ ) is critically damped, and smaller values of $Q$ (or $b > 2$ ) do not oscillate, but they rapidly approach the origin.

It seems that it should be possible to eliminate the $b$ in above by a linear rescaling of $x$ , $v$ , and $t$ , but that cannot be done as simple algebra shows.

This system is special in that it has two distinct time-scales, the damping rate and the frequency of oscillation, and the parameter $b$ controls their ratio.

Write a code to plot $x(t)$ , $v(t)$ and phase-space trajectories for different values of the quality factor (demonstrating underdamped, critically damped and overdamped behaviour of this oscillator). Use both the Euler and Euler-Cromer algorithms. Write down an expression for the determinant of the Jacobian $J$ that maps $(x_n, v_n) \rightarrow (x_{n+1}, v_{n+1})$ in both the cases.

The damped harmonic oscillator is an example of a dissipative system. For $b>0$ is it time-reversible? Explain your answer.

SOLUTION:#

We are solving $\frac{d^2x}{dt^2}= −\frac{k}{m}x -bv$
We can vectorize it in the form

$\frac{d \vec y}{dt} = \vec f$

Explicitly:

\frac{d}{dt} \begin{bmatrix} x \\ v \end{bmatrix} = \begin{bmatrix} v \\ -\frac{k}{m}x - bv \end{bmatrix}

For Damped SHO, the tranformation in Euler method is given by the set of below equations:

$x_{k+1} = x_k + v_{k}\Delta t$
$v_{k+1} = v_k -x_{k} \frac{k}{m} \Delta t -bv_k \Delta t$

The Jacobian of this transformation is given by:

$= 1- b\Delta t+\frac{k}{m} \Delta t^2$

The transformation equations for SHO in the case of Euler-Cromer method are given by:

$x_{k+1} = x_k + v_{k+1}\Delta t$
$v_{k+1} = v_k -x_{k} \frac{k}{m} \Delta t-bv_k \Delta t$

The Jacobian of this transformation is given by:

$J = \begin{pmatrix} {\frac{\partial x_{n+1}}{\partial x_{n}}} & {\frac{\partial x_{n+1}}{\partial v_{n}}}\\ {\frac{\partial v_{n+1}}{\partial x_{n}}} & {\frac{\partial v_{n+1}}{\partial v_{n}}} \end{pmatrix}$ $= \begin{pmatrix} 1-\frac{k}{m} \Delta t^2 & \Delta t-b\Delta t^2 \\ -\frac{k}{m} \Delta t & 1-b\Delta t\end{pmatrix}$

$= 1-b\Delta t$

Hence we can say that Euler-Cromer does not preserve the energy for damped oscillator system due to which we need to use other numerical methods to get a more accurate result.

The solution is an exponentially decaying with time. Reversing time will make it an exponentially growing function.

Hence we can say that it is not time reversible.

1
# This is the defination for SHO with drag.
2
def a(xk,vk,t,):
3
  mass = 1 #kg
4
  k = 1 # spring constant
5
  b = 1/quality_factor # ve are varing value of b
6
  return -1/mass* k* xk - b * vk
7

8
# just the same definitions above again to make each solution self compiles
9

10
def rhs(xk,vk,tk):
11

12
  return [vk, a(xk,vk,tk)]
13

14

15
def EulerODE(initCondn, tRange, rhs):
16

17
  dt = tRange[2]
18
  time = np.arange(tRange[0], tRange[1] + dt , dt)
19

20
  position = np.zeros(len(time))
21
  velocity = np.zeros(len(time))
22

23
  position[0] = initCondn[0]
24
  velocity[0] = initCondn[1]
25

26
  for i in range(0, len(time)-1):
27
    position[i+1] = position[i] + rhs(position[i],velocity[i], None )[0] * dt
28
    velocity[i+1] = velocity[i] + rhs(position[i],velocity[i], None )[1] * dt
29

30
  return position, velocity
31

32

33
def EulerCromerODE(initCondn, tRange, rhs):
34

35
  dt = tRange[2]
36
  time = np.arange(tRange[0], tRange[1] + dt , dt)
37
  position = np.zeros(len(time))
38
  velocity = np.zeros(len(time))
39

40
  position[0] = initCondn[0]
41
  velocity[0] = initCondn[1]
42

43
  for i in range(0, len(time)-1):
44
    velocity[i+1] = velocity[i] + rhs(position[i],velocity[i], None )[1] * dt
45
    position[i+1] = position[i] + rhs(position[i+1],velocity[i+1], None )[0] * dt
46

47
  return position, velocity

1
time_0 = 0
2
time_max = 50
3
time_step = 1/10 # 10th part of a second
4

5
position_0 = 5 # in meters difference from rest position
6
velocity_0 = 0        # angular velocity
7

8
initCondn = [position_0, velocity_0]
9
tRange = [time_0, time_max, time_step]
10

11
dt = tRange[2]
12
time = np.arange(tRange[0], tRange[1] + dt , dt)
13

14
quality_factor = 1
15
position1, velocity1 = EulerODE(initCondn, tRange, rhs)
16

17
quality_factor = 2
18
position2, velocity2 = EulerODE(initCondn, tRange, rhs)
19

20
quality_factor = 3
21
position3, velocity3 = EulerODE(initCondn, tRange, rhs)
22

23
lable = label('Q = 1','angle(rad)','Angular velocity (rad/s)','--')
24
g3.newPlot(position1, velocity1,lable)
25
lable = label('Q = 2','angle(rad)','Angular velocity (rad/s)','--')
26
g3.multiPlot(position2, velocity2,lable)
27
lable = label('Q = 3','angle(rad)','Angular velocity (rad/s)','--')
28
g3.multiPlot(position3, velocity3,lable)
29
plt.plot(position_0,velocity_0,'*r')
30

31
lable = label('Q = 1','time(s)','Angular velocity (rad/s)','--')
32
g3.newPlot(time, velocity1,lable)
33
lable = label('Q = 2','time(s)','Angular velocity (rad/s)','--')
34
g3.multiPlot(time, velocity2,lable)
35
lable = label('Q = 3','time(s)','Angular velocity (rad/s)','--')
36
g3.multiPlot(time, velocity3,lable)
37

38
lable = label('Q = 1','time(s)','angle(rad)','--')
39
g3.newPlot(time, position1,lable)
40
lable = label('Q = 2','time(s)','angle(rad)','--')
41
g3.multiPlot(time, position2,lable)
42
lable = label('Q = 3','time(s)','angle(rad)','--')
43
g3.multiPlot(time, position3,lable)

png

#####Interactive solutions

1
def a(xk,vk,t):
2
  mass = mass_i.value #kg
3
  k = spring_const_i.value # spring constant
4
  b = 1/quality_factor_i.value
5
  return -1/mass* k* xk - b * vk
6

7

8
time_0 = 0
9
time_max = 20
10
time_step = 1/10 # 10th part of a second
11
tRange = [time_0, time_max, time_step]
12
dt = tRange[2]
13
time = np.arange(tRange[0], tRange[1] + dt , dt)
14

15
#--
16
# Defining the variables.
17

18
#definations made for interactive plot.
19
def update_plot(mass,spring_const,quality_factor,position_0):
20

21
  initCondn = [position_0, velocity_0]
22
  position, velocity = EulerODE(initCondn, tRange, rhs)
23
  position_EC, velocity_EC = EulerCromerODE(initCondn, tRange, rhs)
24

25
  L = label('for EulerODE Q = {0}'.format(quality_factor),'Position (meters)','velocity (meters/sec)','-r')
26
  g3.multiPlot(position, velocity, L )
27
  plt.plot(position[0], velocity[0], '*r' )
28

29
  L = label('for EulerCromerODE Q = {0}'.format(quality_factor),'Position (meters)','velocity (meters/sec)','-k')
30
  g3.multiPlot(position_EC, velocity_EC, L )
31

32
#making slider
33
mass_i = widgets.FloatSlider(min=0, max=10, value=1,step = 0.5 , description = 'Mass (kg) : ')
34
spring_const_i = widgets.FloatSlider(min=0, max=10, value=1,step = 0.5 , description = 'Spring Constant : ')
35
quality_factor_i = widgets.FloatSlider(min=0.1, max=3, value=1,step = 0.2 , description = 'Quality factor : ')
36
position_0_i = widgets.FloatSlider(min=0, max=10, value=1,step = 0.5 , description = 'Initial position from rest (meters) : ')
37
#Starting the intractive session
38
widgets.interactive(update_plot,mass = mass_i, spring_const =spring_const_i ,quality_factor = quality_factor_i,position_0=position_0_i)

1
interactive(children=(FloatSlider(value=1.0, description='Mass (kg) : ', max=10.0, step=0.5), FloatSlider(valu…

The End#

To Open Lab 3 https://github.com/ashwin-r-k/Computational-Physics_3labs/blob/main/LAB_3.ipynb