Solving second order ODEs

Introduction:

A simple pendulum consists of a mass m hanging from a string of length L and fixed at a pivot point P. When displaced to an initial angle and released, the pendulum will swing back and forth with periodic motion.

Aim:

To write a program in MATLAB that will simulate the motion of a simple pendulum.

Objective:

To write a Matlab program to solve an ODE that represents the equation of motion of a simple pendulum with damping.

Formulae & Explaination:

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Procedure:

1. Note down the inputs (b,g,l,m).
2. Create a function ode_func to find dtheta/dt.
3. Input the variables in the main program such as the damping coefficient, gravity constant, mass and Length.
4. Set the range for the time span.
5. Use solver (ode45) and call the function ode_func to solve the ODE.
6. Plot figure() Time vs Plot.
7. Assign the counter to 1 and initialize the first coordinates.
8. Calculate the final coordinates
9. Plot the figure for the motion of the pendulum.
10. Save the animation file by entering the code to save it as movie.

Program code:

%Program to simulate the motion of a simple pendulum using MATLAB
clear all
close all
clc

%inputs b,g,l,m:
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%Initial condition
theta_0 = [0;3];

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

%solve ODE by using inbuilt solver and a function created
[t, results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span,theta_0 );

%Plotting the graph for time vs plot
 figure()
 hold on
 plot(t,results(:,1))
 plot(t,results(:,2))
 ylabel('Plot') % Labelling the y axis
 xlabel('Time')      % Labelling the x axis
 title("Characteristics of simple pendulum") 
 grid on

 %Plotting the motion for the simple pendulum
 %First coordinates
 x0 = 0;
 y0 = 0;
 ct = 2;
for i = 1:length(results(:,1))
    %First coordinates
     x0 = 0;
     y0 = 0;
    %final coordinates
    x1 = 1*sin(results(i,1)); %final coordinates x axis
    y1 = -1*cos(results(i,1));  %final coordinates y axis
    
    %PLot for the motion of pendulum
    figure(ct)
    hold on
    plot([-1 1],[0 0],'linewidth',2,'color','black')
    
    line([x0 x1], [y0 y1], 'linewidth',2,'color','blue')
    hold on
    
    %Bob of the pendulum
    plot(x1,y1,'-0','markersize',15,'markerfacecolor','r')
    grid on
    
    %Setting up the axis
    axis([-1.5 1.5 -1.5 1.5])
   
    M(ct)=getframe(gcf);
    ct = ct+1;
end

%Saving the animation file
movie(M)
videofile = VideoWriter('pendulum.avi','uncompressed avi')
open(videofile,M)
writeVideo(videofile,M)
close(videofile)

Function Code:

%Function code to solve an ODE.
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

Explanation of the Code:

The given details are the inputs for the program. Now we set the range for the time span, we then use the ode45 solver and the function ode_func to solve the equation.

We now Plot the figure Time vs Plot, after this we use the looping statements to obtain the simulation of the pendulum. x0 ,y0 are the initial coordinates and x1,y1 are the final coordinates.

The bob of the pendulum is shown in the plot by using Marker command. The videoWriter command creates a video file of the animation.

Errors Occurred:

An error had occurred while calling the function, I solved it by entering another variable ‘theta_0’ for the function ode_func.

I entered ‘sine’ command in an equation but, matlab uses the command ‘sin’ to find sine of the angle.

Result:

The solution of an ODE that represents the equation of motion of a simple pendulum with damping is obtained and is shown as follows:

Simulation output:

Conclusion:

From the plots and simulation we an conclude that, the relationship between the periods of oscillation of a simple pendulum with its length from the Centre, which increases with the length. The time also varies as a square of the length of the pendulum.