Week 3 - Solving second order ODEs

Objective:

write a program that solves the following 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.

Input values:

l=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Solution of ODE:

Step-1: convert 2nd order ODE to 1st order

theta=theta1

dtheta/dt=dtheta1/dt=theta2

d^2theta/dt^2=d^2theata1/dt^2=dtheata2/dt

dtheta2/dt=(-b/m)theta2-(g/l)sintheta

Step-2: creating a function for ODE

function[dtheta_dt]=ode_func(t,theta,b,g,l,m)%creating function named ode_func

theta1=theta(1); %Assign value theta(1) to theta 1
theta2=theta(2); %Assign value theta(2) to theta 2

%find derivation of theta1 and theta2
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);

%creating array of dtheta_dt for function which is column matrix
dtheta_dt=[dtheta1_dt;dtheta2_dt];

end

by this function we can write main program:

% solving 2nd order ode for pendulum
clear all
close all
clc

%input data

b=0.05;
g=9.81;% m/s^2
l=1; %meter
m=1; %kg

%inital condition 1st column represent displacement while 2nd column
%represent velocity
theta_0=[0;3]; 

%time span for pendulum
t_span=linspace(0,20,100);

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

%plotting for analysing pendulum transient and annimation behaviour
figure(1) %creating figure window for plot
plot(t,results(:,1))%ploting graph for time and displacement
hold on  %holding figure window for future plotting
plot(t,results(:,2))% ploting graph for time and velocity
ylabel('Displacement/velocity')
xlabel('time')
ct=1;%assign value 1 to counter
for i=1:length(results(:,1)) %assign value to i

    %origin or base point of pandulum which is fix
    x0=0;
    y0=0;
    %varied value for end point of pendulum
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
      figure(2)% creating figure window for ploting an animation for pendulum
    hold on %hold the figure 2

    plot([-1 1],[0 0],'linewidth',8); %ploting fix roof for pendulum fix point
    plot([x0 x1],[y0 y1],'linewidth',4);%plotting ling or thread for which pendulum ball is hang
    plot(x1,y1,'marker','o','markersize',30,'markerfacecolor','b');%ploting the ball of pendulum
    plot(x0,y0,'marker','o','markersize',9,'markerfacecolor',[1 1 1]); % plotting the fix point of pendulum
    line([x0 x0],[y0 -1.5],'linestyle','--','color','r','linewidth',3);%plotting center line
    
    axis('tight'); %to adject axis
    pause(0.0003);%figure intervel time

    text=sprintf('test%d.png',ct); %creat .png file for each image
    saveas(gca,text);%save all .png file

    clf(); %clear the current graphic image
    d(ct)=getframe(gcf);

    ct=ct+1;%increament in fihure count

  

end
movie(d)
videofile=VideoWriter("animation_pendulum.avi","uncompressed AVI");
open(videofile)
writeVideo(videofile,d)
close(videofile)