Week 3 - Solving second order ODEs

Aim- to solve the second-order differential equation using MATLAB to find out the transient behaviour of a simple pendulum.

Objective- to write a program that solves the following second order differential equation.

Equation-

d^2 theta/dt^2 + (b/m)(d theta/dt) + (g/L)sin theta = 0

In the above equation,

g = gravity (m/s^2)

L = length of the pendulum (m)

m = mass of the ball (kg)

b = damping coefficient

Fixed Parameters-

Length (L) = 1m

Mass (m) = 1kg

Damping coefficient (b) = 0.05

Acceleration due to gravity (g) = 9.81 m/s^2

For Time period = 0-20 sec

Matlab code-

Main code:

Clear all

close all

clc

b = .05; % b=damping coefficient

l = 1; % L = length of the pendulum in m,

g = 9.81; % g = gravity in m/s^2

m = 1; % m = mass of the ball in kg,

% initial condition

theta_0 = [0.5];

% time point

t_span = linspace(0,20,100);

% solve 0DE

[t, result] = ode45(@(t, theta) 0DE_func(t,theta,b,g,l,m), t_span, theta_0);

figure (1);

plot(t,[result(:,1) result(:,2)]);

grid on

ylabel(‘Plot’)

xlabel(‘Time’)

legend(‘Displacement (m) ‘,’Velocity (m/s^2)’);

saveas(figure(1), ‘SHM.jpg’);

theta = result(:,1);

i = 1;

For j = 1:length(theta)

Theta = theta(j);

%initial points

x0 = 0;

y0 = 0;

%movement of pendulum

x1 = l*sin(Theta);

y1 = -l*cos(Theta);

figure (2);

%Plotting string

plot([x0 x1], [y0 y1], ‘color’, ‘g’);

hold on

%Plotting pendulum bob

plot(x1,y1, ‘0’, ‘markersize’, 10, ‘color’, ‘r’, ‘MarkerFaceColor’, ‘r’);

hold off

axis([-1.2 1.2 -1.2 1.2]);

pause(0.001);

M(i) = getframe(gcf);

i = i+1;

end

movie(M);

videofile = VideoWriter(‘pendulum.avi’, ‘Uncompressed avi’);

open(videofile)

writeVideo(videofile,M)

close(videofile);

explanation:-

• In the above program, the fixed parameters are length, mass, damping coefficient, and acceleration due to gravity.
• Time and angle are given as the initial condition.
• The the t_span command, 100 values are taken in between 0 to 20 times period.
• Using the ode 45 command ODE function is introduced.
• Using the plot command velocity and displacement are plotted in y-axis against the time on x-axis.
• Using the for loop command the mass at various points and various times are marked in the graph and are plotted.
• Using the movie command the plot is saved as an compressed .avi file.

function code:ODE_func

function [dtheta_dt] = 0DE_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

• In this function program, the second order ODE is simplified to the first order ODE by introducing the term theta1 & theta2.
• Then theta1 and theta2 are differentiated and the dtheta_dt gives a column matrix.