Week 4.1 - Solving second order ODEs

% Pendulum Simulation

% Parameters
g = 9.81;       % Gravity (m/s^2)
L = 1;          % Length of the pendulum (m)
m = 1;          % Mass of the ball (kg)
b = 0.05;       % Damping coefficient

% Initial conditions
theta1 = 0;             % Initial angular displacement (rad)
theta2 = 3;             % Initial angular velocity (rad/s)

% Time span for simulation
t_span = [0, 20];        % Simulation time interval (s)
dt = 0.1;               % Time step for animation (s)

% Define the ODE function
odefunc = @(t, y) [y(2); -(g/L)*sin(y(1)) - (b/m)*y(2)];

% Solve the ODE
[t, results] = ode45(odefunc, t_span, [theta1; theta2]);

% Extract the angular displacement and calculate the x-y coordinates
theta = results(:, 1);
x = L * sin(theta);
results = -L * cos(theta);

% Create an animation
figure;
for i = 1:length(t)
    % Plot the pendulum
    plot([0, x(i)], [0, results(i)], 'k', 'LineWidth', 2);
    hold on;
    plot(x(i), results(i), 'ro', 'MarkerSize', 20, 'MarkerFaceColor', 'r');
    hold off;
    
    % Set plot properties
    xlim([-L-0.5, L+0.5]);
    ylim([-L-0.5, 0.5]);
    axis equal;
    title('Pendulum Motion');
    xlabel('x (m)');
    ylabel('y (m)');
    
    % Pause for animation
    pause(dt);
    
    % Create movie frames
    frames(i) = getframe(gcf);
end

% Create the animation movie
animation = VideoWriter('pendulum_animation.mp4', 'MPEG-4');
animation.FrameRate = 10;
open(animation);
writeVideo(animation, frames);
close(animation);