Week 4.1 - Solving second order ODEs

SOLVING SECOND ORDER ODEs

Aim:-

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

Introduction:-

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a 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 for a period of time. The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the ball with respect to time.

A simple pendulum

The figure at the right shows an idealized pendulum, with a massless string of length L and a bob of mass m. 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.

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

where, 'b' is the damping coefficient, 'm' is the mass of the bob, 'g' is the acceleratiion due to gravity, 'L' is the length of the pendulum and 'theta' is the angular displacement.

Similarly, there are hundreds of ODEs that describe key phenomenas , which means ODEs have great importants in engineering.

Function Definition Code:-

function [dtheta_dt] = ode_funct(t,theta,l,m,b,g)

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

Code:-

% SOLVING SECOND ORDER ODE OF A SIMPLE PENDULAM

clear all 
close all 
clc

% inputs for solving second order ODE

l = 1;                         % length of the pendulam in m
m = 1;                         % mass of the ball in kg
b = 0.05;                      % damping coefficient
g = 9.81;                      % acceleration due to gravity in m/s2
t_span = linspace(0,20,400);   % time for which how far the integration is need to be done

% initial conditions
theta_0 = [0;3];               % initial angular dispacement in rad &  angular velocity in rad/s

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

%plotting 
hold on

plot(t,results(:,1));
plot(t,results(:,2));
title('Damping Motion of a Pendulum')
legend('Angular displacement','Angular velocity')
ylabel('Plot')
xlabel('Time(seconds)')

%for loop is introduced as a first step of simple pendulam animation
ct = 1;

for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1))
    y1 = -l*cos(results(i,1))

    % plotting the animation
    figure(2)
 plot([-1 1], [0 0], 'linewidth', 5)
 hold on
 plot([x0 x1],[y0 y1],'linewidth',3)
 hold on
 plot(x1,y1,'o','markersize',20,'markerFacecolor','g')
 title('Animation of Damping Motion of a Pendulum')
 axis([-1.5 1.5 -1.5 1])
 pause(0.01)
 hold off
 
 M(ct) = getframe(gcf)
 ct = ct+1
 
 end
 
% creating an avi file 
movie(M)
videofile = VideoWriter('Animation of simple pendulam.avi', 'Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Code Explantion:-

• Before writing the code, it's important to start with clear all, close all and clc commands.
• Next step is to define the inputs and variables. That are, the damping coefficient, acceleration due to gravity, length of the string and mass of the bob.
• Next step is to define the initial conditions, the angular displacement angular velocity andTimespan have been defined.
• Next step is to solve the ODE using an inbuild matlab function 'ode45'. The matlab function solves an ODE using numerical method. 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.
• The graph has been plotted using the result and time.
• A for loop is introduced at the second phase of the code for finding the coordinate point of the pendulum at its motion. A counter is provided along with the for loop for the execution of 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.
• plotting the support for the pendulum and each position of the bob with respect to time.
• The bob is plotted using the markersize command.
• Axis ranges are provided for the clear view of the animation.
• A pause command was used for pausing the command for a certain time(0.01).
• 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.

Results:-

The figure shows the graph obtained based on the the input values provided including the period of oscillation. 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 animation of the simple pendulum can watch through the below link.

https://youtu.be/ao0zKhsdRWk

Errors Faced:-

This program does not shown any kind of errors in the command window and showing the graph of angular displacement and angular velocity vs time but not the animated motion at the same time, but the value of x1 and y1 are empty vector. then i am so confused.

After a long time checking i realized that the error is in the for loop, becouse the program shows the above mentioned graph. Then i finded the mistakes in the 37 th and 38th line. The mistake i was done is provided the range operator instead of comma in the both lines. After the changing the range operator with comma the error was solved.