Week 3 - Solving second order ODEs

AIM: - SIMULATION OF TRANSIENT BEHAVIOUR OF A SIMPLE PENDULAM AND SAVING THE ANIMATION VIDEO

OBJECTIVES:

To write a program in MATLAB to simulate the motion of a simple pendulum, and make a graph of displacement and angular velocity with respect to time by using physical parameters.

l = 1 m length of pendulum

m = 1 kg mass of the bob

b = 0.05 damping coefficient

g = 9.81 m/s² gravity

Simulation of the motion between 0-20 sec. for initial displacement is 0 and initial angular velocity is 3 rad/sec.

DESCRIPTION:

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth.

In engineering systems, ODE (ordinary differential equation) is used to describe the transient behaviour of a system. A very basic example is simple pendulum. The second order ordinary differential equation for motion of simple pendulum is given by.

`(d^2 θ)/(dt^2 )+b/m dθ/dt+g/l sinθ=0`

Since ode can solve only first order derivatives so we have to change above second order derivatives into two first order derivatives.

Let, Ɵ= Ɵ_1, above equation will be

`(d^2 θ_1)/(dt^2 )+b/m (dθ_1)/dt+g/l sinθ_1=0`

Let, `(dθ_1)/dt= θ_2` -------eqⁿ-1

So, `(dθ_2)/dt+b/m θ_2+g/l sinθ_1=0`

`(dθ_2)/dt=-b/m θ_2-g/l sinθ_1` ----------- eqⁿ-2

_{`dθ/dt=d/dt [[θ_1],[θ_2 ]]=[[θ_2],[-b/m θ_2-g/l sinθ_1 ]]`}

STEPS:

1. Start with basic clear all, close all, and clc.
2. Write all input physical data i.e. length of pendulum, mass, damping coefficient, and gravity.
3. Write initial condition of displacement and angular velocity.
4. Make separate calling function to solve ODE i.e. @ode_func.
5. Now use function in main body to solve ode.
6. Plot the graph between time and results.
7. For animation find the coordinates of pendulum.
8. Use for loop command for animation purpose.
9. Get all frame of animation.
10. Make a movie of animation.

CODE USED IN MATLAB

% using second order ODE to calculate pendulam movement with time
clear all
close all
clc
 
% input given
l = 1;              % length of pendulam in m
m = 1;              % mass of the ball in kg
b = 0.05;           % damping coefficient
g = 9.81;           % gravity in m/s^2
 
% initial condition
theta_0 = [0;3];
t_span = linspace(0,20,500);        % time span in s
 
% solve ode
 
[t,result] =ode45(@(t,theta)ode_func(t,theta,l,m,b,g),t_span,theta_0);
figure(1)
sol = result(:,1);
 
%plotting curve
 
plot(t,result(:,1),'linewidth',2)
hold on
plot(t,result(:,2),'linewidth',2)
xlabel('timme in s')
ylabel('displacemment and angular velocity')
legend('displacement','angular velocity')
title('transient behaviour of pendulam')
grid on
% for animation of pendulam
% initial coordinates of string
hold off
ct=1
for j = 1:length(t_span)
    theta3 = result(j,1)
x0 = 0;
y0 = 0;
 
x1 = l*sin(theta3);
y1 = -l*cos(theta3);
 
x2 = -0.5;
y2 = 0;
x3 = 0.5;
y3 = 0;
 
figure(2)
plot([x2 x3],[y2 y3], 'color','b','linewidth',7)
hold on
plot([x0 x1], [y0 y1],'color','k','linewidth',3)
plot(x1,y1,'o','markersize',20,'markerFacecolor','r')
axis([-1 1.5 -1.5 1])
title(['simple pendulam similuation        \theta = ' num2str(sol(ct))])
hold off
pause(0.03)
M(ct) = getframe(gcf)
ct = ct+1
 
end
% creating movie and saving animation
movie(M)
videofile = VideoWriter('pendulam.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

FUNCTION CODE

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

ERROR

ERROR- forgot hold off command and getting wrong animation.

Solved, by using hold off command

OUTPUT

PLOT OF ANGULAR DISPALCEMENT AND ANGULAR VELOCITY WITH RESPECT TO TIME

SYNTAX USED FOR ODE45

[t,result] = ode45(@(t,theta)ode_func(t,theta,l,m,b,g),t_span,theta_0), where t_span = [t0 tf], integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions theta_0. Each row in the solution array result corresponds to a value returned in column vector t.

Where l,m,b,g predefined value to ode_func

FOR ANIMATION FOR LOOP IS USED

for index = values

statements

end

M(ct) = getframe(gcf), for getting all frame for animation

Movie(M), creating movie for M

for animation of pendulum see video below