Week 4.1 - Solving second order ODEs

% 2nd-Order Ordinary Differential Equation for a pendulum

clear
close
clc

b = 0.05;       % Damping Coefficient 
g = 9.81;       % Acceleration due to gravity (in m/s^2)
l = 1;          % Length of the pendulum (in m)
m = 1;          % Mass of the pendulum (in kg)

% Initial condition
theta0 = [0; 3];

% Time points
time_span = linspace(0,20,500);

% Solving ODE
[t, results] = ode45(@(t,theta) ODE_function(t,theta,b,g,l,m), time_span, theta0);

% Plotting the graph
figure(1)
plot(t, results(:,1))
hold on
plot(t, results(:,2))
xlabel('Time')
ylabel('Plot')
legend("Angular displacement","Angular velocity")

% Writing a for loop to make the animation
ct = 1;
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 1;
    
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    % Plotting the animation
    figure(2)
    plot([-1 1],[1 1],'linewidth',5)
    hold on
    plot([x0 x1],[y0 y1],'linewidth',3)
    hold on
    plot(x1,y1,'o','markersize',20,'markerfacecolor','y')
    axis([-2 2 -2 2])
    pause(0.03)
    hold off
    M(ct) = getframe(gcf);
    ct = ct+1;
end

% Converting the animation into a movie
movie(M)
videofile = VideoWriter('Oscillating_Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)





% Solving 2nd order ODE

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

REPORT

Aim : To solve a 2nd Order ordinary Differential Equation of a simple pendulum using MATLAB and create an animation.

Theory :

A second-order ODE is of the form -

A pendulum’s governing equation can be written as -

Acceleration due to gravity, g = 9.81 m/s2

Length of the pendulum, l = 1 m,

mass of the ball, m = 1 kg,

damping coefficient, b = 0.05

θ = θ₁

d(θ)/dt = d(θ₁)/dt = θ₂

d²(θ)/dt² = d²(θ₁)/dt² = d(θ₂)/dt

Thus, we can write the main equation as:

d²(θ₁)/dt² + (b/m)θ₂ + (g/l)sin(θ₁) = 0

Therefore:

D(θ₁)/dt = θ₂

D(θ₂)/dt = -(b/m)θ₂ - (g/l)sin(θ₁)

Code:

Function Program code -

Main program code -

Code to generate and save frame-by-frame images to create video animation -

Code to make and save the video file -

Errors faced were zero.

Outputs :

Animated pendulum, file link -

https://drive.google.com/file/d/1GEscLqKQo-o_WLzI-VTdootv1GRgyhf-/view?usp=sharing

The plotted graph -

Conclusion :

Successfully solves the 2^nd order ODE of a pendulum.

Graph was plotted and animations were created.