Project

Solving Second Order ODEs using Matlab

MMurtadha Zainy7 Aug 2020

Aim:-

To solve an ordinary differential equation (ODE) and use the calculated information to animate the motion of the pendulum via MatLab.

Introduction:-

Pendulums have many applications and were utilized often before the digital age. They are used in clocks and metronomes due to the regularity of their period, in wrecking balls and playground swings, due to their simple way of building up and keeping energy.

A simple pendulum consists of a point mass suspended on a string or wire that has negligible mass. If the pendulum weight or bob is pulled to a relatively small angle from the vertical and let go, it will swing back and forth at a regular period and frequency. Figure 1 shows a simple pendulum model

The figure at the right shows an idealized pendulum, with a "massless" string or rod of length L and a bob of mass m. The open circle shows the rest position of the bob. When the bob is moved from its rest position and let go, it swings back and forth. The time it takes the pendulum to swing from its farthest right position to its farthest left position and back to its next farthest right position is the period of the pendulum.

The primary forces acting on the bob are the gravitational force that makes it move in the first place and the force exerted by the string to keep it moving along a circular path. In addition, there may be a damping force from friction at the pivot or air resistance or both. We will construct a model to describe how the angle theta of the pendulum varies as a function of time t.

Let s(t) be the distance along the arc from the lowest point to the position of the bob at time t, with displacement to the right considered positive. Let theta(t) be the corresponding angle with respect to the vertical. The figure shows tangential and radial components of the gravitational force on the pendulum bob. The radial component is exactly balanced by the force exerted by the string, so the only relevant force producing the motion is the tangential component of the gravitational force. For the moment, we ignore the damping force, if any.

An ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The ODE equation and the calculations are shown below in Figure 2.

Code Explanation:-

The function code is shown in figure 3

% A function program which will be called in the main program to solve the
% ODE


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

A function is a group of statements that together perform a task. In MATLAB, functions are defined in separate files. The name of the file and of the function should be the same. Functions operate on variables within their own workspace, which is also called the local workspace, separate from the workspace you access at the MATLAB command prompt which is called the base workspace.

Function handles are simple to function codes that can be created without using the regular function definition using the keywords ‘function’ and ‘end’. They do not require any separate ‘.m’ files or not needed to be declared separately below/above the main code. Function handles are directly declared in the main code itself and also accessible from any point in that particular code, i.e they are localized to a particular conditional. The declaration of a function handle is very similar to that of a function declaration. It can be distinguished by the means of ‘@’ symbol before your value/expression allocation. The syntax for creating a function handle is. A function was defined in order to establish a relation between all the given input values, the function code is based on the equation derived in figure 2.

Figure 4 shows the main code of the project.

% Objective: Solving second order ODEs of a penedulum


clear all
close all
clc


b = 0.05 % damping coefficient
g = 9.81 % gravity (m/s^2)
l = 1   % length of the spring (m)
m = 0.1  % mass of the bob (kg)

% initial condition
% angular displacement is 0 and angular velocity is 5

theta_0 = [0;5]

% time span
% Pendulum oscillates between 0-20 seconds

t_span = linspace(0,20,600)

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

%Plotting
figure(1)
hold on
plot(t,results(:,1))
plot(t,results(:,2))
title('Damped Harmonic motion of a Pendulum')
legend('Angular displacement', 'Angular velocity')
ylabel('Plot')
xlabel('Time(seconds)')



% For loop 
% Before starting the for loop, a counter value was defined

ct =1;
 for i = 1:length(results(:,1))
 x_start = 0     
 y_start = 0
 x_end = l*sin(results(i,1))
 y_end = -l*cos(results(i,1))
 
 % plotting the animation
 figure(2)
 plot([-1 1], [0 0], 'linewidth', 5)
 hold on
 plot([x_start x_end],[y_start y_end],'linewidth',3)
 hold on
 
 plot(x_end,y_end,'o','markersize',20,'markerFacecolor','y')
 title('Animation of transient behaviour of a Pendulum')
 axis([-1.5 1.5 -1.5 1])
 pause(0.03)
 hold off 
 
 M(ct) = getframe(gcf)
 ct = ct+1
 
 end
 
 
% Animation of the pendulum

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

Before writing the code, it's important to start with clear all, close all and clc commands since

Clear all - removes data from the workspace, freeing up system memory

Close all- deletes the current figure or the specified figure

clc - clears all input and output from the Command Window display, giving you a "clean screen

Next step is to define the inputs and variables. The damping coefficient, gravity, length of the spring and mass of the bob have all been defined.

Next step is to define the initial conditions, the angular displacement and angular velocity have been defined. Timespan has also been defined.

Next step is to solve the ODE, Ode45 was used, it is a function in Matlab that solves ordinary or partial differential equations using a numerical method. The numerical method in question here is the 4th order Runge Kutta method. The function ode45 takes 3 inputs. The 3 inputs are the function to be integrated, the time span and the initial conditions. [t, results] means that the results that have been calculated with respect to time.

After calculating the results, a graph has to be plotted, the graph contained two plots. One for angular displacement with respect to time and one for angular velocity with respect for time. The legend for the graph and the axis' labels have been defined. Figure 5 shows the graph obtained

The second part of the code is to establish a for loop, A for loop is a repetition control structure that allows you to efficiently write a loop that needs to execute a specific number of times. a counter loop was added before the For loop. the start and endpoints for x and y-axis have been defined. The loop is used where all the values of timespan are assigned to a variable 'i' one at a time. The coordinates of

[ xo x1 y0 y1] and their corresponding relations are entered to form a 'string' of the pendulum. Markersize command was used to form the bob of the pendulum.

A pause command was used for pausing the command for a certain time, for this project a pause of 0.03 s was used. The M command was used to capture the current frame every time the for loop runs. The final step of the loop is to add a +1 to the ct variable; this will increment the frame every time the for loop run. End command was used for the plot and for the for loop to stop both from running.

The final step of the code is to create a movie of the animation of the robotic arm, so the movie command was used. The VideoWriter command was used to assign the video format, for this project .avi was used as the format. The open command was used to open the video. writevideo command was used to write the video with plot frames. Finally, the close command was used to close the video file.

Conclusion:-

The x_start,y_start, x_end y_end values are shown below in figure 6

x_start =

     0


y_start =

     0


x_end =

   -0.0113


y_end =

   -0.9999

Animation Video:-

Here is the link for the video of the animation that I got for this project.

https://www.youtube.com/watch?v=JnsGTYADL6I

References:-

Solving Second Order ODEs using Matlab

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.