Week 4.1 - Solving second order ODEs
https://youtu.be/1z9F7ivEVeU
Code:
clc; clear; close all;
% Given Parameters
g = 9.81; % Acceleration due to gravity (m/s^2)
L = 1; % Length of the pendulum (m)
m = 1; % Mass of the pendulum bob (kg)
b = 0.05; % Damping coefficient
% Time span for simulation
tspan = [0 20];
% Initial conditions: [theta(0), omega(0)]
theta0 = 0; % Initial angular displacement (rad)
omega0 = 3; % Initial angular velocity (rad/sec)
y0 = [theta0; omega0];
% Solve the ODE using ode45
[t, y] = ode45(@(t, y) pendulumODE(t, y, g, L, b, m), tspan, y0);
% Extract angular displacement
theta = y(:, 1);
% Convert to Cartesian coordinates for animation
x = L * sin(theta);
y = -L * cos(theta);
% Create figure for animation
figure;
axis equal;
xlim([-1.2 1.2]);
ylim([-1.2 1.2]);
hold on;
grid on;
% Plot pivot point
plot(0, 0, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 8);
% Create pendulum components
pendulumLine = plot([0 x(1)], [0 y(1)], 'k-', 'LineWidth', 2);
pendulumBob = plot(x(1), y(1), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 10);
% Animate the pendulum motion
for k = 1:length(t)
set(pendulumLine, 'XData', [0 x(k)], 'YData', [0 y(k)]);
set(pendulumBob, 'XData', x(k), 'YData', y(k));
pause(0.05);
end
function dydt = pendulumODE(~, y, g, L, b, m)
% y(1) = theta, y(2) = omega
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -(b/m)*y(2) - (g/L)*sin(y(1));
end