Solving Second Order ODE (Simple Pendulum)

Problem Statement :-

Write a Program to solve following equation.

`frac {d^2 theta}{dt^2} + frac {b}{m} frac {d theta}{dt} + frac{g}{L} sin theta = 0`

This ODE represents the equation of motion of simple pendulum with damping.

MATLAB program should contain simulation the motion of pendulum between 0 to 20 seconds.

Data Given :-

Length of pendulum L = 1 metre

Mass of the ball m = 1 kg.

Damping Coefficient b = 0.05

Gravitational Acceleration g = 9.81 m/`s^2`

time span t = 0 to 20 seconds

at t = 0 angular displacement = 0 and angular velocity = 3 rad/sec. [Initial Conditions]

Theory :-

ODEs and ODEs syntax :-

• An ODE is an equation that caontains one independent variable(e.g. time) and one or more derivatives with repect to that independent variable.
• In the time domain, ODEs are initial-value problems, so all the conditions are specified at the initial time t=0.
• Matlab has several different in-built functions for numrical solution of ODEs.
• These solvers can be used with the following syntax:

[output] = function_handle(inputs)

[t,state] = solver (@dstate,tspan,Initial_conditions,options)

• Where State is an array, which represent the solution of ODE at time t.
• Solvers are MATLAB algorithm which are use to solve the ODE. e.g ode45, ode23 etc.
• @dstate handle for function that containing the derivatives.
• tspan is a vector that specifies the interval of the solution.
• Initial condition is a vector of the initial conditions for the system.

Pre-work for defining the given ODE :-

• Given ODE is in the form of Second order differntial equation. [order is represented by the power of the highest order derivative.
• To solve this equation first convert this equation from second order to standard form : `theta_1` = `theta` and `theta_2 = frac {d(theta)}{dt}`
• from this we obtain

`theta_2 = frac {d theta}{dt} = frac {d theta_1}{dt}`

• Putting both value in objective equation :-

`frac {d^2 theta}{dt^2} + frac {b}{m} frac {d theta}{dt} + frac{g}{L} sin theta = 0`

or `frac {d^2 theta_1}{dt^2} + frac {b}{m} frac {d theta_1}{dt} + frac{g}{L} sin theta_1 = 0`

or `frac {d}{dt}(frac {d theta_1}{dt}) + frac {b}{m} frac {d theta_1}{dt} + frac{g}{L} sin theta_1 = 0`

or `frac {d theta_2}{dt} + frac {b}{m} theta_2 + frac{g}{L} sin theta_1 = 0`

or `frac {d theta_2}{dt} = - frac {b}{m} theta_2 - frac{g}{L} sin theta_1`

• No we have two differential equation as follow:-
• First is

`frac {d theta_1}{dt} = theta_2`

• Second equation is : -

`frac {d theta_2}{dt} = - frac {b}{m} theta_2 - frac{g}{L} sin theta_1`

• Now we can define both equations in following form

`frac {d}{dt}[[theta_1],[theta_2]] = [[theta_2],[- frac {b}{m} theta_2 - frac{g}{L} sin theta_1]]`

• L.H.S matrix is known as solution vector and this entire vector is called by `theta`.
• By the value of `theta_1` we get value of displacement and by `theta_2` we get value of velocity for pendulum.

MATLAB Code :-

1. For function :-

First we define function of equation which is to be solve. Name of function is very important because we have to save the file with same name we define in function and same function name is use in calling the function. Here we define the function name as "ode_func" and save the file as "ode_func.m". {Important point}

Code for function : -

% Define equation
% Syntax  function [output] = function_name(Input)
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

% Define solution array theta
    theta1 = theta(1);
    theta2 = theta(2);
% define both solution equations
    dtheta1_dt = theta2;
    dtheta2_dt = -(b/m)*theta2 - (g/1)*sin(theta1);
% define final solution vector
    dtheta_dt = [dtheta1_dt ; dtheta2_dt];
    
end

2. Code for program :-

Following codes are use for function calling and create the animation of pendulum.

clear all   % Clear the Workspace
close all   % Clear all figures and plot
clc         % Clear the Command Window

% Input Data
b = 0.05;   % Damping Coefficient
g = 9.81;   % Acceleration Due to gravity in m/s^2
l = 1;      % Length of Pendulum String in meter
m = 1;      % Mass in Kg

% Initial condition zero for displacement and 3 for velocity
theta_0 = [0 ; 3];

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

% Solve ODE
% ode45 is using to solve the problem.
% ode45 syntax is in the form of as follow.
% [t,y] = ode45(odefun,tspan,y0)
% In this 't' is known as Evaluation points and 'y' is known as solution
% 'odefun' contain the function which is to be solve
% 'tspan' it is known as Interval of integration  
% 'y0' is known as Initial condition {here it is define by theta_0}

[t, result] = ode45(@(t,theta) ode_func(t,theta, b, g, l, m), t_span, theta_0);

% Give command to generate figure in full screen 
figure('units','normalized','outerposition',[0 0 1 1])
% By using subplot command we generate number of plots in a single figure window
% Here by using subplot(6,6,[...]) we divide the figure window in 6 by 6 parts
% In square bracket we allocate the position of particular figure in the figure window 
subplot(6,6,[25,26,27,31,32,33]);
plot(t,result(:,1),'linewidth',1.5,'color','b');  % Time and displacement plot
hold on;
plot(t,result(:,2),'linewidth',1.5,'color','r');  % Time and velocity plot
hold off;
ylabel('Displacement and velocity')
xlabel('Time')
legend('Displacement','Velocity')
title('Change in dispacement and velocity with repect to time')

subplot(6,6,[29,30,35,36]);   % Other subplot and it's position
plot(result(:,1),result(:,2));  % Angular displacement and Angular velocity plot
xlabel('Angular Displacement')
ylabel('Angular Velocity')
title('Relation between Angular Displacement and Angular Velocity')

k=1;  % Constant define to get value for movie in loop
for i=1:length(result(:,1))
    x0 = 0;
    y0 = 0;
    x1 = 1*sin(result(i,1)); % Position of Ball in x-direction
    y1 = -1*cos(result(i,1)); % Position of Ball in y-direction
    subplot(6,6,[2,3,4,5,8,9,10,11,14,15,16,17]); % Subplot for Pendulum
    plot([-1 1],[0 0],'linewidth',2,'color','k'); % Upper line where pendulum is fix
    hold on;
    line([x0 x1],[y0 y1],'color','b','linewidth',1.5); % Pendulum length string
    plot(x1,y1,'-o','markers',20,'markerfacecolor','m'); % Ball of Pendulum
    grid on;
    title('Motion of Pendulum')
    axis([-1.5 1.5 -1.3 0.2]) % Fix axis x and y
    axis equal % Equal distribution of Axis
    pause(0.0003);
    hold off
     
    % getframe command returns the movie frame. 
    M(k)= getframe(gcf); 
    k=k+1;   
end

% Play the movie in matrix M 
movie(M) 
% Assign animation file name and file format
videofile= VideoWriter('Pendulum_all.mp4','MPEG-4' );
% Open the video file for writing
open(videofile)
% Write the video file with the plot frames in the M array.
writeVideo(videofile,M)
% Close the video file after writing.
close(videofile)

After run the program following output is displayed.

Output :-

Animation Link :-

Following animation is result of this program :

https://youtu.be/knTmmsp1o_M

Error during programming :-

Main error arise during programming is 'cdata fields in frames must be the same size'. As shown in figure :

This error is arise due to manual maximize the figure window after run the program. By manual maximization the 'cdata' of frame is get changes and error is occur. To overcome from this problem I use figure command in which I give the command to open figure window in maximization or full screen mode.

END