Week 3 - Solving second order ODEs

Aim:

To write a program to simulate the simple pendulum with damping.

Governing equations:

θ′′+b/m θ′+ gLθ=0

this equation is reduced to two first-order ODE for solving in MATLAB.

Body of content:

MATLAB is used to carry out the simulation. The code is as follows:

clear all
close all
clc

b = 0.05;
g = 9.81;
l = 1;
m = 1;

% Initial condition 
theta_0 = [0;3];

%timespan
t_span = linspace(0,20,500);

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

% figure(1)
% plot(t,results)
% ylabel('plots')
% xlabel('time')
% legend('displacement','velocity')

x0 = 0;
y0 = 0;
ct = 1
for i = 1:length(results)
   
    x1 = -l*cos(results(i,1)); 
    y1 = l*sin(results(i,1));
    
    plot([y0,y1],[x0,x1],'linewidth',2)
    hold on
    plot(y1,x1,'o','MarkerFaceColor','red')
    axis([-2 2 -2 0]) 
    hold off
    pause(0.3)
    M(ct)= getframe(gcf)
    ct = ct+1
   

   
end
movie(M)
videofile = VideoWriter('pendulum_animation.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

The following parameters are taken as constants for a simple pendulum

L = 1 meter.

m = 1 kg.

b = 0.05.

g = 9.81 m/s

The initial conditions are set as displacement = 0 m; angular velocity = 3 rad/s

timelimit is set as 20 seconds.

the result is found by inbuilt ODE solver ode45 where the pendulum is defined in a function as follows:

 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

the ODE solver uses the function with the input parameters of time, initial conditions of displacement, and angular velocity from the array, theta_0, and solves to give a resultant matrix with 2 columns and time span.

since the pendulum is fixed at one point, hence we define two points

x0,y0 = origin.

and the moving potion x1,y1

For loop is initiated from 0 to the length of elements in the result array.

the points x1,y1 defines the position of the pendulum at any time in the time interval.

the values of x1 and y1 are defined as

x1 = -L*cos(results(t,1)), which gives the vertical position of the pendulum at any time t.

y1 = L*sin(results(t,1)), which gives the horizontal position of the pendulum at any time t.

The values are plotted in the 3rd and 4th quadrants and the axis limit is set from x-axis: -2 to 2, y-axis: -2 to 0.

Each fame is saved to make a video simulation of the pendulum motion using the commands, getframe(gcf), and movie respectively.

The simulation video file can be found here:

https://drive.google.com/file/d/1CjbYRAESucCBz8WpNkT6KxbC-XVOwddj/view?usp=sharing