Solving second order ODEs

AIM

To model and create the animation of a simple pendulum-bob system with damping whose motion is described by the 2nd order ODE -

`(d^2theta)/dt^2+b/m*(d theta)/dt+g/l*sintheta=0`

where

theta - angular displacement of the bob from the mean position

d theta/dt - angular velocity of the bob

t - time

m - mass of the bob

L - length of the string

b - damping coefficient

g - acceleration due to gravity

INTRODUCTION

An ordinary pendulum-bob system consists of a bob of mass m and string of length l. The forces that act on the system are the gravitational forces that act on the bob due to its mass and cause its motion which is countered by the tension force developed in the string that keeps the bob moving in a circular path. There is also a damping force without which the pendulum would keep swinging indefinitely in time. This may be developed due to air resistance or friction at the pivot of the string. This is also factored into the differential equation above which will be used to model the system.

METHODOLOGY

The code starts with the process by first solving the 2nd order ODE to obtain the angular displacement positions (theta) of the pendulum for different instances in time. This is done by the Runge-Kutta method using the ode45 function. This is done as follows

The 2nd order ODE is broken down into 2 single order ODEs through the substitution of new variables

`theta_1 = theta`

`theta_2 = (d theta)/dt`

Therefore the 2 new single order ODEs are

`(d theta_1)/dt = theta_2`

`(d theta_2)/dt = -b/m*theta_2 - g/l*sintheta_1` (bu substitution into the main equation)

These equations are stored in a column vector called dtheta/dt. This is done in a seperate function called ode_func.

The ode45 function takes these equations and along with the time-span and the initial values (i.e, the known values when time = 0), it will provide the results matrix containing 2 columns for angular displacement and angular velocity of the bob at particular time steps within the time span.

STEPS

1. The first step is to input all of the known data values such as the
   l = 1m
   m = 1kg
   b=0.05
   g = 9.81 m/s^2
   
   
2. The initial conditions (i.e, the initial values of theta = 0 and angular velocity = 3 rad/s at t = 0) are stored in theta0 as a colum vector.
   
   
3. The position of the pivot of the pendulum (where the string is attached) is described as (0,0)
   
   
4. The time span of the system's oscillation is stored in the t_span row vector and generated using the linspace command between t = 0 and the given end time of t = 20s.
   
   t_span = linspace(0,20,500)
   
   The time interval is split into 500 steps.
   
   
5. Then the above mentioned ode_func function is created which executes the steps mentioned above that break down the 2nd order ODE into 2 single order ODEs and stores them in the dtheta/dt column vector.
   
   
6. The 2nd order ODE is solved using the ode45 function that is shown below
   
   [t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
   
   The LHS of the equation is the output of the ode45 MATLAB function, which are 2 matrices. t is a column vector that will contain the time span and, results is a 2x500 matrix that contains the angular displacement calculated at appropriate time steps calculated by the ode45 function. It contains 2 columns, the first containing the angular displacements and the 2nd containing the angular velocity for each of the time step instances.
   
   The RHS of the equation contains the inputs to the ode45 function. The first one @(t,theta) which describes the fields involved in the differential equation which are time and theta (angular displacement. The second one is the output of the ode_func function which is the dtheta/dt column vector containing the 2 single order ODEs. The third is the time span and the fourth is the initial values.
   
   
7. The angular displacement and angular velocity values at progressing time steps from the initial condition are stored in the results matrix. In order to develop an animation of the pendulum oscillating, the angular displacement is plotted using the plot function. The string is plotted as a line with one end fixed at (x0,y0) and the other end varying with the angular displacement.
   These varying positions of the other end of the string that is attached to the bob are calculated using a for loop.
   
   The for loop runs from i = 1 to the no of elements in the time span (i.e the number of time steps) and within the loop the position of the bob is calculated using the angular displacement data from the results matrix which is accessed using
   
   ang_disp = results(i,1);
   
   and the following trigonometric relations
   
   x1 = l*sin(ang_disp);
   y1 = -l*cos(ang_disp);
   
   Therefore for each angular displacement at a particular time step, the position of the bob is calculated.
   
   
8. The plot function is used to plot the pendulum with a line representing the string and the bob represented by a circular marker.
   
   plot([x0 x1],[y0 y1],x1,y1,'o','MarkerEdgeColor','r','MarkerFaceColor','r','linewidth',1);
   
   The axis of the plot is also adjusted to allow the proper visibility of the pendulum's oscillation by
   
   axis([-1 1 -1 0]);
   
   After this, a pause of 0.03 seconds is included to ensure each position is clearly visible.
   
   This plot function is also placed within the loop so that it will plot the position of the pendulum for each time step.
   
   
9. In order to include this plot as a frame in the animation video file, it is stored using the getframe command in vector M's first position using the ct (counter) variable. Then in order to ensure that the subsequent frames are stored in adjacent positions on the M matrix, the ct variable is increased by 1 (ct = ct + 1). As the for loop runs for each time step counted out by the i variable, the position is plotted by the plot function and stored within the M vector by getframe until the end of the time span leading to M being a 1x500 vector.
   
   
10. After the end of the for loop, the M matrix is utilized to create the animation video. The VideoWriter command is used to create the name (here, ode_45_animation) and type of the videofile and then the frames in the M matrix are written into the video using writeVideo command.
11. The plots showing the variation of both the angular displacement and the angular velocity are also created using the plot function at the end. This shows the effect of damping on the oscillation of the pendulum system.

The code is shown below

clear all
close all
clc

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

%Initial Condition
theta_0 = [0;3];

%Fixing the base of the string carrying the pendulum
x0 = 0;
y0 = 0;

%Time-points
t_span = linspace(0,20,500);

%Solve ODE

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

%Initiate the counter
ct = 1;

%To create the animation
for i = 1:length(t_span)
    ang_disp = results(i,1);
    
    x1 = l*sin(ang_disp);
    y1 = -l*cos(ang_disp);
    
    %Plot
    plot([x0 x1],[y0 y1],x1,y1,'o','MarkerEdgeColor','r','MarkerFaceColor','r','linewidth',1);
    axis([-1 1 -1 0]);
    grid on
    pause(0.03)
    M(ct)= getframe(gcf);
    ct = ct+1;
end

%To create the animation videofile

movie(M)
videofile = VideoWriter('ode_45_animation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

%To plot the variation of angular displacement with time

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Plot')
xlabel('time')

The code for the function ode_func is shown 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

For the animation video of the pendulum click here