Week 3 - Solving second order ODEs

Program for solving second order ODE:

ODE function program:

5. Figure(1) is going to tell how damping is done for pendulum. the angular velocity and displacement will come to rest at end of function.

6. Plotting the values for angular velocity and displacement by using plot command.

7. The label command is used to label x and y axis.

8. Legend command is used to label the angular velocity and displacement.

9. To animate the motion of simple pendulum, a forloop is used to create specified set of points to keep bob at different locations.

10. Movie(M) is used to create a film of stored data of the motion of pendulum.

11. VideoWriter command specifies name of the file and type of data is stored in it.

12. Open command is used to open file where the video is going to store.

13. Writevideo is used to create video animation by M file.

14. Close command is used to close the video file.

Output:

1. A plot is drawn between angular velocity and displacement of simple pendulum.

2. Blue curve is displacement where it starts from zero and angular velocity is starting at 3(rad/sec).

Conclusion:

1. We have animated simple pendulum using ode45 solver where the second order is reduced to simple differential equation where it is easy to solve given equation.

2. The animation is made using different points of the bob and motion of simple pendulum is created.

clear all
close all
clc

%solving the second order ODE and show the animation

b=0.05;  %damping coefficient
l=1;  %length of pendulum in metre
m=1;  %mass of bob in kg
g=9.81;  %gravity in m/s^2

%initial condition 
theta_0=[0;3];

%time points
t_span=linspace(0,20,500);

%solve ODE
[t,results]=ode45('ode_function',t_span,theta_0);
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Velocity')
legend('Displacement','Angular velocity')
grid on 
hold off
displacement=results(:,1);
velocity=results(:,2);

%fixed coordinates
x0=1;
y0=1;
ct=1;
for i=1:length(displacement)
    %coordinates
    x1=1-l*sin(displacement(i));
    y1=1-l*cos(displacement(i));
    figure(2)
    %plot fixed link
    plot([ 2.5 -1],[1 1],'linewidth',2)
    hold on 
    plot([x0 x1],[y0 y1],'linewidth',3,'color','b')
    axis([-1 2 -1 2.5])
    plot(x1,y1,'o','Markersize',20,'color','r')
    pause(0.0005)

    M(ct)=getframe(gcf);  %to change the frame of movie
    ct=ct+1;
    hold off

end

movie(M)  %for creating animation video
videofile=VideoWriter('simple pendulum.avi','uncompressed AVI');
open(videofile)  %open the videofile while running the code
writeVideo(videofile,M)  %for writing the videofile
close(videofile)  %for closing the videofile after code run

Program function:
 
function[dtheta_dt] = ode_function(t,theta,b,g,l,m)
b=0.05;  %damping coefficient
l=1;  %length of pendulum in metre
m=1;  %mass of bob in kg
g=9.81;  %gravity in m/s^2
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

Youtube link :

https://youtu.be/UgZ8_jYO2CQ