Week 4.1 - Solving second order ODEs

AIM:-

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

where,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

t= 0 to 20 sec

theta=[angular_displacement angular_velocity]

at t=0, theta=theta0=[0 3]

PROCESS:-

• Here I have written 2 programs one is for the ode function of pendulum and other is for the animation of pendulum
• I have given step by step explanation of both the program in side the program in comment format.
• Also attached animation of pendulum, snapsorts of the plots and program for your reference.

1st Program :- (ODE function of pendulum)

% Writing a function for ODE which is used to describe the transient
% behavior of a simple pendulum.
function [dtheta_dt]= ode_func_pendulum(t,theta)

% fixed values
b=0.05;
m=1;
L=1;
g=9.81;

% function statements
theta1=theta(1);
theta2=theta(2);
dtheta_dt=zeros(2,1);
dtheta_dt(1)=theta2;
dtheta_dt(2)=(-b/m)*theta2-(g/L)*sin(theta1);
end

2nd Program:- (Animation of Pendulum)

clear all
close all
clc
% g = gravity in m/s2,
% L = length of the pendulum in m,
% m = mass of the ball in kg,
% b=damping coefficient.
L=1;
% time span from 0 to 20
t=linspace(0,20,150);
% theta array contains angular displacement and angular velocity
% at t=0,theta=theta0,angular displacement=0 & angular velocity=3 rad/sec
theta0=[0 3];
% calling the function containing ODE function for pendulum
[t theta]=ode45('ode_func_pendulum',t,theta0);
% Subplot for angular displacement
subplot(2,1,1)
plot(t,theta(:,1),'-b','linewidth',3);
xlabel('time t in sec');
ylabel('Angular displacement');
% subplot for angular velocity
subplot(2,1,2) 
hold on
plot(t,theta(:,2),'-k','linewidth',3);
xlabel('time t in sec');
ylabel('Angular velocity');
grid on

X0=0;
Y0=0;

% creating variable ct to store each position of pendulum
ct=1;
% using for loop to calcute aall the positions of pendulum from 0 to 20 sec
for i=1:length(t)
    X1=L*sin(theta(i,1));
    Y1=-L*cos(theta(i,1));
    
    figure(2)
    % Plotting the line or bar containing pendulum string
    plot([-1 1],[0 0],'linewidth',6,'color','g')
    axis([-3 3 -3 3]);
    % Plotting the string containing pendulum
    line([X0 X1],[Y0 Y1],'linewidth',2,'color','b');
    hold on
    % plotting the joint of string with bar
    plot(X0,Y0,'*','markersize',5,'markerfacecolor','y');
    % Creating the pendulum
    plot(X1,Y1,'O','markersize',18,'markerfacecolor','r');
    % Giving title to the plot
    title('Damping of simple pendulum');
    % legend is used here to distinguish the parts of the pendulum
    legend('support rod','Length of pendulum','Fixed point','Ball of Pendulum');
    grid on;
    hold off;
    % storing each ct or each frame in 'M' array
    M(ct)=getframe(gcf);
    % Counter weight increment
    ct=ct+1;
end
% below options are used to create a movie
movie(M);
%Create a video file
videofile=VideoWriter('Simple Pendulum.avi','Uncompressed AVI');
% To open the videofile
open(videofile)
% To write the movie frame on videofile
writeVideo(videofile,M)
close(videofile)

Error Faced by me:-

• I did't have any idea about how to create the bar, string and pendulum and to make them in their position
• Also have some confusion about the theory behind it

Solution:-

• Fot these above problems I have gone through previous video lecture on ODE and used youtube to know how to solve above problems using MATLAB.