Week 3 - Solving second order ODEs

The oscilation of the pendulum with respect to the angular displacement and the velocity is been calculated using the ODE45 solver where the result is a form of an array that is plotted and also the position of the pendulum with respect to the time is been animatd using the required code as stated below.

The ode function was developed based on the eqn

there by the equation was broke into smaller single degree eqns for simplicity as follow

theta=theta1

dtheta1/dt=theta2

d^2theta/dt=dtheta2/dt

thereby the the above equation gt reduced to the form of dtheta2/dt=-(b/m)theta2-(g/l)sin(theta1)

this was represented as an ode_func as follows

function [dtheta_dt]= ode_func(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

there by as the main function to implement the ode45 solver was ready,now the main program to find the values of the angular displacement and the angular velocity based on the time as per the requirement.

here ang disp starts from 0 and ang vel from 3rad/sec and time interval for 0-20 sec with 500 values. the following is the program to find the values.

clear all
close all
clc
%b-damping coefficient; g-acceleration due to gravity;
%l-length of the pendulum;m-mass of the ball
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition in the form of initial displacement and velocity
theta_0=[0;3];

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

%solve ODE
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0)

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('displacemnt/velocity')
xlabel('time')
legend 'displacement' 'velocity'

The above code gives the values of the displacement and velocity of the pendulum with respect to the time and is been plotted as follow.

The above plot shows the variation of the displacement and the velocity with respect to time which is being damped at a value of 0.05 damping coefficient, gradually in the given interval.

Now to find the position of the pendulum with respect to time the above code is used along with the video animator code to make it as an animation. Here each value is to be found with respect to the each value of the angular displacemnt which is radians.

x0=0,y0=0 are the fixed positions of the pendulum and x1=l*sin(results(i,1)) and y1=-l*cos(results(i,1)) are the end positions of the free end pendulum as a function of the results obtained. y1 is negative as the fixed point is considered as the origin. The following code will give the entire code for the required animation where the for loop is used to initiate till the end of the length of the array t inorder to access all the angular displacement values .

clear all
close all
clc
%b-damping coefficient; g-acceleration due to gravity;
%l-length of the pendulum;m-mass of the ball
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition in the form of initial displacement and velocity
theta_0=[0;3];

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

%solve ODE
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0)

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('displacemnt/velocity')
xlabel('time')
legend 'displacement' 'velocity'

%now to find the position of the pendulum as per the displacent using the
%moive animator
%ct is a counter used for the movie making with an array M
ct=1;

for i=1:length(t)
    %position of the pendulum fixed at
    x0=0;
    y0=0;
    
    %position of the free end of the pendulum that oscillates 
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));

    %to get the plot in a new window use figure to differentiate.
    figure(2)
 
    plot([-3 3],[0 0],'color','r');

    line([x0 x1],[y0 y1],'linewidth',3);
    axis([-2 2 -2 2])
    hold on
    plot(x1,y1,'.','markers',10,'markerfacecolor','g');
    hold off   
    pause(0.03);
     
    M(ct)=getframe(gcf)
     
    ct=ct+1;
end
%following the syntax for the movie generation using the values in the
%array M
movie(M)
videofile=VideoWriter('ode_main.avi','uncompressed avi')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Here a for loop is used to access all the values in an array and each iteration is counted into the new array M that is used to animate all the values to depict the motion of the pendulum.

Errors faced

-syntax error for the movie making

-getting the 2 plots to in 1 figure

These errors were rectified by refering the videos in the content to create a movie and also the figure command to differentiate the plots.

The above screenshot shows that the both the plots were seen in one.

The above shows the error in the video format code which is rectified

The following is a Youtube link for the animated video that is generate by the above code.

https://youtu.be/gjWV5dRStGQ