Week 3 - Solving second order ODEs

Simple Pendulum:

A simple pendulum has a small-diameter bob and a string that has a very small mass but is strong enough not to stretch appreciably. The linear displacement from pendulum is s, the length of the arc. Also shown are the forces on the bob, which result in a net force of -mg sinθ toward the equilibrium position that is, a restoring force.

  

A simple pendulum is defined to have an object that has a small mass, also known as the pendulum bob, which is suspended from a light wire or string such as shown in figure.

We being by defining the displacement to be the arc length s. We see from figure that the net force on the bob is tangent to the arc and equals -mg sinθ

F = -mg sinθ

Tension in the string exactly cancels the component mg cosθ parallel to the string. This leaves a net restoring force back toward the equilibrium position at θ = 0.

T = mg cosθ

Now if we can show that the restoring force is directly proportional to the displacement, then we have a simple harmonic oscillator. In trying to determine if we have a simple harmonic oscillator, we should note that for small angles (less than about 15) sinθ ≈ θ. Thus for angles less than 15 the restoring force F is,

F ≈ -mgθ

The displacement s is directly proportional to θ. When θ is expressed in radians, the arc length in the circle is related to its radius (L in this instance) by:

s = Lθ

so that,

θ = s/L

For small angles the expression for restoring force is:

F ≈ -(mg/L)s

This expression is of the form,

F = -kx

Where the force constant is given by k = mg/L and the displacement is given by x = s. For angle less than about 15, the restoring force is directly proportional to the displacement, and the simple pendulum is a simple harmonic oscillator.

Using this equation, we can find the period of pendulum for amplitudes less than 15. For the simple pendulum:

T = 2π(m/k)^(1/2) = 2π(m/(mg/L))^(1/2)

Thus

T = 2π(L/g)^(1/2)

For the period of simple pendulum.

Ordinary differential equation:

In mathematics, ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions.

Standard form of second order ODE:

d²y/dx² + b(dy/dx) +cy = r(x)

In engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum, the way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the ball w.r.t. time.

Second order differential equation for a simple pendulum:

Position and velocity of the bob with respect to the time can be determined by solving the second order differential equation of simple pendulum using matlab. In order to solve the second order ODE, we first break the higher order differentials into lower order or simply convert them into first order differentials.

d²θ/dt² + (b/m)(dθ/dt) + (g/L) sinθ = 0

where,

b = damping coefficient

m = mass of the bob (kg)

l = length of the string (m)

g = acceleration of gravity (m/s²)

consider,

θ = θ₁

dθ/dt = dθ₁/dt = θ₂

then,

d²θ/dt² = d²θ₁/dt² = (d/dt)(d(θ₁/dt)) = d(θ₂)

from this we apply the above in the simple pendulum ODE,

d²θ/dt² + (b/m) (dθ₁/dt) + (g/l) sinθ₁ = 0

dθ₂/dt + (b/m)θ₂ + (g/l)sin(θ₁) = 0

dθ₂/dt = -(b/m)(θ₂) – (g/l) sinθ₁

and now forming the matrix,

dθ₁/dt = θ₂

dθ₂/dt = -(b/m)(θ₂) – (g/l)sinθ₁

(d/dt)] = []

Solving the matrix will leave us with theta1 and theta which can be used to solve ODE of simple pendulum.

Function program:

function [dtheta_dt] = ode_function(t,theta,b,m,l,g)

 

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

 

ODE Actual program:

clear all

close all

clc

import math.*

 

%inputs

b = 0.05;

m = 1;

l = 1;

g = 9.81;

 

%initial conditions

theta_0 = [0;3];

 

%time span

time = linspace(0,20,500);

 

%solving ODE

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

 

%plotting

figure(1)

plot(t,result(:,1),'color','r')

hold on

plot(t,result(:,2),'color','b')

xlabel('time')

ylabel('Aungular velocity and Displacement')

legend('change in Displacement','change in Aungular Velocity')

hold off

figure(2)

plot(result(:,1),result(:,2))

xlabel('position')

ylabel('angular velocity')

 

%animation

c = 1

for i = 1:length(result(:,1))

    x0 = 0

    y0 = 0

    x1 = l*sin(result(i,1))

    y1 = -l*cos(result(i,1))

    plot([-1,1],[0,0],'linewidth',4,'color','k')

    hold on

    plot([x0 x1],[y0 y1],'linewidth',2,'color','b')

    hold on

    plot(x1, y1,'o','markers',20,'markerfacecolor','r')

    grid on

    axis([-2 2 -2 2])

    pause(0.03)

    hold off

    M(c) = getframe(gcf)

    c = c + 1

end

 

movie(M)

videofile = VideoWriter('pendulum_ode.avi','UNCOMPRESSED AVI')

open(videofile)

writeVideo(videofile,M)

close(videofile)

 

 

 

Code explanation:

clear all – is used to clear the previously used variables

close all – is used to clear all the previously plotted graph

clc – is used to clear the command window

Initial conditions:

theta_0 = [0;3]

In order to solve second order ODE, we have to be given with 2 initial conditions to proceed with. In this case the initial conditions are in the form of displacement and angular velocity. For the particular problem the initial state of the displacement is given as 0 and the angular velocity of the bob at the initial state is given as 3. The initial condition is represented as a column matrix in the name “theta_0”.

TIME POINTS:

Time = linspace(0,20,500)

Time points are the total time interval for the particular problem for which the pendulum swings. In the given problem, the interval is from 0-20 sec. using the linspace command, 500 different values for t have been created between 0-20 which are stored in the “time” array and are used to plot graph and simulate the movement of the swing.

SOLVING THE ODE:

[t,result] = ode45(@(t,theta) ode_function(t,theta,b,m,l,g), time, theta_0)

The command ODE45 is and inbuild MATLAB ode solver, the solution for the ode is saved in a array [t,result] where ‘t’ represents the time and ‘result’ represents different value for displacement and angular velocity as time changes.

 ode45(@(t,theta)) – the ode has ‘t’ as its independent variable and ‘theta’(an array with theta1 and theta2 in it) as its dependent variable.

ode_function(t,theta,b,m,l,g) – invokes the ‘ode_function’ function which solves the ODE for theta 1 and theta 2 and stores it in a array ‘dthetta_dt’

here ODE45 uses the ‘time’ array to get all 500 different time values and act upon the function with the initial condition being retrieved from the ‘theta_0’

ode_function(t,theta,b,m,l,g):

the function contains the input parameters used to solve the ODE i.e. time, theta, b, m, l, g which are basic and are substituted in the below formulas to determine ‘dtheta_dt’

plotting:

figure(1)

plot(t,result(:,1),’color’,’r’)

hold on

plot(t,result(:,2),’color’,’b’)

xlabel(‘time’)

ylabel(‘Angular velocity and Displacement’)

legend(‘change in displacement’, ‘change in angular velocity’)

hold off

after solving the ODE we now move on to plot the displacement and the angular velocity for different position as time varies in order to visualize how angular velocity and displacement change with respect to each other.

figure() – this command is used to create new figure

plot(t,result(:,1),’color’,’r’) – this plot is used to create a 2D plot of varying displacement for different points, the result(:,1) is used because result matrix comprises of 2 column, first one for displacement and second one for velocity hence in order to access all the values for displacement in that matrix we use it.

hold on – this is used to display different plot in a same graph as without this plot only displays the latest used plot command

 

plot(t,result(:,2),’color’,’b’) – this plot is used to create a 2D plot of varying angular velocity for different time points, the result(:,2) is used because the result matrix comprises of 2 column, first one for displacement and second one for velocity hence in order to access all the values for the angular velocity in that matrix we use it.

xlabel() – is used to label the x axis with the name mentioned within the brackets

ylabel() – is used to label the y axis with the name mentioned within the brackets

hold off – this command is used to hold state to off hence clearing all the plots, axes and their properties.

figure(2)

plot(result(:,1), result(:,2))

xlabel(‘position’)

ylabel(‘angular velocity’)

figure() – this command Is used to create a new figure window

plot(result(:,1), result(:,2)) – creates a 2D plot of varying angular velocity for different displacement. result(:,1) comprises of all the values of displacement while result(:,2) comprises of all angular velocity for different position.

ANIMATION

for I = 1:length(result(:,1))

x0 = 0

y0 = 0

x1 = l*sin(result(i,1))

y1 = -l*cos(result(i,1))

For loops:

For loops are used to repeat a set of code for required number of times according to the program requirements. For loop uses a dynamic variable in its syntax which can be used for increment until it reaches the max value and hence the for loop end

In this case for loop start with i value 1 and ends at size of ‘result’ array

X0 and y0 – represents the origin point or the starting point of the string in the pendulum

X1 and y1 – represents the end point or the point where the bob is connected to the pendulum

l*sin(result(i,1)) – gives the x coordinate of the bob

-l*cos(result(i,1)) – gives the y coordinate of the bob

figure (3)

plot([-1,1],[0,0],’linewidth’,4,’color’,’k’)

hold on

plot([x0 x1],[y0 y1],linewidth,’2’,’color’,’b’)

hold on

plot(x1,y1,’o’,’markers,’20’,’markerfacecolor’,’r’)

grid on

axis ([-2 2 -2 2])

pause (0.03)

hold off

plot([-1,1],[0,0],’linewidth’,4,’color’,’k’) – creates a plot that as a black horizontal line from -2 to 2 which act an wall to which the simple pendulum is tied to.

plot([x0 x1],[y0 y1],linewidth,’2’,’color’,’b’) – creates the string of a simple pendulum, as the line connect the starting point (x0,y0) and the end point (x1,y1)

plot(x1,y1,’o’,’markers,’20’,’markerfacecolor’,’r’) – create a red color bob at the end point of the string of a pendulum with size as mentioned in the command

grid on – adds grid line to the plot

axis([-2 2 -2 2]) – confines both x and y axis

M(c) = getframe(gcf)

c = c + 1

M(c) is a array used to store all the 500 plots that was made during the program, the plot are received with the ‘getframe()’ command

c – a variable which is used to allocate the space inside M so it needs to be incremented regularly as a new plot is formed hence ct = ct + 1 is done

end – this command is used to end the for loop

PLOTS:

Plot for angular velocity, displacement with different time:

  

Plot for angular velocity for varying displacement(position):

 

 Animation video link:

 https://youtu.be/n61Mmlv6Qwg