Week 3 - Solving second order ODEs

OBJECTIVE: To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping and also write a program in Matlab that will simulate the pendulum motion.

CODE 1:

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

CODE 2:

b=0.05
g=9.81
m=1
l=1
theta_0=[0;3]
t_span=linspace(0,20,150)
[t,result]= ode45(@(t,theta)ode_func(theta,b,g,l,m),t_span,theta_0);
plot(t,result(:,1))
hold on
plot(t,result(:,2))

ang_disp = result(:,1)
ct=1
hold off
for i=1:length(t_span)
    x0=0;
    y0=0;
    D=ang_disp(i)
    x1=l*sind(D);
    y1=-l*cosd(D);
    
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',3,'color','g');
    hold on
    plot(x1,y1,'-O','markersize',25,'color','r');
    axis([-0.2 0.2 -1 1])
    
    M(ct)=getframe(gcf);
    ct=ct+1;
    hold off
    
end

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

EXPLANATION: The first step is to change the second order derivative to first order derivative .

d/dt[ theta1 ; theta2 ] = [ theta2 ; -b/m*theta2 -g/l*sin(theta1) ]

where theta will take values of theta1 and theta2.

The ode45 command is used in this program is to perform the first order derivative calculation.

The first step is to assign the function. where b is the damping co-efficient, g is the gravity, l is the length, and m is the mass.

The initial value of angular displacement is assigned as 0 and that of angular velocity as 3rad/sec.

The simulation of the motion is taken between 0-20 secs.

The first order derivative values for various time intervals are been calculated using the ode45 command and is being stored in the array [t, result].

In results the first column gives the value of angular displacement and the second column gives the value of angular velocity with respect to time t.

The values obtained are the plotted in figure1, where the blue curves represent the angular displacement with respect time to time and the red curves represent angular velocity with respect to time.

ANIMATION : The animation of the above plot is being done using the for loop and capturing each values of angular displacement.

Here the bob of the pendulum is being created using the marker command and the size of the bob is being mentioned. The range of the axis is also mentioned.

The obtained frames are being stored in an array called M. The movie command is used in Matlab to create the movie of the obtained frames.

ERRORS OBTAINED : Using the hold on command earlier, which result in multiple images of the pendulum in single figure. It was avoided by using the hold on command after plotting the line, thus by making only the marker to be in motion.

Video Link : https://youtu.be/617aWUDF_iw