Week 4.1 - Solving second order ODEs

Aim: Solving a simple pendulum's second-order ordinary differential equation and animating the damping effect.

Objective: 

• To change an ordinary differential equation of second order into an equation of first order.
• Using the built-in MATLAB function for solving ordinary differential equations to create the function.
• To create a simple pendulum animation.

FACTORS ENTAILED:

• A simple pendulum is a mechanical arrangement that demonstrates periodic motion. The simple pendulum comprises a small bob of mass ‘m’ suspended by a thin string secured to a platform at its upper end of length L.
• The simple pendulum is a mechanical system that sways or moves in an oscillatory motion. This motion occurs in a vertical plane and is mainly driven by gravitational force. Interestingly, the bob that is suspended at the end of a thread is very light somewhat we can say it is even massless. The period of a simple pendulum can be made extended by increasing the length string while taking the measurements from the point of suspension to the middle of the bob. However, it should be noted that if the mass of the bob is changed will the period remains unchanged. The period is influenced mainly by the position of the pendulum in relation to Earth as the strength of the gravitational field is not uniform everywhere.
• By applying Newton’s second law of motion we can obtain the equation of motion for a simple pendulum which is given as:
• $\frac{d^2\theta }{dt}+\frac{b}{m}d\frac{\theta }{dt}+\frac{g}{L}\sin \theta =0$

• The equation must be changed from a second-order partial differential equation to a first-order differential equation in order to be solved in MATLAB using the "ode45" function. we get,

$\theta {\text{'}}_1={\theta }_2$

and $\theta {\text{'}}_2=-\left(\frac{b}{m}{\theta }_2\right)-\frac{g}{L}\sin \theta 1$

Code Explanation:

• The input data from the problems represent the code's initial inputs. The time domain for which the entire code will run is then created. The command "linspace" was used, with a range of 0 to 30 seconds and a value of 100 in between.
• Then, because we need to calculate both values from the equation of motion for plotting purposes, we give the initial conditions for the angular displacement and velocity in column vector form.
• The code includes a function that is used inside the main code to solve the second-order differential equation. Another function in the main code computes the array of time and displacement data.
• After calculation plotting is done in which the subplot command is used, Then, after saving the movie, we created a for loop for which we are creating the frame motion for the straightforward pendulum.

Matlab Program: 

%given data
%damping coefficient
b = 0.05; 
%acceleration due to gravity  
g = 9.81;
%length of Pendulum in meters   
l = 1;
%mass      
m = 1;       
%time duration for simulation 
t_s = linspace(0,30,100);
%angular displacement & velocity
y_starting = [0;3];

%using ODE45
a = @(t,y)ode_second_order_pendulum(y,b,m,l,g);
[t,y]= ode45(a,t_s,y_starting);

%Plot of angular displacement & velocity w.r.t time function.
figure(1);
grid on
subplot(1,2,1)
plot(t,y(:,1));
title('Displacement vs Time');
xlabel('Time in seconds');
ylabel('Angular displacement in rad');
subplot(1,2,2)
plot(t,y(:,2));
title('Angular velocity Vs Time')
xlabel('Time in seconds');
ylabel('Angular velocity in rad/sec');
grid off
saveas(figure(1),'subplot','jpg')

%polar co-ordinate to cartesian co-ordinate
y1 = -1*cos(y(:,1));
x1 = 1*sin(y(:,1));

%video for pendulum motion
for i = 1:100
    figure(2)
    plot([0 x1(i)],[0 y1(i)],'linewidth',3,'color','r');
    hold on
    plot(x1(i),y1(i),'marker','o','linewidth',10,'markerfacecolor','g');
    plot([-1.2 1.2],[0 0],'linewidth',3,'color','b');
    axis([-1.2 1.2 -1.2 1.2]);
    title(['Time: ', num2str(t_span(i))]);
    M(i) = getframe(gcf);
    close(2)
end

movie(M)
videofile = VideoWriter('simplependulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile)

The Function is given as follows;

function ydot = ode_second_order_pendulum(y,b,m,l,g)

ydot = [y(2);-(b/m)*y(2)-(g/l)*sin(y(1))];

end

Output:

• The graph illustrates that the damping factor causes the displacement to reduce as we advance in time. The plot makes it clear that displacement starts at 0 as the initial condition. The angular velocity vs. time plot shows that it began at 3 rad/sec and has since dropped significantly.

• A simple pendulum's second-order ordinary differential equation and animating the damping effect.

Conclusion:

• Using the MATLAB function ordinary differential equations have been created.
• Hence the simple pendulum animation is created.