Week 3 - Solving second order ODEs

Aim:

The main aim is to simulate the transient behavior of the Simple Pendulum and the animation is created.

Theory:

Ordinary Differential Equations:

In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions. The term ordinary is used in contrast with the term partial differential equation which may be with respect to more than one independent variable.

Simple Pendulum:

A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging it back and forth. The mathematics of pendulums are in general quite complicated. Simplifying assumptions can be made, which in the case of a simple pendulum allow the equations of motion to be solved analytically for small-angle oscillations.

Working Principle:

In mechanics and physics, simple harmonic motion is a special type of periodic motion or oscillation where the restoring force is directly proportional to the displacement and acts in the direction opposite to that of displacement. The best example of the Simple Harmonic motion is a Simple Pendulum motion.

A pendulum consists of a massless rod, a Frictionless pivot, a bob. The simple gravity pendulum is an idealized mathematical model of a pendulum. This is a weight (or bob) on the end of a massless cord suspended from a pivot, without friction. When given an initial push, it will swing back and forth at a constant amplitude. Real pendulums are subject to friction and air drag, so the amplitude of their swings declines. Pendulum works on Newton Second Law and a Second Order Ordinary Differential Equations (ODE) is formed that describes the position of the pendulum with respect to time.

Equations Used:

The second order ODE for a simple pendulum:

`(d^2 θ)/(dt^2 )+b/m (dθ)/dt+g/L sin⁡θ=0`……. (1)

Where,

b= Damping Coefficient

m = Mass of the bob (kg)

g = Acceleration due to gravity ( )

L = Length of the rod (m)

To change the ODE from 2^nd order to 1^st order,

Let `θ=θ_1`

Therefore, Equation (1) becomes

`(dθ)/dt=(dθ_1)/dt=θ_2` ……. (2)

`(d^2 θ)/(dt^2 )=(d^2 θ_1)/(dt^2 )=d/dt ((dθ_1)/dt)=(dθ_2)/dt`…… (3)

Substitute Equation (2) and (3) in Equation (1)

Therefore The Equation becomes

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

`(dθ_2)/dt+b/m θ_2+g/L sin⁡θ_1 =0`

`(dθ_2)/dt=-b/m θ_2-g/L sin⁡θ_1`

To arrange it in the matrix form

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

Procedure:

• Initially, the pre-input values are given such as the damping coefficient, Gravitational force of an Earth, Length of a string, and Mass of a ball.

%Initial input Values
b= 0.05;    % Damping Coefficient
g= 9.81;    % Gravity of the Earth in m/s^2
l=1;        % Length of a String in m
m=1;        % Mass of a ball in kg

• The Angular Displacement and Angular Velocity values are taken as inputs for the Initial angles.

%initial condition
Angular_Displacement = input('Enter the Angular Displacement: ');
Angular_Velocity = input('Enter the Angular Velocity: ');
theta_0 = [Angular_Displacement; Angular_Velocity];

• Time is given from 0 to 10 seconds count and 500 frames in between it. It means that the time is equispaced up to 500 values between 0 and 10.

%time Points
time = linspace(0,20,500);

• The ODE function is created using the formula mentioned the ‘Equations Used’

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

• The ODE function is called along with the time and the `theta`

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

• The graph is plotted for the Angular Displacement, Angular Velocity Vs Time.

%Plotting the graph
figure(1)
plot(t,results(:,1))
hold on; 
plot (t, results(:,2));
title('Simple Pendulum Motion');
xlabel('Time (sec)');
ylabel('Angular Displacement and Angular Velocity');
legend('Angular Displacement ($\frac{d\theta_1}{dt}= \theta_2$)', 'Angular Velocity ($\frac{d\theta_2}{dt}= -\frac{b}{m}\theta_2- \frac{g}{l}sin\theta_1$ )', 'Interpreter', 'latex')
saveas(1, 'Simple Pendulum.png');

• To show the simulation of the Pendulum a ‘for’ loop is used to continuously plot the graph showing the pendulum.

%Simulation of Pendulum
count =1; 
for i =1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1));
    y1 = -1*cos(results(i,1));
    figure(2)
    
    plot([-0.5 0.5], [0,0],'linewidth', 5, 'color', 'k');
    line([x0 x1],[y0 y1], 'linewidth', 4, 'color','#A2142F');
    hold on;
    
    plot(x1,y1,'o', 'markers', 20, 'markerfacecolor', 'b');
    grid on;
    title ('Simple Pendulum Motion');
    axis([-1.5 1.5 -1.5 0.5]);
    pause(0.01);
    hold off; 
    M(count) = getframe(gcf);
    text= sprintf("pendulum %d.png",count);
    saveas(2, text);
    count =count +1;
end

• To create the movie, we used VideoWriter to create the animation in a particular format.
• The animation is opened, written, and closed to provide an animation continuously.

%creating movie
movie(M);
videofile = VideoWriter('pendulum_animation.avi', 'uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

Code:

For Function:

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

For Main Code:

% Solving Second-order ODE 
clear
close 
clc

%Initial input Values
b= 0.05;    % Damping Coefficient
g= 9.81;    % Gravity of the Earth in m/s^2
l=1;        % Length of a String in m
m=1;        % Mass of a ball in kg

%initial condition
Angular_Displacement = input('Enter the Angular Displacement: ');
Angular_Velocity = input('Enter the Angular Velocity: ');
theta_0 = [Angular_Displacement; Angular_Velocity];

%time Points
time = linspace(0,20,500);

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

%Plotting the graph
figure(1)
plot(t,results(:,1))
hold on; 
plot (t, results(:,2));
title('Simple Pendulum Motion');
xlabel('Time (sec)');
ylabel('Angular Displacement and Angular Velocity');
legend('Angular Displacement ($\frac{d\theta_1}{dt}= \theta_2$)', 'Angular Velocity ($\frac{d\theta_2}{dt}= -\frac{b}{m}\theta_2- \frac{g}{l}sin\theta_1$ )', 'Interpreter', 'latex')
saveas(1, 'Simple Pendulum.png');

%Simulation of Pendulum
count =1; 
for i =1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1));
    y1 = -1*cos(results(i,1));
    figure(2)
    
    plot([-0.5 0.5], [0,0],'linewidth', 5, 'color', 'k');
    line([x0 x1],[y0 y1], 'linewidth', 4, 'color','#A2142F');
    hold on;
    
    plot(x1,y1,'o', 'markers', 20, 'markerfacecolor', 'b');
    grid on;
    title ('Simple Pendulum Motion');
    axis([-1.5 1.5 -1.5 0.5]);
    pause(0.01);
    hold off; 
    M(count) = getframe(gcf);
    text= sprintf("pendulum %d.png",count);
    saveas(2, text);
    count =count +1;
end

%creating movie
movie(M);
videofile = VideoWriter('pendulum_animation.avi', 'uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);

Terms used:

Ode45:

ode45 is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problem

To solve an ODE, we have to integrate and solve for the variable. Assuming we have an ode du, we would solve for u by integrating du provided some limits of integration. However, rather than direct integration, we can use a numerical solver, by discretizing the equation, that is, we convert the ODE into a linear equation and then solve this equation to approximate the solution.

This function must then be input to ode45 as

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

‘@’(Anonymous Functions):

An anonymous function is a function that is not stored in a program file but is associated with a variable whose data type is function_handle. Anonymous functions can accept multiple inputs and return one output. They can contain only a single executable statement.

The declaration of a function handle is very similar to that of a function declaration. It can be distinguished by the means of the ‘@’ symbol before your value/expression allocation. The syntax for creating a function handle is,

Func_handle=@(input_argument) func_definition

Error:

• The error is faced while creating an animation which is clearly just a manual error.
• The error is rectified by changing the saveas(1, Text) to saveas(2, Text).
• Since there are 2 plots created and the figure number decides which image is needed to be saved.
• The error is overcome by checking the ‘saveas’ term in search documentation in MATLAB.

Result:

The plot is generated between Angular Displacement, Angular Velocity, and time. The animation is generated for the simple pendulum in the required format.

The output is shown below:

The Graph between Angular Displacement, Angular Velocity Vs time is plotted as shown below.

The Animation video is attached below.

or the link to open the Video is Simple Pendulum Animation

Link:

Introduction to ODE: Ordinary Differential Equation Introduction

Pendulum:

• Pendulum- Wikipedia
• Pendulum MathWorks

ode45:

• Understanding ode45
• ode45 MathWorks

Anonymous Function:

• Understanding Anonymous Functions
• MathWorks Anonymous Functions