Week 4.1 - Solving second order ODEs

Aim: To write a program that solves the following ODE.

Mathematical calculation

`(d theta)/(dt) = theta_1`

`(d theta_1)/(dt) = theta_2`

`(d^2theta_2)/(dt^2) =(d theta_1)/(dt)`

substituting these values in the main equation of Pendulum we wil get

`(d theta_1)/(dt)+b/mtheta_1+g/Lsintheta = 0`

Matlab code

Function code

function [xdot] = Damping_pendulum(t,x);%using x as theta
% inputs
L = 1;
m=1;
g=9.81;
b=0.05;
xdot = zeros(2,1);%assigning nyll matrix for xdot
%rate of change of displacement= velocity
xdot(1) = x(2);%first element is angular velocity which is at y(2)position
%
xdot(2) = (-b/m)*x(2)-(g/L)*sin(x(1));%this will obtaian after conversion of 2nd order ode into 1st order
end

Main code

%solving a second order ode
clear all 
close all
clc

L = 1;%string length
t =linspace(0,20,150);%time array
x0 = [0 3];%position array
%caling function code to solve using ode45
[t,x] = ode45('Damping_pendulum',t,x0);
%Plotting for angular displacement and angular velocity 
subplot(2,1,1)
plot(t,x(:,1),'b','linewidth',3)
xlabel('time t in sec')
ylabel('angular displacement')
subplot(2,1,2)
hold on;
plot(t,x(:,2),'k','linewidth',3)
xlabel('time t in sec')
ylabel('angular velocity')
grid on; 
X0=0;%initial co-ordinate of x for string
Y0=0;%initial co-ordinate of y  for string

ct=1;%counter of every frame
for i=1:length(t) %for lop for each iteration of time
    X1=L*sin(x(i,1));%co-ordinate of x for  string end  at every  iteration
    Y1=-L*cos(x(i,1));%co-ordinate of y for  string end  at every  iteration
%     Plotting the output of each teration
    figure(2)
    plot([-1 1],[0 0],'linewidth',6,'color','black');%plotting support rod
    axis([-2 2 -2 2]);%defiig axis
    line([X0 X1],[Y0 Y1],'linewidth',2,'color','blue')%string
    hold on;
    plot(X0,Y0,'O','markersize',3,'markerfacecolor','yellow');%fixed point
    plot(X1,Y1,'O','markersize',20,'markerfacecolor','red');%bob of pendulum
    title('damping of simple pendulum');%title
    legend('support rod','length of Pendulum','Fixed point','Bob of pendulum')%legends
    grid on;
    hold off;
    pause(0.01)
    m(ct)=getframe(gcf);%geting each frame if iteration for creating animation
    ct = ct+1;%updating frame countend%end of loop
end
%creating animation
movie(m)
videofile = VideoWriter('Damping_of_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,m)
close(videofile)

Output plot

The above plot shows the time vs angular displacement and angular velocity respectively.

Animation link.

https://youtu.be/eciSTckjfUI