Week 3 - Solving second order ODEs

% Solving second order ODEs
l = 1;
m = 1;
b = 0.05;
g = 9.81;

% initial value of theta at t = 0
theta_0 = [0;3];

t_span = linspace(0,20,500);

% to solve ode
[t, results] = ode45(@(t,theta) ode_func(theta,b,g,l,m), t_span, theta_0);
plot(t, results(:,1));
hold on
plot(t, results(:,2));
legend('displacement', 'velocity');
pause(2);
hold off

%to plot the pendulum animation
x0 = 0;
y0 = 0;
c = 1;
for i = 1:length(t_span)
   plot([-1.1 1.1],[0 0] ,'LineWidth', 20);
   axis([-1.5 1.5 -1.5 0.5]);
   hold on
   x1 = l*sin(results(i, 1));
   y1 = -l*cos(results(i, 1));
   plot([x0 x1],[y0 y1], 'LineWidth', 5);
   plot(x1,y1,'o',MarkerSize = 20,MarkerFaceColor = 'b');
   plot(x0,y0,'o',MarkerSize = 10,MarkerFaceColor = 'g');
   pause(0.05);
   M(c) = getframe(gcf);
   c = c+1;
   hold off
end
movie(M);
video = VideoWriter('pendulum animatioin', 'Uncompressed AVI');
open(video);
writeVideo(video,M);
close(video);

function [dtheta_dt]  = ode_func(theta, b, g, l, m)
theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end

output:

pendulum animation link

Explanation-

initial theta value  = 0;

inital velocity = 3m/s

t_span is a time matrix that has 500 values between 0 and 20

in results the first column has the values of angle theta and the second column has the values of velocity.

ode45 is an in-built MATLAB function. ode_func is a user made function to get the equation in a suitable form for solving in MATLAB.(the single second order equation is coverted into two first order equations)

results(:,1) represents the entire first column. Both the first column and the second are plotted with respect to time.

'hold on' is used so that both the displacement and velocity plots come in a single plot and then hold off is used to release the plot for next plot.

the tip of the pendulum is plaved at (x0,y0)=(0,0). c is the counter variable for frames. i is a loop control variable.

an arbitrary base is plotted using the some arbitrary coordinates which represent a horizontal line.

hold on is used to keep the base intact in the next plots. x1 and y1 are the other end points of the string. a bob is made at x1,y1 by using the marker 'o'.

M is a matrix that contains all the frames. c is incremented each time at the end to change the frame. 'hold off' will release the present plot so that the next plot does not get overwritten.

a movie is made with the stored frames in M. The video name is pendulum animation.

theta(1) is an angle but theta(2) is velocity so it is a differential of theta(1). So, dtheta1_dt is stores theta2. so different variables are used to convert one second order ODE into two first order ODEs. 

the values of dtheta1_dt , dtheta2_dt is passed to the main function and ode45 solves the converted equatioins.