Solving second order ODEs and Pendulum Animation

Aim: To solve a second order ODE and create an animation simulating the motion of the pendulum

Theory:

ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on Newton's second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena, and with the ability to solve these equations any system can be simulated.

MATLAB’s standard solver for ordinary differential equations (ODEs) is the function
of ode45. This function implements a Runge-Kutta method with a variable time step for
efficient computation.

Syntax for ode45 function

[t,y] = ode45(odefun ,t_span ,yθ)

Data Given:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2,

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

where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b= damping coefficient.

Formulae:

Program:

Code for the main program

% Clearing  the screen, command line, work space and closing all plots
clear all
close all
clc

% Inputs

b = 0.05; % Damping co-efficient
l = 1; % Length of the String in m
g = 9.81; % Gravity in m/s^2
m = 1; % Mass in Kg

% Initial Conditions
theta_0 =[0;3];

% Time
t_span = linspace(0,20,300);

% [t,Results] = ode45('Challenge_2nd_ode',t_span,theta_0)

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

figure(1)
plot(t,Results(:,1));
hold on
plot(t,Results(:,2));
title('Variation of Displacement and Velocity with time');
ylabel('Plot');
xlabel('Time');
legend('Displacement','Velocity');

figure(2)
hold on
ct=1;
for i=1:length(t)
   x0=0;
   y0=0;
   x1=l*sin(Results(i,1));
   y1=-l*cos(Results(i,1));
   plot([-1.5,1.5],[0,0],'linewidth',6,'color','k') % Fixed End
   hold on
   plot([x0,x1],[y0,y1],'linewidth',3,'color','r') % String
   hold on
   plot(x1,y1,'o','markers',20,'markerfacecolor','c') % Pendulum Bob
   axis([-1.5 1.5 -1.5 1.5])
   grid on
   
   pause(0.1)
   hold off
   M(ct)=getframe(gcf);
   ct = ct+1;
   
end

movie(M)
videofile = VideoWriter('Pendulum_Motion.avi','uncompressed avi')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Steps:

1. We define the inputs which are passed onto the function later.
2. Then the initial condition for the angle and time is created in the form of variables as mentioned above.
3. Here we are introduced with the ode45 function whose syntax is mentioned in the theory, where the parameters such as odefun,t_span, and yθ (i.e initial condition) are entered.
4. Then the results are plotted for Displacement vs Velocity by using the code mentioned above.
5. This step is a key one because here we try to create an animation for the pendulum, by using a 'for' loop.
6. The code as mentioned above is written and the plotting is done based on the solution of the ode function.
7. The formatting can be done within such as the string color, thickness(i.e linewidth) and the size of the bob can be adjusted.
8. Then the pause function is used to take the snap of the frame at each position of the pendulum.
9. And after the loop ends with the help of the VideoWriter function we create a video of the pendulum movements.
10. Note: It is very important that this file should be run and not the file where the function is stored.

Code for the function ode_func

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/1)*sin(theta1);
    dtheta_dt = [dtheta1_dt;dtheta2_dt];

end

Steps:

1. The function ode_func is defined with parameters t,theta,b,g,l, and m and then assigned to dtheta_dt.
2. Then the following code is written as shown above.
3. Note: Do not run this code because the parameters are not defined here so it'll give an error.

Errors:

This error has occurred because the file has been run from the ode_func.m file and there are no input arguments which are called so the MatLab shouts an error. To avoid this run the file from the main file where the primary code is written.

Result:

This is the plot that describes the variation of displacement and velocity with time.

Animation Video Link:

https://www.youtube.com/watch?v=ssvKsHRV0_Q&feature=youtu.be