Executive Programs

Workshops

Projects

Blogs

Careers

Placements

Student Reviews

For Business

Academic Training

Informative Articles

Find Jobs

We are Hiring!

All Courses

Choose a category

All Courses

Executive Programs

Workshops

For Business

Success Stories

Placements

Student Reviews

Projects

Blogs

Academic Training

Find Jobs

Informative Articles

We're Hiring!

+91 9342691281 Log in

Week 4.1 - Solving second order ODEs

In this project, I will be writing code in MATLAB to solve the motion of a simple pendulum. A simple pendulum motion's depends on Newton's Second Law. The equation which governs the motion of a simple pendulum is (with damping) $\frac{d^{2} θ}{{d t}^{2}} + \frac{b}{m} \frac{d θ}{d t} + \frac{g}{L} \sin θ = 0$ $(d^2theta)/(dt^2) + b/m(d theta)/(dt) + g/L sin theta = 0$ Where, $θ$ $theta$ is the angular displacement…

MATLAB

Dushyanth Srinivasan
updated on 23 Aug 2022

In this project, I will be writing code in MATLAB to solve the motion of a simple pendulum.

A simple pendulum motion's depends on Newton's Second Law. The equation which governs the motion of a simple pendulum is (with damping)

$\frac{d^{2} θ}{{d t}^{2}} + \frac{b}{m} \frac{d θ}{d t} + \frac{g}{L} \sin θ = 0$ $(d^2theta)/(dt^2) + b/m(d theta)/(dt) + g/L sin theta = 0$

Where, $θ$ $theta$ is the angular displacement of the pendulum ( $r a d$ $rad$ ), $b$ $b$ is the damping coefficient ( $N s / m$ $Ns//m$ ), $m$ $m$ is the mass of the pendulum ( $k g$ $kg$ ), $g$ $g$ is the gravitational acceleration ( $m / s^{2}$ $m//s^2$ ) and $L$ $L$ is the length of the pendulum's string ( $m$ $m$ ).

This equation is a second order differential equation, these equations cannot be directly solved in MATLAB and hence will be subjected to a few changes as follows:

$L e t \frac{d θ_{1}}{d t} = θ_{2} and θ = θ_{1}$ $Let (d theta_1)/(dt) = theta_2 and theta = theta_1$

Substituting $θ = θ_{1}$ $theta = theta_1$ in $\frac{d^{2} θ}{{d t}^{2}} + \frac{b}{m} \frac{d θ}{d t} + \frac{g}{L} \sin θ = 0$ $(d^2theta)/(dt^2) + b/m(d theta)/(dt) + g/L sin theta = 0$

$\frac{d^{2} θ_{1}}{{d t}^{2}} + \frac{b}{m} \frac{d θ_{1}}{d t} + \frac{g}{L} {\sin θ}_{1} = 0$ $(d^2theta_1)/(dt^2) + b/m(d theta_1)/(dt) + g/L sin theta_1 = 0$

$\frac{d (\frac{d θ_{1}}{d t})}{d t} + \frac{b}{m} \frac{d θ_{1}}{d t} + \frac{g}{L} {\sin θ}_{1} = 0$ $(d((d theta_1)/(dt)))/(dt) + b/m (d theta_1)/(dt) + g/L sin theta_1 =0$

Substituting $\frac{d θ_{1}}{d t} = θ_{2}$ $(d theta_1)/(dt) = theta_2$ ,

$\frac{d θ_{2}}{d t} + \frac{b}{m} θ_{2} + \frac{g}{L} {\sin θ}_{1} = 0$ $(d theta_2)/(dt) + b/m theta_2 + g/L sin theta_1 =0$

Rearranging,

$\frac{d θ_{2}}{d t} = - \frac{b}{m} θ_{2} - \frac{g}{L} {\sin θ}_{1} = 0$ $(d theta_2)/(dt) =-b/m theta_2 - g/L sin theta_1 =0$

$\frac{d θ_{1}}{d t} = θ_{2}$ $(d theta_1)/(dt) = theta_2$

The second order differential equation has been converted into a set of two first order differential equations, which can be solved in MATLAB easily.

The conditions of the simulation are as follows:

L = 1 $m$ $m$ , m = 1 kg, b = 0.05 $N s / m$ $Ns//m$ , g = 9.81 $m / s^{2}$ $m//s^2$ .

The intial conditions are as follows:

Angular displacement = 0 $r a d$ $rad$ , angular velocity = 0.5 $r a d / s$ $rad//s$ . (Note: angular velocity of 3 $r a d / s$ $rad//s$ meant that the pendulum did multiple full revolutions around the center.

Code (with explanation)

Function

The purpose of this function is to solve the set of ODEs and return the solution. The input paramets are the time span and initial conditions while the angular displacement and angular velocity of the pendulum are the output parameters.

filename: ode_pendulum.m

% purpose of this function is to solve the set of ODEs and return the solution

function [thetadot] = ode_pendulum (t, theta)

    % initialising input values
    b = 0.05;
    m = 1;
    l = 1;
    g = 9.81;
    
    thetadot = zeros (2,1); 
    
    thetadot(1) = theta(2);  % equation 1
    
    thetadot(2) = (-1*(b*theta(2)/m)-(g*sind(theta(1))/l)); % equation 2

end

Main Program

The main program will call the function, and receive the solution. The received solution will be plotted and the animation will be generated as well.

clear all
close all
clc

% initialising input values
b = 0.05; % damping coefficient (Ns/m)
m = 1;    % mass of pendulum's ball (kg)
l = 0.2;  % length of string of pendulum (m)
g = 9.81; % gravitational acceleration (m/s^2)

% the timespan/time period of oscillation
t = linspace(0,20,100);

% initial conditions. initial angular displacement is 0 rad and initial angular
% velocity is 0.5 rad/s
theta0 = [0;0.5];

% ode function time period in theta1, ode45 solves for displacement and velocity
% in theta(:,1) and theta (:,2) respectively
[t ,theta] = ode45('ode_pendulum', t, theta0);

ct = 1; % counter varaible for animation frames

% for loop to generate animation of the pendulum
for i=1:length(theta)

    x=l*sin(theta(i,1));  
    y=-1*l*cos(theta(i,1));  % calculating coordinates of center of pendulum ball

    plot([-0.1 0 0.1],zeros(1,3),'linewidth',3) % plotting horizontal line
    hold on
    plot(linspace(x,0,5) ,linspace(y,0,5),'color','bla') % plotting string of pendulum
    hold on
    plot(x,y,'o','markersize',20,'color','r') % plotting the pendulum
    hold on
    axis([-0.25 0.25 -0.3 0.1]) % making the pendulum centered in the animation
    axis off % removing uneseccary clutter
    M(ct)=getframe(gcf); % storing each frame in structure M
    ct = ct + 1;
    clf 
end  

close all
% converting frames into a video
videofile=VideoWriter('test.avi','Uncompressed AVI');
open (videofile);
writeVideo (videofile,M);
close(videofile);

% figure showing plot of time vs angular displacement

figure (2)
plot(t,theta(:,1),'color','r','linewidth',1) 
xlabel('Time (s)')
ylabel('Angular Displacement (rad)')

% figure showing plot of time vs angular velocity

figure (3)
plot(t,theta(:,2),'color','b','linewidth',1)
xlabel('Time (s)')
ylabel('Angular Velocity (rad/s)')

Output

Angular Displacement vs Time

Graph of angular displacement vs time is generated, the displacement eventually tends to zero due to damping but this is not visible due to limited simulation time.

Angular Velocity vs Time

Graph of angular velocity vs time is generated, the velocity eventually tends to zero due to damping but this is not visible in this graph due to limited simulation time.

Animation of Pendulum

https://youtu.be/HOCVqsgXsg0

Errors Faced

The function ode_pendulum couldn't take more than two input parameters.

Command Window Screenshot