Solving second order Ordinary Differential Equations

AIM: Solving second order ODEs

OBJECTIVES OF THE PROJECT:

• objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping

THEORITICAL BACKGROUND

• The figure at the right shows an idealized pendulum, with a "massless" string or rod of lengthL and a bob of mass m.
• The open circle shows the rest position of the bob. When the bob is moved from its rest position and let go, it swings back and forth.
• The primary forces acting on the bob are the gravitational force that makes it move in the first place and the force exerted by the string to keep it moving along a circular path.
• In addition, there may be a dampingforce from friction at the pivot or air resistance or both.
• We will construct a model to describe how the angle thetaof the pendulum varies as a function of time t.

For the system shown :

According to second law of motion

`F=m×a`

`τ=Iα`

`⇒-mgsinθ=mL^2 (d^2 θ)/dt^2`

Assuming damping force is directly proportional to velocity & including it in our equation:

`-b (dθ)/dt -mgsinθ=mL^2 (d^2 θ)/dt^2`

`⇒mL^2 (d^2 θ)/dt^2 +b dθ/dt+mgsinθ=0`

Assuming `θ=θ_1` so , `(dθ)/dt=(dθ_1)/dt=θ_2`

`(d^2 θ)/dt^2 +b/m dθ/dt+g/L sinθ=0`

`(dθ_2)/dt+b/m θ_2+g/L sinθ_1=0`

On resolving this and converting it into matrix form we get :

`d/dt [(theta_1),(theta_2)]=[[θ_2,0],[-b/m θ_2,(-g)/L sinθ_1 ]]`

This is how we can reduce a 2^nd order differential equation into a 1^st order differential equation.

Workflow:

Analytical Solution

MATLAB CODE:

clear all
close all
clc

%INPUTS
l=1; % length of the pendulum in m
m=1; %mass of the ball in kg
b=0.05; %damping coefficient
g=9.81; %gravity in m/s2

%Initial Conditions

theta_0=0; %Angular displacement
omega_0=3; %angular velocity

%time_points
t_span=linspace(0,20,500); %Observing the motion for 20 Secs

%ODE solver
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,[0 3]);

ode_func code:

function [dtheta_dt]=ode_func(t,theta,b,g,l,m)

theta1=theta(1)
theta2=theta(2)

dtheta1_dt=theta2
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt]


end

figure(1)
A=plot(t,results(:,1)); %plotting angular displcement  
legend([A],'Angular Displacement');
hold on
B=plot(t,results(:,2))  %plotting angular velocity 
legend([B],'Angular Displacement');
hold off

disp=results(:,1); %Storing the angular displacement in array
ct=1               % Initializing the counter 

for i=1:length(t_span)
 x0=0;                              %Hinge Points 
 y0=0;                              %Hinge Points 
 x1=l*sin(disp(i));                 %Bob 
 y1=-l*cos(disp(i));
 figure(2)
 plot([-2 2],[0 0],'linewidth',5)
 hold on
 
 plot([x0  x1],[y0 y1],'-o','markersize',3,'linewidth',3,'color','y');
 plot(x1,y1,'-o','markersize',15,'markerfacecolor','r');
 axis([-1.5 1.5 -1.5 1.5]) %Defining Axis
 pause(0.03)
 hold off
 M(ct)=getframe(gcf)
 
 ct=ct+1
end

%To create a movie from the array
movie(M); %Defining the movie storing array M
videofile=VideoWriter('Pendulum.avi','Uncompressed AVI'); %Creating videofile object
open(videofile); %opening the videofile object
writeVideo(videofile,M); %writing the frames into a video
close(videofile); %Closing of Video object

Code Explanation:

1. Define the Input parameters as g = gravity in m/s2,,L = length of the pendulum in m,m = mass of the ball in kg, b=damping coefficient.
2. To use a ode45 solver, create a ode_func (Anonymous function) which will return a column matrix for the solver with the reduced 1^st order equations.
3. Provide time limits and initial conditions according to the ode45 syntax
4. The results will contain two separate columns, access the result array to plot the angular velocity and angular displacement.
5. Using FOR loop capture the variation of time and plot the location of pendulum with co-ordinates varying as a function of angular displacement
6. Generate the animation using movie command.

OBSERVATION:

The damped harmonic motion of the pendulum can be observed from the graph where we can observe both the velocity and the instantaneous displacement dying off due to the damping effect.

YouTube Animation

CONCLUSION:

Solution for a 2^nd order differential equation for a pendulum can be solved in MATALAB using ode45 solver by reducing the 2^nd order equation into a 1^st order equation.