Week 3 - Solving second order ODEs

AIM :

To write a program that solves second order differntial equation and to animate the motion of a simple pendulum based on the second order differential equation.

THEORY :

The second order differntial equation is as follows ;

d^2(theta)/dt^2 + (b/m)dtheta/dt + (g/l)sin(theta) = 0 -----------------------------------------(1)

In this equation ;

b= damping coefficient 

m = mass 

= accelaration due to gravity in m/s^2`

l= length of the string in meter

To solve this equation  , the second order differntial equation is broken down to first order differntial equations as follows ;

Let , 

theta= theta1-----------------------------------------------------------------------------------------------(2)

differntiating both sides with respect to time ;

d(theta)/dt = d(theta1)/dt

and this is assigned to another theta value theta2 

theta2 = d(theta1)/dt = d(theta)/dt-------------------------------------------------------------------------------(3)`

Differntiating again with respect to t :

d^2(theta)/dt^2 = d/dt(d(theta1)/dt)= d^2(theta1)/dt^2 = d(theta2)/dt -----------------------------------------(4)

 Putting the values of equation 2 , equation 3 and equation 4 in equation 1 we get ;

d(theta2)/dt + (b/m)d(theta2)/dt +(g/l)sin(theta1) =0 -----------------------------------------------------------(5)

Now we write the two equations of theta1 and theta2 in an array form:

matrix[[theta1],[theta2]] = matrices[[theta2],[ -(b/m)theta2 - (g/l)sin(theta1)]]

The left hand side is theta array , and this is the solver matrix

CODE EXPLANATION :

The different inputs are provided that are provided in the problem statement . 'b' is the damping coefficient, 'g' is the accelaration due to gravity , 'l' is the length of the string .'m' is the mass of the pendulum bob.

As it is a second order differntial equation , it should have two initial conditions ; one for intial displacement and one for initial velocity .Initial displacement is 0 and initial velocity is 3 rad/s.

As we integrate the equation , it should have a time span . the time points between which the equation is integrated. Here it is given as between 0 and 20 seconds. The value is calculated at 500 points.

The ODE function is called now. ode45 is the ODE solver that matlab provides. ode_func is the ODE function that is a separate matlab code written. The ode45 solver takes into account the t and theta values. These are the dependent and indeendent variables. After that the ode_func and the various arguments are provided. This is stored as an array form .

The ODE function is written as 

The first value of theta is theta1 and the second value of theta is theta2. dtheta1_dt is our first ode and dtheta2_dt is our second ode. Both the values are clubbed inside dtheta_dt  and it is stored as two rows separated by a semicolon . This function is stored as ode_func file name and called in the main program.

The plotting is done using the matplotlib.pyplot . the displacement and velocity is plotted . 

A counter variable is assigned as ct=1; A for loop is created for tracing the pendulum for every theta value .the base of the pendulum is created by assigning the coordinates and the bob is made by markers command . A movie is created using the matlab command MOVIE. a avi file is created in the name seconorder . the animation is uploaded in youtube. the youtube link is provided in the report.

MATLAB CODE:

The Matlab code for the ode function is 

function [dtheta_dt] = ODEfunction(t_span,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

The matlab code for the program is :

clear all 
close all 
clc 

%inputs 
b = 0.05;%damping coefficient
g = 9.81;%accelaration due to gravity
l = 1;%length of pendulum string
m = 0.1;%mass of bob

%initial conditions
theta_0 = [0;3];%initial displacement is 0 and initial velocity is 3
%time span
t_span = linspace(0,20,500);%between 0-20s 
%ODE solver 
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

% results
figure(1)
plot(t,results(:,1),'linewidth',2,'color','r')%plotting the displacement
hold on
plot(t,results(:,2),'linewidth',2,'color','b')%plotting the velocity
hold off
legend('displacement','velocity')
ylabel('plot')
xlabel('time')

ct = 1;%assigning a counter variable
h = results(:,1);
%a for loop to create different images for different theta values
for i = 1:length(h)
    theta_h = h(i);
    x0 = 0;
    y0 = 0;
    x1 = l*sin(theta_h);
    y1 =-l*cos(theta_h);
    plot([x0 x1],[y0 y1],'linewidth',2,'color','r')
    hold on
    plot(x1,y1,'o','markers',20,'markerfacecolor','g')
    hold on
    plot([-1 1],[0 0],'linewidth',3,'color','black')%creating the base of pendulum
    grid on
    hold off
    axis([-3,3,-2,2])
    pause(0.005)
    M(ct)= getframe(gcf);
    ct = ct+1;
end
%creating a movie
movie(M)
videofile = VideoWriter('secondorder.avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

OUTPUT :

The plot obtained is pasted below

 The youtube link for the animation is pasted below :

https://youtu.be/rMNp2d0504k

CONCLUSION:

The animation of the pendulum can be simulated between differnt time points and for various values damping coefficients .

The initial displacement and velocities can also be altered .

This finds use in multibody dynamics of automotive applications.