Week 3 - Solving second order ODEs

Aim: Solving second order differential equation

Objective: To write a program that solves following second ODE by using ode45 solver . This ODE represents the equation of motion of a simple pendulum.

Concept:

This is how we solve 2nd order differential equation in to a 1st order differential equation.

where,

g = gravity in (m/s2),

L = length of the pendulum in (m),

m = mass of the ball in (kg),

b=damping coefficient(N-m/s)

Main Program for the motion of a Simple Pendulum:

% program for solving Second order ODE using ode45 solver

%input variables for second order ODE of a simple pendulum
b=0.05;
g=9.81;
l=1;
m=0.1;

%intial condition
theta_0=[0; 3]; %angular displacement and velocity

%timepoints
t_span=linspace(0,20,500); %motion for 20secs

%solve ode
[t,results]=ode45(@(t,theta)odefunc(t,theta,b,g,l,m),t_span,theta_0);

%plotting the displacement and velocity plot
figure(1)
plot(t,results(:,1)); %plotting angular displacements
hold on
plot(t,results(:,2)); %plotting angular velocity
xlabel("time in sec")
ylabel("displacement(rad)/angular velocity(rad/s)")
legend('displacement','angular velocity')
title("Simple pendulum motion")

angle = results(:,1); %storing the angular displacement in array

ct=1; %intialising the counter

for i = 1 : length(angle)
    omega = angle(i);
    figure(2)
    title('motion of pendulum')
    
    %fixed coordinates of the bob
    x0=0;
    y0=0;
    
    %moving coordinates of the bob
    x1=1*cos((3*pi/2)+omega);
    y1=1*sin((3*pi/2)+omega);
    
    %plotting
    plot([x0 x1],[y0 y1],'linewidth',3,'color','g');
    hold on
    plot(x1,y1,'o','markers',15,'markerfacecolor','k');
    hold on
    plot([-3 3],[0 0],'linewidth',2,'color','b');
    hold off
    axis([-3 3 -2 2])
    grid on
    pause(0.1)
    M(ct) = getframe(gcf);
    ct=ct+1;
end

%creating a movie
movie(M)
videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function Program:

%function program for solving second ode
function [dtheta_dt]=odefunc(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]; %column matrix
end

Explanation:

• Give the input parameter to solve the second ODE: g = gravity in (m/s2),L = length of the pendulum in (m),m = mass of the ball in (kg),b=damping coefficient(N-m/s)
• To use a ode45 solver, create a function program which will return a column matrix for the solver with the reduced 1st order differential equation.
• Provide time span and intial condition(displacement and angular velocity) according to the ode45 syntax.
• Finally, it will give the results in a two separate matrix which will be temporarily stored in workspace and then we use plot command to plot the angular displacement and angular velocity.
• Inorder to capture the variation of time and plot the location of pendulum, the FOR loop command has been used.
• Atlast,create a animation using movie command and the animation video will be stored in the primary location of the saved code file.

OUTPUT:

• The Above picture shows the plot between the angular displacement and angular velocity.
• The Blue curve shows the angular displacement of the bob where it starts from 0 and the other curve is angular velocity which starts from rad/sec.
• The Harmonic motion of the pendulum can be observed where it comes to zero at the end of time span which means it comes to rest due its damping effect.
  

ANIMATION VIDEO:

https://youtu.be/0KFfJ47WBN0

ERRORS OCCURED:

ERROR 1: couldn't able to find what the error was and then matlab itself automatically debugged and it displays the function program of 'update legend menu tool bar' and after that the error was not occured.

ERROR 2: animation was unable to create.

CORRECTION : After running the matlab in adminstrator mode the animation was created