Week 3 - Solving second order ODEs

AIM : 

 to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping and to write a program to simulate the pendulam motion.

THEORY : 

what is ODE ?

It is a type of differential equation that involves a function of a single variable and its derivatives. ODEs are commonly used to describe relationships between variables in various scientific, engineering, and mathematical contexts. They play a fundamental role in modeling dynamic systems and processes, as they describe how quantities change with respect to time or another independent variable.

why it is useful ?

• predicting the future
• understand the nature
• computer simulation
• modelling tools
• modelling dynamic systems

simple pendulum :

the simple pendulum is another mechanical system that moves in an oscillatory motion.  It consists of a point mass ‘m’ suspended by means of light inextensible string of length L from a fixed support. The motion occurs in a vertical plane and is driven by a gravitational force. The forces which are acting on the mass are shown. The tangential component of the gravitational force, mg sin θ, always acts towards the mean position θ = 0 opposite to the displacement, restoring force acting tangent to the arc.

 Basic format is derived from F = ma, the equation of the simple pendulum will be

total force applied to the body (rotational force) = motion of the body (rotaional motion)

torque = moment of inertia * angular acceleration 

τ = I * α

where Torque is a measure of the force that can cause an object to rotate about an axis. 

τ = -mg L sinθ

m = mass of the object

g = acceleration due to gravity

L = lenght of the string in the pendulum

sinθ = It is a tangential component. This component is in the direction perpendicular to the string and is responsible for providing the restoring force for the pendulum's motion.

The moment of inertia (I) is a measure of an object's resistance to changes in its rotational motion. It depends on both the mass distribution of the object and the axis of rotation. For a simple pendulum, the pivot point is the axis of rotation, and the moment of inertia of the pendulum bob is given by

I = mL^2

and the angular acceleration is given by the second order derviative, d2θ/dt2. where is equation is differentiated by the first order ODE. it represent the angular velocity 

hence the overall above equation will be,

-mg L sinθ = mL^2 * d2θ/dt2

g sinθ = L * d2∂/dt2

d2θ/dt2 = - g sinθ / L

the equation shown above is no damping, (air resistance , or any other friction). 

Now let's look at the case where the damping gets involved

where 

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

for writing a simple pendulum program in python, we must convert second order ODE to the first order ODE, 

θ = θ1

dθ / dt = dθ1 / dt = θ2 ...................................equation 1

d2θ/dt2 = d2θ1/dt2 = d/dt (θ2)

dθ2/dt = -b/m  θ2 - g/L * sinθ1....................equation 2 

from equation 1 and 2 , we get

d/dt(θ1,θ2) = [θ2, -b/m  θ2 - g/L * sinθ1] 

The position of the bob can be determined in cartesian co-ordiantes as

x = l sinθ

y = -l cos θ

PYTHON CODE FOR THE ABOVE DERIVATION OF SIMPLE PENDULUM : 

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from PIL import Image

def model(theta, t, b, m, g, l):
           theta1 = theta [0]
           theta2 = theta [1]
           dtheta1_dt = theta2
           dtheta2_dt = -(b/m)*theta2 - (g/L)*math.sin(theta1)
           dtheta_dt = (dtheta1_dt,dtheta2_dt)
           return dtheta_dt

#lenght of the pendulum in m
L = 1

#gravity in m/s2
g = 9.8

#mass of the ball in kg
m = 1

#damping co-efficient 
b = 0.05

#initial condition
theta_0 = [0,3]

#time interval
t = np.linspace(0,30,125)

#ODE equation
theta = odeint(model,theta_0,t,args=(b,m,g,L))

#plotting for velocity and displacemet curve
plt.plot(t,theta[:,0],'-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'--',label=r'$\frac{d\theta2}{dt} = -\frac{b}{m}\theta_2-\frac{g}{L}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

# Create the images for the animation
image_files = []
ct = 1
for i in theta[:, 0]:
    x1 = 0
    y1 = 0
    x2 = L * math.sin(i)
    y2 = -L * math.cos(i)  # Negative sign to flip the y-axis for plotting
    file_name = f"{ct:03d}.png"  # Added leading zeros for filename formatting
    ct = ct + 1
    image_files.append(file_name)

    plt.figure()
    plt.plot([-2, 2], [0, 0])
    plt.plot([x1, x2], [y1, y2])
    plt.plot(x2, y2, 'o', markersize=15)
    plt.xlim(-2, 2)
    plt.ylim(-2, 2)
    plt.savefig(file_name)
    plt.close()

    # Add a brief pause before saving each frame
    plt.pause(0.1)  # Adjust the pause duration as needed

# Convert the images to a GIF
duration = 200  # Duration in milliseconds for each frame
loop = 0        # 0 means loop forever

# Open and add each image to the image list
image_list = []
for image_file in image_files:
    image = Image.open(image_file)
    image_list.append(image)

# Save the GIF
image_list[0].save("simple_pendulum.gif", save_all=True, append_images=image_list[1:], duration=duration, loop=loop)

CODE EXPLANATION : 

• import math
  import matplotlib.pyplot as plt
  import numpy as np
  from scipy.integrate import odeint
  from PIL import Image
  

Here, the code is importing necessary libraries and modules. math is used for mathematical operations, matplotlib.pyplot for plotting, numpy for numerical operations, scipy.integrate.odeint for solving ordinary differential equations, and PIL.Image for working with images.

• def model(theta, t, b, m, g, l):
      theta1 = theta[0]
      theta2 = theta[1]
      dtheta1_dt = theta2
      dtheta2_dt = -(b/m)*theta2 - (g/L)*math.sin(theta1)
      dtheta_dt = (dtheta1_dt, dtheta2_dt)
      return dtheta_dt

This function defines the system of differential equations that describe the motion of the damped pendulum. It takes the current state theta, time t, and parameters b, m, g, and l as arguments. It calculates the derivatives of angles (dtheta1_dt and dtheta2_dt) and returns them as a tuple dtheta_dt.

dtheta1_dt represents the rate of change of the first angle (theta1) with respect to time. It's the angular velocity of the pendulum bob.
dtheta2_dt represents the rate of change of the second angle (theta2) with respect to time. It's the angular acceleration of the pendulum bob.

• L = 1
  g = 9.8
  m = 1
  b = 0.05
  theta_0 = [0, 3]
  t = np.linspace(0, 30, 125)

Here, the parameters for the pendulum are set: L is the length of the pendulum, g is the acceleration due to gravity, m is the mass of the pendulum bob, b is the damping coefficient, theta_0 are the initial conditions (initial angle and angular velocity), and t is an array of time points.

• theta = odeint(model, theta_0, t, args=(b, m, g, L))

Using odeint from SciPy, this line numerically integrates the system of differential equations over the specified time points t using the defined model function and the given parameters.

• plt.plot(t, theta[:, 0], '-', label=r'$\frac{d\theta_1}{dt}=\theta2$')

plt.plot(t, theta[:, 0], '-'): This part of the line creates a plot of the angular displacement (theta[:, 0]) against time (t). The '-' specifies that a solid line should be used for the plot.
label=r'$\frac{d\theta_1}{dt}=\theta2$': This part adds a label to the plot. The label is a LaTeX-formatted string that represents the equation dθ1 / dt = θ2 ​which indicates the angular velocity.

• plt.plot(t, theta[:, 1], '--', label=r'$\frac{d\theta2}{dt} = -\frac{b}{m}\theta_2-\frac{g}{L}\sin\theta_1$')

This line is similar to the previous one, but it plots the angular velocity (theta[:, 1]) against time (t) using a dashed line ('--').
The label is a LaTeX-formatted string that represents the equation dθ2/dt = -b/m  θ2 - g/L * sinθ1​, which indicates the angular acceleration with damping and gravitational terms.

These lines plot the angular displacement and velocity over time using the results obtained from odeint. The plt.plot function is used to create the plots, and plt.show() displays the plot.

• plt.legend(loc='best')

plt.legend() adds a legend to the plot, which explains the different lines in the plot.
loc='best' specifies that the legend should be placed at the location that obstructs the plot the least.

image_files = []
ct = 1
for i in theta[:, 0]:
    x1 = 0
    y1 = 0
    x2 = L * math.sin(i)
    y2 = -L * math.cos(i)
    file_name = f"{ct:03d}.png"
    ct = ct + 1
    image_files.append(file_name)

    plt.figure()
    plt.plot([-2, 2], [0, 0])
    plt.plot([x1, x2], [y1, y2])
    plt.plot(x2, y2, 'o', markersize=15)
    plt.xlim(-2, 2)
    plt.ylim(-2, 2)
    plt.savefig(file_name)
    plt.close()

    plt.pause(0.1)

This loop creates a series of images to visualize the motion of the pendulum over time. It calculates the coordinates of the pendulum bob and saves a plot as an image for each time step. The images are saved with a numbered filename (e.g., "001.png", "002.png", etc.).

• image_list = []
  for image_file in image_files:
      image = Image.open(image_file)
      image_list.append(image)
      image_list[0].save("simple_pendulum.gif", save_all=True, append_images=image_list[1:], duration=200, loop=0)

This part of the code opens the saved images, compiles them into a GIF animation, and saves the animation as "simple_pendulum.gif". The duration parameter specifies the time duration of each frame in milliseconds, and the loop parameter specifies whether the animation should loop indefinitely (0 means loop forever).

Overall, this code calculates and visualizes the motion of a damped simple pendulum using numerical integration and creates a GIF animation to show the pendulum's oscillation over time.

ERROR WHILE PROGRAMMING : 

• the pendulum is graphed in upwards which is meaning less , this is due to the y axis assigned in the positive co-ordinates

• this is rectified using adding negative symbol in y2 and thus gives the appropriate results.

OUTPUT GRAPH :

• the blue line indicates the velocity curve which starts from zero
• the orange line indicates the acceleration curve which starts from 3 m/s2
• as its a damped controlled system, there will be a gradually decrease of the oscillation curve.
• the damped controlled system is acts by means of air resistance or friction
• this attains a equilibrium position after certain period of time

ANIMATION OF SIMPLE PENDULUM :

• this graph is developed byy using the below program line 
•  for i in theta[:, 0]:
      x1 = 0
      y1 = 0
      x2 = L * math.sin(i)
      y2 = -L * math.cos(i)
• this indicates that x1 and y1 will be in 0 and corresponding x2 and y2 are performed by the trignometry function, this used to be calculate the one position.
• but the for loop determines the runs for 125 times.
• thus the result is iterated, stored as separate file and it has been saved as the animation

CONCLUSION : 

Hence ODE of simple pendulum is solved using python and respective plots were plotted with respect to time.