clear; clc; close all; % Problem parameters L = 22.5; % Length of pendulum (m) g = 9.81; % Accel. due to gravity (m/s^2) theta0 = pi/2.5; % Initial angle (rad) theta_dot0 = 0; % Initial angular velocity tspan = [0 150]; % Time interval for simulation %--- CASE 1: NO DRAG (Cd = 0) --- Cd = 0; % ODE options for high and low tolerance opts1 = odeset('RelTol',1e-3,'AbsTol',1e-6); opts2 = odeset('RelTol',1e-6,'AbsTol',1e-9); % ODE function pendulumODE = @(t, y) [y(2); - (g/L)*sin(y(1)) - Cd*abs(y(2))*y(2)]; % Solve at lower tolerance [t1, sol1] = ode45(pendulumODE, tspan, [theta0 theta_dot0], opts1); % Solve at higher tolerance [t2, sol2] = ode45(pendulumODE, tspan, [theta0 theta_dot0], opts2); % Interpolate solution at specific times (e.g., t=22.5, t=50) requested_t = 22.5; theta_interp = interp1(t2, sol2(:,1), requested_t, 'linear'); theta_dot_interp = interp1(t2, sol2(:,2), requested_t, 'linear'); fprintf('Case 1 (No Drag): At t=%.1f s, theta = %.5f rad, dtheta/dt = %.5f rad/s ', ... requested_t, theta_interp, theta_dot_interp); % Plot (t) and (t) figure; subplot(2,1,1); plot(t2, sol2(:,1),'b'); xlabel('t (s)'); ylabel('\theta (rad)'); title('Nonlinear Pendulum Angle, Case 1 (No Drag)'); grid on; subplot(2,1,2); plot(t2, sol2(:,2),'r'); xlabel('t (s)'); ylabel('d\theta/dt (rad/s)'); title('Angular Velocity, Case 1 (No Drag)'); grid on; % Phase plot (Poincar/Phase Map) figure; plot(sol2(:,1), sol2(:,2)); xlabel('\theta (rad)'); y