Week 3 - Solving second order ODEs using MATLAB

Solving second order Ordinary Differential Equation using the ODE 45 solver in MATLAB

The scope of this project is to solve the second order ODE system of a simple pendulum and represent it as an animation for a specified range of time and other known input parameters.

To start with the following inputs are known:

Length of the pendulum l = 1 m

Mass of the bob m = 1 kg

Damping Coeffecient b= 0.5

Acceleration due to gravity g = 9.81 m/s^2

Time of simulation t = 0 - 20 s

Any differential euqation can be solved only if its values at the boundaries are known. In this case, at time t=0, angular displacement=0,angular velocity=3 rad/sec. This is represented as a row vector theta_0.

theta_0=[0;3];

Then, the time range is split into 200 segments of equal value. ( Each time step has a value of 0.1 s). In order to obtain a better resolution of the results, the size of the time step can be further reduced.

The next major command in this program is to extract the values of displacement at specific time steps suing the ode45 solver. This command is an inbuilt function which takes in the dependant variables, a function to represent the differential equation, time span and the boundary values of theta at time 0 s.

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

Function to represent the differential equation:

The differential equation is modeled as a seperate function so as to provide it as an input to the ode45 command.

% Function to represent the equations of the simple pendulum

function [dtheta_dt]= pendulum_ode(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

Another major part of this program is to interpret the results graphically in the form of an animation. In order to do this, the following steps are used:

1. A counter variable is initialized first (ct=1). A for loop is then declared which iterates for as many number of times as the number of values in the time span.

2. Appropriate expressions which represent the position of the pendulum with respect to time are then declared.

3. Using the 'line', 'markers' and plot commands, the string, the bob and the fixed point of the simple pendulum are represented.

4. The 'hold on' command makes sure that each values of each time step is plotted on the same figure. The 'hold off' command ensures to remove the previous plot after a time step in order to show the next plot.

5. Then these frames are individually captured using the getframe command.

M(ct)=getframe(gcf)

6. These individual frames are then written to a videofile . The following snippet is used to accomplish this:

% Storing the captured frames in a videofile
movie(M)
videofile=VideoWriter('simple_pendulum_animation')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Finally, an AVI file with the animation of the simple pendulum is stored in the current working directory.

Actual code:

% A program to solve the ODE of a simple pendulum and represent it as an
% animation

clear all
close all
clc

% Known inputs

l=1; % Length in m
m=1; % Mass in kg
b=0.5; % Damping coeffecient
g=9.81; % Acceleration due to gravity in m/s2

%Initial condition
theta_0=[0;3];

% Time values
time_span=linspace(0,20,200);

%Using inbuilt ODE solver
[t,results]=ode45(@(t,theta) pendulum_ode(t, theta, b,g,l,m),time_span,theta_0);

%plotting the damping 
figure(1)
plot(t,results(:,1),'color','b')
hold on
plot(t,results(:,2),'color','r')
xlabel('Time')
ylabel('Displacement and Velocity')
legend('Displacement','Velocity')


% Plotting the pendulum figure using markers and getting the frames using
% loop

ct=1
for i=1:length(t)
    
x0=0;
y0=0;
x1=l*sin(results(i,1));
y1=-l*cos(results(i,1));
figure(2)
 plot([-1 1],[0 0], 'linewidth',2)
line([x0 x1],[y0 y1],'linewidth',4)
hold on
plot(x0,y0,'o','markers',7,'markerfacecolor','y')
plot(x1,y1,'o','markers',20,'markerfacecolor','r')
axis([-2 2 -2 2])
hold off
M(ct)=getframe(gcf)
ct=ct+1

end

% Storing the captured frames in a videofile
movie(M)
videofile=VideoWriter('simple_pendulum_animation')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Results:

A plot of displaced and velocity against time was also obtained in order to visualize the damping.

The animation was uploaded to YouTube:

https://youtu.be/-kPdNOKXFYM

Errors faced while programming:

One of the major erros face was due to the improper naming of the function used with the ode45 command. When calling a seperate function within the current program, the file containing the function must have the same name of the function. Otherwise errors would occur.