Week 4.1 - Solving second order ODEs

Code with Explanation :

1. Simple_pendulum.m

clear all

close all

clc

% Define initial conditions and parameters

l = 1; % Length of the pendulum

m = 0.1; % Mass of the pendulum bob

b = 0.05; % Damping coefficient

g = 9.81; % Acceleration due to gravity

theta_0 = [0; 5]; % Initial angle and angular velocity

t_span = linspace(0, 10, 500); % Time span for ODE integration

% Solve the ODE using ode45 and ode_func function

[t, results] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);

% Plot the angular displacement over time

plot(t, results(:, 1))

hold on

plot(t, results(:, 2))

pause(0.03)

ylabel('Angular Displacement')

xlabel('Time')

% Iteration using a for loop to create an animation of the simple pendulum

ct = 1; % Counter for frames

theta_1 = results(:, 1);

for i = 1:length(results(:, 1))

% Initial coordinates of the pendulum

x0 = 0;

y0 = 0;

theta_2 = theta_1(i);

% Final coordinates of the pendulum bob

x1 = l * sin(theta_2);

y1 = -l * cos(theta_2);

% Plotting for animation

figure(2)

% Base of the pendulum

plot([-1, 1], [0, 0], 'linewidth', 4, 'color', 'black')

hold on

% String of the pendulum

line([x0, x1], [y0, y1], 'linewidth', 6, 'color', 'b')

hold on

% Bob of the pendulum

plot(x1, y1, 'marker', 'o', 'markersize', 28, 'markerfacecolor', 'r')

axis([-2, 2, -2, 2])

pause(0.03)

hold off

% Capture the frame for the animation

M(ct) = getframe(gcf);

ct = ct + 1;

end

% Create a movie from the captured frames

movie(M)

% Save the animation as an uncompressed AVI file

videofile = VideoWriter('simple_pendulum.avi', 'uncompressed AVI');

open(videofile);

writeVideo(videofile, M);

close(videofile);

2. ode_func.m

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

Output video :

https://drive.google.com/file/d/1fZ-yokaHjGJIZ9I71Rd_tOSVAVDh8Oy5/view?usp=drive_link