Solving second order ODEs

Aim:- To write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping

`(d^2theta)/(dt^2) + (b/m).(d theta)/(dt) + g/L.sintheta = 0`

Objective:- To simulate the transient behaviour of a simple pendulum and creating an animation of it\'s motion using Matlab.

Simple Pendulum:-

A pendulum is a weight suspended from a pivot so that it can swing freely. 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 back and forth. The time for one complete cycle, a left swing and a right swing is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude (or the width of the pendulum\'s swing

The simple 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, the bob swings back and forth at a constant amplitude. This is not the case with real pendulums. Unlike simple pendulum, real pendulums are subjected to friction and air drag which reduces the amplitude of their swing with time.

The motion of a pendulum is called a simple harmonic motion and the differential equation which represents the motion of a simple pendulum is :

`(d^2theta)/(dt^2) + g/L . sin theta = 0`

where,

g - acceleration due to gravity in (`m/s^2`),

L - length of the pendulum in (m), and

θ - angular displacement.

As mentioned before, real pendulums experience resistive forces and hence their oscillations die out with time. This effect is called damping. The differential equation for the motion of a simple pendulum with damping is as follows :

`(d^2theta)/(dt^2) + (b/m).(d theta)/(dt) + g/L.sintheta = 0`

where,

θ - angular displacement ,

b - damping coefficient,

m - mass of bob in (kg),

g - acceleration due to gravity in (`m/s^2`),

L - length of the pendulum in (m),

In this project, we will be trying to solve this ODE using a python program and then create an animation of the motion of the pendulum.

Differential equations:-

Differential equation is an equation that relates one or more functions and their derivatives. In applications, the functions represent physical quantities, the derivatives represent their rates of change, and the differential equation defines a relationship between the two. Meaning, one dependent and other independent variable and their relationship.The highest degree to which any equation is subjected, determines it\'s type. To solve this pendulum motion we need to solve a second order differential equation.

ODE45:-

ODE45 is a inbuilt solver in Matlab, that solves differential equations. In this case we provide it with the two inputs arguments and the function which contains the first order differential eqution.

`(d^2theta)/(dt^2) + (b/m).(d theta)/(dt) + g/L.sintheta = 0` 2nd order ODE

The above equation is the second order differential that solves the pendulum motion. We need to convert it into first order equation so that ode45 can solve it.

Let, `theta = theta _1`

`(dθ)/(dt) = (dθ_1)/(dt) = θ_2`

`(d^2θ)/(dt^2)=(d^2θ_1)/(dt^2) = d/dt((dθ_1)/(dt))=d/dt(theta_2)`

`(d^2theta)/(dt^2) + (b/m).(d theta)/(dt) + g/L.sintheta = 0`

`(d^2theta_1)/(dt^2)+(b/m)*(d theta_1)/(dt)+g/L*sintheta_1=0`

`(d theta_2)/(dt)+b/m*theta_2+g/L*sintheta_1=0`

Hence, 2nd order derivative is converted into \'two\' 1st order differential equations.

1. `(d theta_1)/(dt) = theta_2`, and 2. `(d theta_2)/(dt) = -b/m*theta_2 - g/L*sintheta_1`

Code:-

Function Code:-

Function to calculate the ODE\'s

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

Main Code:-

code to calculate the displcement and velocity with respect to time using ode45 solver and above function code also plotting the results frame by frame to get an animation.

t = 20 sec (simulation time);

% Solving second order differential equation and creating an animation
clear all
close all
clc

% Inputs
b = 0.05; % Damping coefficient
g = 9.81; % Gravity
l = 1;    % Length of string
m = 0.1;  % Mass

% Initial conditions
%initial position 0 radians, initial velocity 3 rad/s
theta_0 = [0;3]; % [displacement,velocity]

% Time
%time span to simulate the motion 0-20 sec.
t_span = linspace(0,20,500);

% Solve ODE
%using ode function to obtain results[displacement;velocity]
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%Output plot
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel(\'time\')
ylabel(\'Angular Displacement and Angular velocity\')
grid on
legend(\'Angular Displacement\',\'Angular velocity\')
hold off

%plotting pendulum
figure
ct=1;
for i=1:length(t_span)
    x_fixed=0;
    y_fixed=0;
    x_freeend=l*sin(results(i,1));
    y_freeend=-l*cos(results(i,1));
    plot([x_fixed x_freeend],[y_fixed y_freeend],[-0.5 0.5],[0 0],\'linewidth\',3)
    hold on
    plot(x_freeend,y_freeend,\'marker\',\'o\',\'markerSize\',30,\'markerFacecolor\',\'r\')
    axis([-1.5 1.5 -1.5 0.5])
    grid on
    pause(0.01)
    M(ct)=getframe(gcf);
    ct=ct+1;
    hold off
end

movie(M)
videofile=VideoWriter(\'ode_func.avi\',\'uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Code description:-

1. clear all; close all; and clc; command is used to clear all the variables in the workspace, to close all the figures and plots, clear the command window of all the previous notations and values respectively.

2. Providing all the input values required for the calculations length, mass, gravitational accl., damping coefficient, initial displacement and velocity etc.

3. Defining initial condition theta_0=[t, theta], where, t is an array of time span for which the motion is to be simulated and theta is an array of initial displacement and velocity. Here, time is independent variable whereas displacement and velocity are dependent variables.

4. ode45 is used to solve ODe\'s which takes two variables, one is dependent and another one is independent. But in this case two independent variables are to be obtained, so create a vector of two dependent variables and pass it to the ode solver as a single dependent variable with time as the independent variable.

5. Function is developed where these values are taken as input and parsing this function to the ode45 solver.

6. The function has too many input arguments in order to solve the equation. But as the solver needs only two variables, we have to specify the two variables which are independent and dependent respectively. The syntax for this is to be paid with heavy attention and importance, otherwise silly mistakes will keep on throwing errors.

7. [t, results] contain the results we need to plot with respect to time. Out of which results(1) is the displacement in angle and result(2) is the velocity. We plot both these quantities with respect to time by obtaining the points of the mass at various angular displacement and at the proposed interval of time. We can calculate the formula for the position of mass at different angular displacements from simple pendulum FBD.

8. Depending upon the length of time and angular displacement,large number of frames are obtained. In order to plot all these frames together we use timer and pause every frame in the same plot for 0.01sec.

9. The method to pause the frames one after another is actually not appropriate. So, we create a movie (video file) of this animation with the help of getframe command and counter.

10. Save the videofile in the Matlab directory with the aid of VideoWriter and writevideo commands and get the appropriate animation of the penduluym motion. This is how we can simulate such motions which requires solving of ODEs of different order.

Output Plot:-

The above graph of angular velocity and displacement with respect to 20 sec time shows that the displacement and velocity reduces with the time due to damping coefficient.

If damping coefficient value is set to 0 it will give never-ending motion theoretically, in practical coding, the program skips the value of b (damping coefficient) and runs the program for the defined value of time.

Animation:-

Results File:-

https://drive.google.com/drive/folders/1JzP8Z3jjuZhuKyAnCyukxevvTT29n4u8

Errors:-

1. Undefined function or variable \'theta_o\'.

The error occured because \'theta_o\' was typed instead of \'theta_0\' so matlab threw an error. It was corrected by typing \'theta_0\'.

2. Undefined function or variable \'l\'.

This error occured because uppercase \'L\' was used in main code and smallcase \'l\' was used in function statement. Due to this mistake it caused error in main code, function code, and in ode45 solver also wherever that variable was used.

therefore, pendulum length is corrected by using same variable \'l\' with defined value in main code and function code.