Week 3 - Solving second order ODEs

SOLVING 2ND ORDER ORDINARY DIFFERENTIAL EQUATIONS (ODE) OF A SIMPLE PENDULUM

AIM - To simulate the motion of a simple pendulum using MATLAB

OBJECTIVE :-

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

• Animate the motion of the pendulum.
• Upload the animation on You tube

THEORY:-

• A pendulum is basically 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.

• A simple pendulum is an ideal pendulum consisting of a point mass suspended by a weightless inextensible perfectly flexible thread and free to vibrate without friction.
• It generally includes a lot of equations, from which one of them is given in the objective above. We are going to solve that one in MATLAB and in order to simulate the pendulum, an animation will be created.

CODES :-

%To simulate the pendulum motion
clear all
close all 
clc

%Inputs (damping coefficient=unitless,gravitational acceleration=m/s^2)
%(length of pendulum=m, mass of ball=kg)
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%initial conditions
theta_0 = [0;3];

%time points
time_span = linspace(0,20,500);

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

%Plotting
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
axis([0 20 -4 4])

ct = 1;
theta3 = results(:,1);
for i=1:length(theta3)
    
    %Inputs required to initiate for loop
    x0 = 0;
    y0 = 0;
    THETA3 = theta3(i);
    
    x1 = l*sind(THETA3);
    y1 = -l*cosd(THETA3);
    
    %Plotting
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',3)
    hold on
    plot(x1,y1,'marker','o','markersize',30,'markerfacecolor','g')
    axis([-0.03 0.03 -2 0])
    pause(0.03)
    hold off
    
    %Getting frames of the animation
    M(ct) = getframe(gcf);
    ct = ct+1;
    
end
%Creating Movie
movie(M)
Videofile = VideoWriter('ode_solver.avi','Uncompressed avi');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)

CODE EXPLANATION:-)

• We started the codes by providing fixed terms as inputs such as damping coefficient (b), gravitational acceleration(g), length of the pendulum(l) and mass of the ball(m).
• The initial conditions(i.e., at time,t = 0) for the simple pendulum are angular displacement = 0 and angular velocity = 3 rad/sec. These conditions are specified in terms of a column matrix 'theta_0'.
• To find solution of the ODE, we need to integrate it for time 0-20 seconds. This is specified in terms of linspace command by the variable 'time_span'.
• Now, in MATLAB, there are lots of solvers from which we are using the solver 'ode45'. This solver needs the differential equations to be solved, which is provided with the help of a function 'ode_func'. Some other codes are written to define this function as below:

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 differential equations to be solved are specified in these codes.
• The integrating limits and initial conditions are also provided to the solver.
• The solutions of the differential equations are plotted in the next step. Angular displacement and angular velocity, both are plotted in a single plot:

• Now, to generate an animation of a pendulum, we need a variation of the angular displacement, obtained by solving the differential equations, with time. A specific value of time corresponds to a particular position of the pendulum, which is basically the angular displacement.
• To fulfill above criteria, 'for' loop is used. Each value of the angular displacement is recorded in variable 'THETA3'.
• (x0,y0) is the fixed end of the pendulum and (x1,y1) is the moving end of the pendulum, which is actually the sine and cosine of the length of the pendulum. Negative sign shows that the pendulum is below x-axis.
• Then, the pendulum is animated on the plot in the next step. After capturing the frames, the movie is generated with the name 'ode_solver.avi'
• https://www.youtube.com/watch?v=LoLt1achsLo and the above plot are the results obtained by executing the program.

ERRORS ENCOUNTERED:-

1)

• Here, the size of the marker is in inverted commas, which are required to be removed.

2)

• This is just a spelling mistake. Actual code is 'getframe'.

3)

• In this error, y1 is taken as l*cos(theta3).
• However, it should be -l*cos(theta3) as it is below x-axis.

CONCLUSION - Thus, we have solved ODE of a simple pendulum as well as created its animation in MATLAB.