Week 3 - Solving second order ODEs

Introduction:-

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position.

Aim:-

To solve the second order ODE function using Matlab that will simulate the pendulum motion.

Objective:-

To write a program that solves the second order ODE.This ODE represents the equation of motion of a simple pendulum with damping

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

The parameters given for the above equation

g = 9.81m/s2

L = 1 metre

m = 1 kg

b = 0.05

Angular displacement = 0 ,Angular velocity = 3 rad/sec at time t=0

Function 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

Explanation:-

• In this function program, the second order ODE is simplified to the first order ODE by introducing the term theta1 & theta 2.
• Then theta1 & theta2 is differentiated and the dtheta_dt gives a column matrix.

Main program:-

clear all
close all
clc

% 2nd order ODE 

b = 0.05;  % b is damping coeff
g = 9.81;  % g is gravity
l = 1;      % l is length pendulum
m = 1;   % m is mass

%intial condition

theta_0 = [0,3];

% time points

t_span = linspace(0,20,200);

% solve ODE

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

%plotting

figure(1)

plot(t, results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('velocity(red),dispalcement(blue)')

% animation

ct = 1;
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1)+pi);
    y1 = l*cos(results(i,1)+pi);
    
    figure(2)
    
    plot([x0 x1],[y0 y1],'-o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r','linewidth',4)
    axis([-2 2 -2 2])
    hold on
    plot([-1 1],[0 0],'linewidth',5,'color','[0 0.3 0]')
    M(ct) = getframe(gcf)
    hold off
    ct=ct+1;
end

% Movie& Videofile

movie(M)
v_f = VideoWriter('Pendulum_ODE.avi','Uncompressed AVI')
open(v_f)
writeVideo(v_f,M)
close(v_f)
    

Explanation:-

• In the above program, the fixed parameters are length,mass,damping coefficient,and acceleration due to gravity.
• Tinme and angle are given as the intial condition.
• In the tspan command,200 values are taken in between 0 to 20 times period.
• Using the ode45 command ODE function is introduced.
• Using the plot command velocity and displacement is plotted in the y-axis against time in the x-axis.
• Using the for loop command the mass at various times are marked in the graph and are plotted.
• Using the movie command the plot is saved as an uncompressed avi file.

Plot:-

Animation:-