Week 4.1 - Solving second order ODEs
clc
clear all
close all
% Parameters
L = 1; % Length of the pendulum (m)
m = 1; % Mass of the ball (kg)
b = 0.05; % Damping coefficient
g = 9.81; % Gravity (m/s^2)
% Initial conditions
theta0 = 0; % Initial angular displacement (rad)
omega0 = 3; % Initial angular velocity (rad/s)
initial_conditions = [theta0; omega0];
% Time span
tspan = [0 20]; % Time range for the simulation (s)
% Define the ODE
ode = @(t, y) [y(2); -b/m * y(2) - g/L * sin(y(1))];
% Solve the ODE
[t, y] = ode45(ode, tspan, initial_conditions);
% Extract the results
theta = y(:, 1);
omega = y(:, 2);
% Plot the results
figure;
subplot(2, 1, 1);
plot(t, theta);
title('Angular Displacement vs Time');
xlabel('Time (s)');
ylabel('Angular Displacement (rad)');
subplot(2, 1, 2);
plot(t, omega);
title('Angular Velocity vs Time');
xlabel('Time (s)');
ylabel('Angular Velocity (rad/s)');
% Animation of the pendulum motion and saving as video
video = VideoWriter('pendulum_motion.avi'); % Create a video file object
open(video); % Open the video file for writing
figure;
for i = 1:length(t)
clf;
x = L * sin(theta(i));
y = -L * cos(theta(i));
plot([0 x], [0 y], 'r-', 'LineWidth', 2);
hold on;
plot(x, y, 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b');
axis equal;
axis([-L L -L L]);
title(sprintf('Time: %.2f s', t(i)));
xlabel('X (m)');
ylabel('Y (m)');
drawnow;
% Capture the frame and write it to the video
frame = getframe(gcf);
writeVideo(video, frame);
end
% Close the video file
close(video);