Week 3 - Solving second order ODEs

AIM:

Write a program in MATLAB that will simulate the Pendulum motion.

GOVERNING EQUATION:

Equation of motion of a simple pendulum with damping:

`(∂^2 θ)/(∂t^2 )+b/m*(∂θ)/(∂t)+g/L sin⁡θ=0`

Solving this equation helps to get displacement (`θ`) and rate at which the displacement changes (velocity).

Braking 2^nd order ode into 1^st order ode:

Assuming, `θ=θ_1`

`( ∂θ)/(∂t)=(∂θ_1)/(∂t)=θ_2` → 1

`(∂^2 θ)/(∂t^2 )=(∂^2 θ_1)/(∂t^2 )`

`∂/(∂t) ((∂θ_1)/(∂t))=(∂θ_2)/(∂t)`→ 2

Substituting equation 1 & 2 in 2^nd order ode:

ie, `(∂θ_2)/(∂t)+b/m θ_2+g/L sin⁡θ_1 =0`

1^st ODE: `(∂θ_1)/(∂t)=θ_2`

2^nd ODE: `(∂θ_2)/(∂t)=-b/m θ_2-g/L sin⁡θ_1`

In the above equations,

• g- Gravity `(m/s^2 )`
• L- Length of the pendulum (`m`)
• m- Mass of the ball (`kg`)
• b- Damping coefficient

OBJECTIVE:

• Solve the ODE.
• Simulate the pendulum motion using ODE.

PROGRAM:

This program helps to find the changes in angular velocity and angular displacement with respect to time. It also animates the pendulum oscillation.

ODE45:

Ode45 is a solver in matlab which solves 1^st order ode using Runge Kutta method.

Syntax of ode45: [t,results] = ode45(@(t,y) ode(t,y,args),tspan,y0)

Ode45 is an inbuilt function with inputs t (time) & y (value of function). Since the t & y are not known another function called ode is introduced which contain same input along with other arguments. Here t is the current time and y is the value of function at previous time step. Then time span (tspan) is given which is nothing but limits for the integral. y0 is the initial condition for the function at t=0second.

Explaining the function code:

Defining the function ode45 with another function called ode_func.

%Function for ODE
function [dtheta_dt]=ode_func(t,theta,b,g,l,m)
         theta1=theta(1);
         theta2=theta(2);
         dtheta1_dt=theta2; %Equation for 1st ode 
         dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1); %Equation for 2nd ode
         dtheta_dt=[dtheta1_dt; dtheta2_dt]; %Creating column matrix

end

• Line 2: Starting the function with function name ode_func, input arguments (t, theta, b, g, l, m) are given.
• Line 3-4: Assigning the values of angular displacement and angular velocity to theta1 and theta2. This is done to avoid error because of neglecting the brackets.
• Line 5: 1^st ode equation.
• Line 6: 2^nd ode equation.
• Line 7: Creating a column matrix using 1^st and 2^nd ode equation.
• Line 9: Function is ended.

Explaining the Main code:

clear all
close all
clc

%Program to slove the 2nd order ode.

%inputs
b=0.05; %Damping coefficient
g=9.81; %Gravity (m/s^2)
l=1;    %Length (m) 
m=1;    %Mass(kg) 


%Initial Condition
theta_0=[0;3]; %initial values of displacement and velocity

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

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

%plotting
figure(1) %1st figure 
hold on   
plot(t,results(:,1))%Plotting time vs angular displacement
plot(t,results(:,2))%Plotting time vs angular velocity

%Naming x and y axis
xlabel('Time')
ylabel('Angular displacement & velocity')
grid on
legend('Angular dispalcement','Angular Velocity')%Labeling the curves

%For Loop
ct=1; %counter
for i=1:length(results(:,1));
    %The pivot point
    x0=5;
    y0=5;
    %The points through which the ball oscillates
    x1=x0+l*sin(results(i,1));
    y1=y0-l*cos(results(i,1));
    
    %plotting
    figure(2)
    plot([4 6],[5 5],'linewidth',2,'color','y');%plotting the line where the pendulum pivoted
    hold on
    plot([x0 x1],[y0 y1],'linewidth',2,'color','m');%Ploting a line from pivot to the ball
    plot(x1,y1,'.','MarkerSize',60,'color','g'); %Plotting the position of ball
    title('Pendulum Motion')
    axis ([4 6 4 6])
    axis equal   
    hold off
    M(ct)=getframe(gcf); %Storing each frame as an array
    ct=ct+1; %Incrementing the counter
    
        
end
%Making video file
movie(M);
videofile=VideoWriter('Pendulum.avi','uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

• Line 8-11: Inputs to the program, which are damping coefficient, gravity, length of pendulum and mass of the ball.
• Line 15: Giving the initial condition for the two 1^st order ode. The initial displacement is 0 and velocity is 3`((rad)/s)`.
• Line 18: Giving the limits to integrate the ode with respect to time. Linspace command is used to give time steps. Smaller the time step gives more accurate approximation of integral.
• Line 21: Calling the inbuilt ode45 function. ‘t’ and ‘theta’ are the inputs for ode45. Ode_func is another function to define ode45. Integral limits and initial conditions are also given.
• Line 24: New figure window is created.
• Line 25: hold on command for holding the current plot and adding the next plot to the existing plot.
• Line 26: Plotting a graph of time vs 1^st column in the results which is nothing but angular displacement.
• Line 27: Plotting a graph of time vs 2^st column in the results which is nothing but angular velocity.
• Line 30-31: Labeling X & Y axis using xlabel & ylabel.
• Line 32: Creating grid for the plot.
• Line 33: Using command legend labeling each curve.
• Line 36: Figure counter.
• Line 37: Starting the For Loop. Here For loop counter (i) is initialized as 1 to length of the array called results.
• Line 39-40: Giving the pivot point of pendulum.
• Line 42-43: Finding the co-ordinates of ball at each position of angular displacement.
• Line 46: Creating new figure.
• Line 47: Plotting a base line where the pivot is positioned.
• Line 48: hold on command for adding the subsequent plots to existing plot(plot created in line 47).
• Line 49: Plotting a line from the pivot to each value of x1 and y1, this forms the string of pendulum. Line width is set to 2 and color of line as magenta.
• Line 50: Plotting marker at each value of x1 & y1, with size set as 60 and color of marker as green.
• Line 51: Title for the plot is given.
• Line 52: Axis limit is set.
• Line 53: Here aspect ratio of the plot is adjusted automatically according to the range of data.
• Line 54: Hold off reset the axis properties.
• Line 55: Getting each frame and storing it as an array inside the variable M.
• Line 56: Incrementing the counter.
• Line 59: Ending the For Loop.
• Line 61: Making a movie from the frames stored in variable M.
• Line 62: Creating a video file using command VideoWriter and giving the file name and profile for video.
• Line 63: Opening the video file.
• Line 64: Writing the frames captured (M) to the video file.
• Line 65: Closing the video file.

Errors faced:

1. Error: Error in pendulum animation

Solution: Marker should be placed before marker size to avoid this error.

2. Error: Bouncing feeling during the oscillation of pendulum. (Length variation)

Solution: Error was due to improper axis limits. When I corrected the axis limits and introduced axis equal command, the error was solved. Axis equal command helps the program to adjust the aspect ratio automatically.

RESULT:

Fig: Amplitude VS Time plot

Here at initial stage amplitude of both velocity and displacement are in maximum eventually due to damping the amplitude is decreasing with respect to time.

Simulation of Pendulum:

https://drive.google.com/file/d/1IPc8u9db4X6d9ZBY3yG7757YnoNawFB3/view?usp=sharing