Week 4.1 - Solving second order ODEs

OBJECTIVE:

The following ODE represents the equation of motion of a simple pendulum with damping:

where,

g = gravity in m/s2,
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.

Write a program in Matlab\Octave that will simulate the pendulum motion.

use,

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.

Use this information to animate the motion of the pendulum.

THEORY:

Given:

`(d^2theta)/(dt^2) + b/m(d theta)/(dt) +g/lsin theta = 0 -----(i)`
`(d theta)/(dt) = 3` rad/`s^2`; when `t = 0`;

Let `(d theta)/(dt) = theta_2 -----(ii)`

and `theta = theta_1`:

`implies (d theta)/(dt) = (d theta_1)/(dt) and (d^2 theta)/(dt) = (d theta_2)/(dt)`

Thus, we get:

`(d theta_2)/(dt) + b/mtheta_2 + g/lsin theta_1 = 0 -----(iii)`

Instead of a second-order ode `(i)`, we have two first-order odes (`ii` and `iii`).

We can now model our program for these two DEs.

CODE:

Equation Solver Function:

function [y] = pendulum(t,theta,g,L,b,m)


theta1 = theta(1);
theta2 = theta(2);

% Angular Displacement
y1 = theta2;

% Angular Velocity
y2 = -(b/m)*(theta2) -(g/L)*sin(theta1);

y = [y1; y2];

end

Calculation and Animation Script:

clear all
close all
clc

# Conditions and parameter definition
tperiod = linspace(0,20,1000);    # 1000 step time-array
theta_init = [0 3];               # Initial Condition
g = 9.81;                         # gravity
b = 0.05;                         # damping
L = 1;                            # length
m=1;                              # mass

# Solving using ode45
[t, theta] = ode45(@(t,theta) pendulum(t,theta,g,L,b,m),tperiod,theta_init);


# plotting
for i=1:length(t)

    x = L*sin(theta(i,1));                 # x-postion of bob at time t
    y = -L*cos(theta(i,1));                # y-position of bob at time t

    plot([-1 1],[0 0],'LineWidth',5,'k')   # Pendulum Base
    hold on;
    plot([0 0],[0 -1],'--k')               # Mean Position
    hold on;
    plot(x,y,'Marker','o','MarkerFaceColor','y','MarkerSize',12)   # bob of pendulum
    hold on;
    plot([0 x],[0 y], 'LineWidth',1,'Color','b')                   # String of pendulum
    hold on;
    pause(0.05)
    hold off;
    filename=sprintf('pendulum/%05d.png',i); # store plot in filename
    print(filename);                         # write plot to storage
end

# Creating Animation
system("ffmpeg -i pendulum/%05d.png -y pendulum/pendulum.mpeg")

OUTPUT: