Week 3 - Solving second order ODEs

MATLAB program to simulate the pendulum motion:

The motion of the pendulum swing can be modeled by the equation:

`frac{d^2theta}{dt^2}+frac{b}{m}frac{d theta}{dt}+frac{g}{l}sin theta=0`

Here `theta : text{angular position},b:text{damping coefficient}, m:text{mass of the ball}`

`L:text{Length of the pendulum}, g:text{gravity}`

Now, the pendulum motion is simulated using MATLAB for a time of 0-20 seconds, initial angular velocity as 3 radian per second at time t=0.

The MATLAB code is given below:

%solving 2nd order ode to simulate the pendulum motion

%input parameters assumed
b=0.05; % damping coefficient
g=9.81; % acceleration due to gravity (m/s^2)
l=1;% length of the pendulum (m)
m=0.1; % mass of the pendulum (kg)

%initial conditions 
theta_0= [0;3]; % first element represents the initial position and second element represents initial velocity.
% time points
t_span=linspace(0,20,500); % Motion simulated for a total of 20 seconds with 500 time points

%solving the ODE
% here ode45 is an inbuilt matlab function to solve ordinary differential
% equation.
% first argument represents the independent variables, here time and theta
% second argument represents the differential equation/function that needs to be solved
% third argument represents the total time span for which the function has
% to be solved
%Fourth argument represents the initial conditions
% the output to this function is a column vector of time points and a matrix of position and velocity points. 
[t,results]=ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0)

pos=results(:,1) % storing the position vectors from the 'results' matrix
ct=1; % creating a counter to record the animation


% To show the animation, we can loop over a vector of time values, i.e for
% each time step, a line segment is drawn joining the fixed point and the free end which is represented by the soltuion of the differnetial equation.
for i=1:length(pos)  %To loop over the time position points
    ang=pos(i)       % extracting the position or angle for each loop count
    x0=0; % the x - coordinate of the fixed end of the pendulum
    y0=1;  % the y - coordinate of the fixed end of the pendulum
    
    x1= l*sin(ang);  %the x -coordinate of the free end of the pendulum
    y1= l*(1-cos(ang)); %the y -coordinate of the free end of the pendulum
    
    
    
    plot([x0 x1],[y0 y1],'linewidth',3) % plotting a line between the end points
    hold on
    
    plot(x1,y1,'o','MarkerSize',15,'MarkerFaceColor',[1 0.2 0.8]) % plotting the ball at the free end of the pendulum
    hold off
    axis([-1 1 -0.5 1.5]) % setting up the axis limits
    pause(0.02) % pausing the time frame
    M(ct)=getframe(gcf); % Capturing the present figure and storing it in a variable M
    ct=ct+1; % giving the counter next value so that the next frame can be stored accordingly
end
%recording the movie
movie(M)
videofile=VideoWriter('pendulum_swing.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

The code for the function representing the differential equation is given below:

function [dtheta_dt]=ode_function(t,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

The animation of the pendulum swing is available through this link:

https://drive.google.com/file/d/1dEPa2mUFk6uCThA4GwgDOCOBu-t0gpil/view?usp=sharing