Week 3 - Solving second order ODEs Using MATLAB

AIM:

To write a program that solves 2^nd order ODE (ordinary differential equation) of the simple pendulum.

GOVERNING EQUATION:

• ODE for simple pendulum:

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

`(d^2 θ)/(dt^2) +b/m (dθ)/dt + g/l sinθ=0`

Here,

b = Damping coefficient,

g = gravity (in m/sec²),

l = length of the pendulum (in m),

m = mass of the ball (in kg).

The way the pendulum moves depends on the Newton’s second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" with respect to time.

OBJECTIVE OF THE PROJECT:

1. To write the program in MATLAB to simulate the pendulum motion for the following conditions-

• l = 1 m.
• m = 1 kg.
• b = 0.05.
• g = 9.81 m/sec²

2. Simulate the motion between 0-20 sec, for angular displacement=0, angular velocity=3 rad/sec at time t=0. Get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum.

BASIC THEORY:

Second-order differential equations play a central role in the physical sciences. They are found, for example, in laws describing mechanical systems, wave motion, electric currents and quantum phenomena.

To take a simple case, consider a particle of mass m that moves in one dimension along the x-axis. At any given time t, the particle’s position is x(t), and its velocity and acceleration are given by the derivatives dx/dt and d²x/dt² .

A very important law for its acceleration: Newton’s second law tells us that

Mass × acceleration = force,

Which implies that,

Where F is the force acting on the particle. The force need not be constant, and may vary with the position x or the velocity dx/dt of the particle. So, depending on the precise details, we get a second-order differential equation for x as a function of t, and the solution of this equation tells us how the particle can move.

The system known as a simple harmonic oscillator provides a good example.

BODY OF THE CONTENT:

• MATLAB Programming language is used in this program.
• Main Code to solve 2^nd order ODE of the simple pendulum:

% A program to solve 2nd order ODE of the simple pendulum
close all
clear all
clc

% Inputs - basic parameters 
% Length of pendulum (in m)
l = 1;

% Mass of the ball (in kg)
m = 1;

% Damping coefficient 
b = 0.05;

% Gravity (in m/sec^2)
g = 9.81;

% Initial condition - displacement & angular velocity
theta_0 = [0;3];

% Time interval(in sec)
t_span = linspace(0,20,500);

% Variation in angular displacement(theta1) and angular velocity (theta2)
% ode45- Standard solver for ordinary differential equation.
[t_span,result] = ode45(@(t,theta) ode_func(t,theta,l,m,b,g),t_span, theta_0);
figure(1)
% Plot of angular displacement(theta1) Vs time
plot(t_span,result(:,1),'color','b')
hold on
% Plot of angular velocity(theta2) Vs time
plot(t_span,result(:,2),'color','r')
xlabel('Time')
ylabel('Angular Displacement & Angular velocity')
hold off 

% Angular displacement stored in theta1
theta1 = result(:,1);
% Count tracker
ct = 1;

% For-loop for the animation of simple pendulum
% Length(theta1) have element equal to the element in the result(:,1)(500)
for i = 1:length(theta1)
    % Displacement array
    theta_v = theta1(i);
   
    % Fixed end coordinates of the pendulum rod
    x0 = 0;
    y0 = 0;
    
    % Movable end coordinates of the pendulum rod
    x1 = l*sin(theta_v);
    y1 = -l*cos(theta_v);
    
    % Plot of Displacement Vs time of pendulum
    figure(2)
    plot([-0.3 0.3],[0 0],0,0,'.r',[x0 x1],[y0 y1],'-k',x1,y1,'om','markerfacecolor','m','MarkerSize',17,'linewidth',3)
    axis([-2 2 -2 2])
    pause(0.02)
    xlabel('Displacement')
    ylabel('Time')
    % getframe command to capture every position of pendulum with respect to time in form of frame
    M(ct) = getframe(gcf);
    ct = ct+1;
    % end of the for-loop
end
% movie command to form a movie with captured frames
movie(M)
% Properties of videofile, video name and video compression method
videofile = VideoWriter('Pendulum','Uncompressed AVI');
% To open videofile
open(videofile)
% To apply all the videofile properties to the formed movie/animation 
writeVideo(videofile,M)
% To close formed movie/animation
close(videofile)

• Function for the ODE:

% function to solve simplified 2nd order ODE of simple pendulum
function [dtheta_dt] = ode_func(t,theta,l,m,b,g);

theta1 = theta(1);
theta2 = theta(2);

% Convert 2nd order ODE into 1st order ODE
% d^2theta_dt^2 + (b/m)*dtheta_dt + (g/l)*sind(theta)
% theta = theta1
% dtheta1_dt = theta2 
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 -(g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];

end

Step wise explanation of the code:

1. First write down the basic commands such as close all, clear all, clc.
2. Provide the given parameters- damping coefficient, length of the pendulum, gravity and mass of the ball.
3. Then the given initial conditions that is- angular displacement = 0 and angular velocity = 3 red/sec as theta_0 at time t = 0.
4. Providing the time interval of 0-20 sec as t_span (time span) in 500 equal elements.
5. Program the function code of the simplified 2^nd order ODE (governing equation) in following steps-

• First let’s consider that theta equals to theta1 that is angular displacement equals to theta1 and first derivative of the theta equals to theta2.that is,

`theta=theta_1` & `( d theta_1)/dt=theta_2`

• After substituting considered values (theta1 and theta2) the 2^nd order equation convert into 1^st order equation as

`(d theta_2)/dt+ b/m×theta_2+ g/l×sin⁡(theta_1)=0`

• Final matrix will be as follows-

`[[(dθ_1)/dt],[(dθ_2)/dt]]= [[θ_2],[-b/m θ_2- g/l sin⁡θ_1]]`

• The values of dtheta1/dt and dtheta2/dt form dtheta/dt array as a matrix of 2x1 order.

`[[(d theta_1)/dt],[(d theta_2)/dt]]`

• Consider dtheta/dt array as a function of- ode_func(t, theta, l, b, m, g).

6. For calculation of angular displacement and angular velocity with respect to given time interval ode45 - standard ODE solver is used.

• About ode45 - Ode45 is a function in matlab that solves ordinary or partial differential equations using a numerical method. The numerical method used here is the 4th order Runge Kutta method.
• ode45 working - The function ode45 takes 3 inputs. The 3 inputs are the function to be integrated, the time span and the initial conditions.
• The RK method is a first order ODE solver and hence, so is ode45. Thus, any ODE to be integrated by ode45 must be provided as a set of first order ODEs.

7. Plot angular displacement and angular velocity with respect to given time interval.
8. For getting every position of pendulum with respect to the time interval to form animation following step are used-

• Store resultant angular displacement values with respect to time as theta1 in form of array.
• For animation, let theta_v = theta1 as an array to theta with respect to time.
• x and y coordinates as a function of theta_v.
• Plot coordinates – x0,y0 (0,0) ,x1 and y1 .

9. Commands for animation that is- getframe, movie, videofile, writeVideo, open and close.

PLOTS:

• Plot of Angular Displacement and Angular Velocity Vs Time

• Plot of Displacement Vs Time

• Screenshot of the plot frame for animation

RESULTS:

• Main Code Final values

ERRORS:

• Because of taking mean position coordinates as (0,-L) length of pendulum is so short. so i changed mean position coordinates as (0,0).

• Now due to the hold on command between both the plots all the positions get plotted with respect to time but all positions remain plotted and made plot dense.

• Final plots, by mentioning all the coordinates in a single plot.

CONCLUSION:

ode45 function makes ODE calculation easy and pendulum motion animation work probably. Conversion of 2nd order ODE into 1st order ODE makes calculation easy and also required as one of input of ode45.

YouTube link for the Animation Video:

https://youtu.be/iiuEYzC4STc