Week 3 - Solving second order ODEs

Program to simluate the simple harmonic motion of a pendulum:

A simple pendulum's motion is based on Newton's second law of motion.

The equation of motion is given as follows

where,

theta = angular displacement(in radians)

b = damping coeffecient

m = mass of the bob(in kg)

g = acceleration due to gravity(in m/s^2)

L = length of the string(in m)

ode45 syntax :

[t,results] = ode45(@(function_variable) function_defenition ,t_span,initial_values)

ode45 is a function in MATLAB that solves ODE(ordinary differnetial equations)or PDE(partial differential equations) using numerical methods.

To solve and ode we need to integrate it or we convert it into a linear equation and solve it by numerical methods.

ode45 uses a function called anonymous function which can be declared within the code itself.It does not need regular function defining keywords such as ' function ' and ' end ' .

ode45 solves the given ode by using Runge kutta method

Animation using For loop:

We give the condition to for loop such that every possible displacement of the pendulum is captured in the frame. So i = 1 : lenght(results(:,1)) means that 1st column of the result array which basically consist all possible values of displacement so the loop runs till all possible values of pendulum displacement is plotted.

Then we define the coordinates of the pendulum as

x = 0 .y = 0

x1 = l*sin(results(i;1)) [if i =1 then first value of displacement is taken as theta]

y1 = -l*cos(results(i;1)) [y1 is negative beacuse the pendulum hangs in the downward direction]

Base of the pendulum is plotted using convineint values

String of the pendulum is plotted using the x,x1,y,y1 values

Bob of the pendulum is plotted using x1 and y1 values

getframe() command is used to capture the frames in an array and movie() command is used to conert the captured frames into a movie

clear all
close all
clc

%Inputs

%Damping coefficient
b = 0.05;
%Gravity(in m/s^2)
g = 9.81;
%Length of the pendulum(in m)
l = 1;
%Mass of the bob(in kg)
m = 1;

%Initial values
theta_0 = [0;3];

%Pendulum oscillates between 0-20 secs
t_span = linspace(0,20,500);

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

%Plotting
subplot(5,4,[13,14,17,18]);
%Time vs Angular displacement
plot(t,results(:,1))
hold on
grid on
%Time vs Angular velocity
plot(t,results(:,2))
hold off
xlabel('plot')
ylabel('time (in sec)')
legend('Angular displacement','Angular velocity') 
subplot(5,4,[15,16,19,20])
%Angular displacement vs Angular velocity
plot(results(:,1),results(:,2))

%Declaring counter variable
ct = 1;
%Looping
for i=1:length(results(:,1))
    
    %Coordinates of the pendulum
    x = 0;
    y = 0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    
    %Plotting
    subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12])
    %Base of pendulum
    plot([-1 1],[0 0],'linewidth',4,'color','b')
    hold on
    %Massless string of the pendulum
    line([x x1],[y y1],'linewidth',2,'color','k')
    hold on
    %Bob of the pendulum
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r')
    grid on
    title('Simulation of oscillating simple pendulum')
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.03)
    hold off
    %Capturing the frames in an array
    M(ct) = getframe(gcf);
    ct = ct + 1;
end

%Turning the captured frames into a video
movie(M);        
videofile = VideoWriter('Oscillation of simple pendulum','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

    
    

The given ode is of the order 2. For MATLAB to solve the equation the given ode should be break down into firast order ode.

Let theta = theta1

dtheta/dt = dtheta1/dt = theta2-------------->Equation 1

d^2theta/dt^2 = d^2theta1/dt^2 = dtheta2/dt---------------->Equation 2

Substituting Eq 1 and 2 in the given equation of motion we get

dtheta2/dt = -(b/m)theta2-(g/L)sin(theta1)--------------->Equation 3

%Function to convert the given ODE into first order ODE's
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

Errors faced while programming:

1. In the line y1 = -l*cos(results(i;1)) did not mention the negative sign resulting in pendulum facing the upward direction.

2. Used l*cos(resutls(i;1)) for x1 leading to the pendulum swinging upwards.

3. Basic errors like missing paranthesis.

Animation: https://youtu.be/gAV0UXpY3-g