Week 3 - Solving second order ODEs

Objective:

To solve a second-order ordinary differential equation of a simple pendulum and show the effect of damping by an animation.

Description:

The ODE represents the equation of motion of a simple pendulum with damping:

d^2θ/dt^2 + b/m * dθ/dt + g/L * sinθ = 0

In this equation;

g = gravity in m/s2.

L = length of the pendulum in m.

m = mass of the ball in kg.

b=damping coefficient.

L=1 metre.

m=1 kg.

b=0.05.

g=9.81 m/s2.

The resultant first-order equation will be:

θ'1 = θ2 and θ'2

= -bmθ2 - gLsinθ1

Code explanation:

• We will provide the inputs to begin the code. Then we make a time domain for which the whole code is going to run. Using the ‘linspace’ command from 0 to 20 seconds with 50 values in between them.
• After that, we provide the initial condition for the angular displacement and velocity, which is column vector form, because we need to calculate both values from the equation of motion for plotting purposes.
• The code comprises of function which solves the second-order differential equation which is taken inside the main code. The main code has another function that calculates the array of data for time and displacement.
• After calculation plotting is done in which subplot command is used, in the first plot first column is plotted and in the second subplot, the second column is plotted with respect to time.
• Then we are changing the polar co-ordinate to cartesian co-ordinate because we have used the ‘sin’ command instead of the ‘sind’ command in function.
• Then we created a for loop for which we are creating the frame motion for the simple pendulum and saved it in the movie.

Function code:

function ydot = ode_second_order_pendulum(y,b,m,l,g)
ydot = [y(2);-(b/m)*y(2)-(g/l)*sin(y(1))];
end

Main code:

close all
clear all
clc

%Given data
b = 0.05;   %damping coefficient
g = 9.81;   %Acceleration due to gravity
l = 1;      %Length of Pendulum in meters
m = 1;      %Mass 

%Time duration for simulation 
t_span = linspace(0,20,50);
%Providing Initial Condition as Angular displacement & Velocity
y_initial = [0;5];

%Calculation using ODE45
a = @(t,y)ode_second_order_pendulum(y,b,m,l,g);
[t,y]= ode45(a,t_span,y_initial);

%Plot of Angular displacement & Velocity with respect to Time function.
figure(1);
grid on
subplot(1,2,1)
plot(t,y(:,1));
title('Displacement vs Time');
xlabel('Time in seconds');
ylabel('Angular displacement in rad');
subplot(1,2,2)
plot(t,y(:,2));
title('Angular velocity Vs Time')
xlabel('Time in seconds');
ylabel('Angular velocity in rad/sec');
grid off
saveas(figure(1),'subplot','jpg')

%Polar Co-ordinate To Cartesian Co-ordinate
y1 = -1*cos(y(:,1));
x1 = 1*sin(y(:,1));

%Video for Pendulum motion
for i = 1:50
    figure(2)
    plot([0 x1(i)],[0 y1(i)],'linewidth',3,'color','r');
    hold on
    plot(x1(i),y1(i),'marker','o','linewidth',10,'markerfacecolor','g');
    plot([-1.2 1.2],[0 0],'linewidth',3,'color','b');
    axis([-1.2 1.2 -1.2 1.2]);
    title(['Time: ', num2str(t_span(i))]);
    M(i) = getframe(gcf);
    hold off
end

movie(M)
videofile = VideoWriter('singlependulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile)

Plotting graph:

Animation:

https://drive.google.com/file/d/1KQCZ1qwiEnmQc_XDQLYFSuaJk2Z6Kp5g/view?usp=sharing

Result:

We have created a code to simulate and animate the second-order ODE.