Week 3 - Solving second order ODEs

ODE is used to describe the transient behavior of a system. A simple example is a pendulum. The way pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

g = gravity in m/s2,

L==r= length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

To Write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with a damping step hence solve it.

Defining the initial values,

Inputs:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Initial conditions: Simulating the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0

When I did this I got the position of the pendulum as a function of time. Thus I wrote the code to plot the graph as well as to simulate the pendulum according to the graph

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

%Second ODE 

clear all
close all
clc

%inputs
b = 0.05;   %Damping Coefficient
g = 9.81;   %Gravity m/s^2
m = 1;      %Mass of the Pendulum in kg
l = 1;      %Length of the string attached to pendulum in meters

%initial conditions

theta_0 = [0;3];  %[angular displacement=0;angular velocity=3]   

%time points
t_span = linspace(0,20,500); %pendulum oscillates between 0-20 seconds

%solve ODE

[t,result] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

%plotting the graph
figure(1)

subplot(2,1,1)
p1 = plot (t,result(:,1),'linewidth',2,'color','b')
xlabel('Time')
ylabel('Theta')
title('Angular Displacement(theta) v/s Time')

hold on 
subplot(2,1,2)
p2 = plot(t,result(:,2),'linewidth',2,'color','r')
xlabel('Time')
ylabel('Velocity')
title('Angular Velocity v/s Time')
legend([p1(1),p2(1)],'angular displacement','angular velocity')

%pendulum simulation using ode
theta = result(:,1)
ct=1;
%%creating a for loop for oscillating the pendulum

for i=1:length(result)
    THETA = result(i,1);
    
    %initial boundary conditon
    x0=0;
    y0=0;

    %Defining x and y coordinates
    x1=  l*sin(THETA);
    y1= -l*cos(THETA);
    
    figure(2)
    
    hold on
    plot([x0 x1],[y0 y1],'linewidth',4);
    plot(x1,y1,'-o','markersize',30,'markeredgecolor','g','markerfacecolor',[0.6 0.6 0.6]);
    plot(x0,y0,'-*','markersize',3,'markeredgecolor','g','markerfacecolor',[0.6 0.6 0.6]);
    axis([-2 2 -2 0]);
    pause(0.005);
    
    
    M(ct)=getframe(gcf)
    ct=ct+1;
    clf()         %clf() code is  used for clearing the current graphic image
end
%%creating the animation and movie(avi)

movie(M);
videofile = VideoWriter('second_ode_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)