Recall the pendulum problem The equation of motion of the sy

Recall the pendulum problem. The equation of motion of the system is given by d^2 theta/dt^2 + g/l sin (theta) = 0 where g is the acceleration due to gravity. We can rewrite second order equations as a system of first order equations and apply solution technique\'s for initial value problems. All we do is introduce a variable to track the first derivative and substitute. So. in this case, introduce theta = alpha obviously, then theta = alpha t. Then, instead of the single second order equation, we have the system theta = alpha alpha = -g/l sin(theta) So now, y = f(l, y) is a vector equation and y = [theta, alpha]^T. But we can still use odc45! Write a MATLAB script, that uses ode45 to solve this system of differential equations for theta(0) = 5 degree, theta(0) = 0. Assume g = 9.81 m/s^2 and l = 1 m. Plot the solution, 0. together with the analytic solution for the small angle approximation, sin(theta) theta. Repeat all steps of the previous part, but now with initial conditions theta(0) = 45 degree, theta(0) = 0.

Solution

pendulum.m file

function dy = pendulum(t,y)
dy = zeros(2,1);
omega = sqrt(9.81);
dy(1) = y(2);
dy(2) = -omega*omega*sin(y(1));

pendulumstc.m file

[t,y] = ode45(\'pendulum\',[0 25], [5 0]);%% change for [t,y] = ode45(\'pendulum\',[0 25], [45 0]) for b%%
plot(t,y(:,1),\'-\');
xlabel(\'time\');
ylabel(\'y_{1}(t)\');
title(\'\\theta (t)\');
figure;
plot(t,y(:,2),\'-\');
xlabel(\'time\');
ylabel(\'y_{2}(t)\');
title(\'d \\theta / dt (t)\');
figure;
plot(y(:,1),y(:,2),\'-\');
xlabel(\'\\theta (t)\');
ylabel(\'d \\theta / dt (t)\');
title(\'Phase Plane Portrait for undamped pendulum\');

run pendulumstc