Week 3 - Solving second order ODEs

OBJECTIVE:

To write a MATLAB program that solves the ODE which represents the equation of motion of a simple pendulum with damping.

Using, 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. Also, animate the motion of the pendulum using the results obtained.

SOLUTION:

The aim of this program is to solve an ordinary 2^nd order differential equation, which represents the motion of a simple pendulum. The length of pendulum arm is 1m, mass = 1kg and damping coefficient = 0.05. The motion is simulated for 20 sec starting from rest, where initial conditions were taken as: angular displacement = 0 rad and the angular velocity as 3 rad/s.

Program:

%Program for solving 2nd order ODE

clear all
close all
clc

%intial boundary conditions
theta_0 = [0;3];  %angular displacement = 0, angular velocity = 5

b = 0.05;   %damping factor
g = 9.81;   %gravity
l = 1;      %length of pendulum
m = 1;    %mass

t_span = linspace(0,20,500);   %definition of time array

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

figure(1)
plot(t,results(:,1)); %plotting time vs angular displacement
hold on 
plot(t, results(:,2)); %plotting time vs angular velocity
xlabel('time');
ylabel('Plot');
legend('Angular dispalcement','Angular velocity');

%animation

k = 1; %counter 
for i= 1:length(results(:,1)) %loop for generating the animation
    
    x0 = 0;   %defining the fixed end of pendulum arm
    y0 = 0;
    
    x1 = l*sin(results(i,1));    %defining the moving end of pendulum arm
    y1 = -l*cos(results(i,1));
    
    figure(2)
    
    plot([-1 1], [0 0], 'linewidth', 5, 'color', 'r'); %plotting the fixed beam
    hold on;
    plot(x0,y0,'o','markersize',10,'markerfacecolor','b'); %plotting the pivot point
    hold off;
    axis([-2 2 -2 2]);
    line([x0 x1], [y0 y1], 'linewidth', 3, 'color', 'b'); %plotting the pendulum arm
    hold on;
    plot(x1,y1,'o','markersize',30,'markerfacecolor','g'); %plotting pendulum bob
    hold off;
    
    M(k)= getframe(); %storing all frames in an array
    k=k+1;
    
end

%creating the animation movie
movie(M);
videofile = VideoWriter('Pendulum2.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

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

Explanation:

In this program, we solve a 2^nd order ordinary differential equation by reducing it into two first order DE’s. The initial boundary conditions was defined at t=0, where angular displacement =0 and angular velocity = 3. This was stored in an array named, theta_0 and another array t_span was created with the time values. In order to solve the problem, the ode45 solver was invoked as follows:

The ode45 function call is explained below:

1. [t, results] : the outputs of the function is stored here
2. ode45(@(t,theta) : the solver was invoked and the variables are t & theta
3. ode_func(t,theta,b,g,l,m) : the function name is ode_func and the parameters inside the brackets are the arguments

4. t_span,theta_0 : the initially defined arrays, with the boundary conditions and time values are passed into the function

Using the results obtained from the function, a plot is generated between the angular displacement and angular velocity against time and this plot is stored as figure (1).

For creating the animation, for loop was defined and the iterations was carried out for each displacement value. A new frame is plotted for each displacement value, using the x & y coordinates of both ends. Each of these frames generated was store in array M, and the it was converted into an animation using the ‘movie’ command.

Result:

1. Angular dispalcement vs time, angular velocity vs time

2. Animation

The animation can be found in the following youtube link:

https://youtu.be/K02JZG5MzLE