Week 3 - Solving second order ODEs

SOLVING SECOND ORDER ODE's

AIM:

Create Simulation of Simple pendulum With MATLAB as solution for second order ODE.

OBJECTIVE:

Write a program that solves second order ODE for the simple pendulum with given conditions like damping coeffiecient, lenght of pendulum, gravity and mass.

GIVEN VALUES:

Simulate the motion between 0-20 sec

for angular displacement = 0

angular velocity=3 rad/sec

at time t = 0

In the above equation,

g = gravity in m/s^2 = 9.81   

L = length of the pendulum in m = 1

m = mass of the ball in kg = 1

b=damping coefficient = 0.05

THEORY:

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. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

Code for ode_func

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

Main Program

clear all
close all
clc

%given values for damping coeff, gravity accleration, length and mass
b=0.05;
g=9.81;
l=1;
m=1;

%initial conditions
%[angle displacement,velocity]
theta_0 = [0;3];   
    
%time spans
t_span = linspace(0,20,1000);

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

figure(1)
plot(t,results(:,1))       %plot all rows from first column
hold on                    %retain all the data from first plot
grid on
plot(t,results(:,2))       %plot all rows from second column
xlabel('Time')
ylabel('Velocity')
legend('Angular Displacement','Velocity')
hold off



disp = results(:,1);      %all rows from first column which are angle values
vel = results(:,2);


x0 = 0;                      %to set a hinge point
y0 = -0.1;
pr = 1;      

for i=1:length(disp)  
figure(2)               %open in new window, new figure
x1 = l*sin(disp(i));    %x coordinate from trigonometry
y1 = -l*cos(disp(i));    %y coordinate from trigonometry
plot([-2 2],[0 0],'linewidth',20, 'color', 'k')          %plotting wall
hold on
plot(x1,y1,'-O','markersize',20,'markerfacecolor','b')      %plotting the bob with size
plot(x0,y0, '-O','markersize',5,'markerfacecolor','k')    %joining the end points of line and fixed end along with size for the knob
plot([x0 x1],[y0 y1],'g','linewidth',2)                     %line connecting the fixed end and bob with line size
axis([-2 2 -2 2])                                      
grid on

K(pr) = getframe(gcf);
pr = pr+1;
hold off

end


movie(K)                        %starting to create a video file
videofile = VideoWriter('simple pendulum.avi','uncompressed AVI');      %creating a video file named pendulum with avi extension
open(videofile)                 %open the video file
writeVideo(videofile,K)         %write the plot images to the fideo file
close(videofile)                %close the video file. 

ERRORS FACED:

This is because we used d (degree opeator) with sin(Sine in trigonometry) and to make it work properly we just need to remove the d. It happens because the sind(x) represents the angle in degrees and sin(x) represends the angle in radians.

Since the ball is in the negative direction it needs to be negative in numbering.

The operator to be used is writeVideo instead whati used is writevideo which popped an error for creating a video.

RESULTS:

The video link for the movie created is 

Drive Link for the Simple Pendulum Movie

WorkSpace: