Week 3 - Solving second order ODEs

Aim : To solve an ordinary differential eqution (ODE) and use the calculation information to animate the motion of the pendulum.

Introduction : Pendulum 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.

rough diagram

Objective :

Write a program that solves the following ODE . This ODE represents the equation of motion of a simple pendulum with damping

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg.

b = damping coefficient

Parameters

Length (L) = 1 meter

Mass (m) = 0.1 kg

Damping coiefficient (b) = 0.02

Acceleration due to gravity (g) 9.81 m/s2

Matlab Code :

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 .

Here is Main program

%Solving second order ODEs 

clear all 
close all
clc

b= 0.02;
g= 9.81;
l= 1;
m= 0.1;

% initial condition 
theta_0=[0;5];

% time points
t_span = linspace(0,20,50);

% solve ODE
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
hold on 
plot(t,results (:,1))
plot (t,results(:,2))
title ('Damping Harmonic motion of a Pendulum')
legend ('Angular displacement','Angular velocity')
ylabel('plot')
xlabel('time')

% 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 
      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 Animatin.avi','Uncompressed avi')
  open(videofile)
  writeVideo(videofile,M)
  close(videofile)

Explanition :

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 equation using a numerical method . The numerical method in question here is the 4th order Runge Kutta method . The function ode45 taken 3 inputs .The three inputs are the function to be integrated, the time span and the initial condition. [ t, results ] means that the results 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 .

The second part of the code is to establish a for loop , A for loop is a repetition control steucture that allow you to efficiently write a loop that 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 [ x0 x1 y0 y1]and their corresponding relation are entered to from 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 projest 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 for loop run. End command was used for the plot and for the 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 used to close the video file .

Error :

Not end the For loop causing error.

Conclution :

The x_start, y_start , x_end , y_end values

x_start = 0

y_start = 0

x_end = 0.9274

y_end = -0.3740

Animation Video link :