Week 3 - Solving second order ODEs

AIM: To Solve the Ordinary Differential Equation of second order and plot displacement and velocity using Matlab/Octave.

Objective:

To create a second order differential equation using matlab/octave.

To use 'ode45' solver to solve the ordinary differential Equation.

To plot the different position of the pendulum and later animate them.

To plot the displacement and velocity of the pendulum in timeline.

Introduction:

A linear differential equation is a differential equation that is defined by a linear polynomial in the unknown function and its derivatives, that is an equation of the form:

where;

, and are arbitrary differentiable functions that do not need to be linear,

and are the successive derivatives of the unknown function y of the variable x.

Among ordinary differential equations, linear differential equations play a prominent role for several reasons. Most elementary and special functions that are encountered in physics and applied mathematics are solutions of linear differential equations (i.e. Holonomic function).

When physical phenomena are modeled with non-linear equations, they are generally approximated by linear differential equations for an easier solution. The few non-linear ODEs that can be solved explicitly are generally solved by transforming the equation into an equivalent linear ODE (for example Riccati equation).

Here, We have taken ordinary differntial equation of second order, for the pendulum under simple harmonic motion with damping.

Where, b = damping coefficient

m = mass of the ball

l = length of the pendulum

g = gravity

Here, the single differential equation of second order is converted into two differntial equation of first order.

Now, equation no.(1) and (2) are the main function to be solved using matlab.

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.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

Explanation/Steps to follow:

Here, I am using Octave Software, but the program's code is similar to the Matlab.

There is only difference in coding at animation of the pendulum.

NOTE: Use Octave 5.2.0 Version, as ode45 solver is not available in earlier version.

So, before writing the main program we will write our ordinary differential Equation in our program.

Here, we will define the function, as we have evaluated in the derivation.

Remember that the function name and the file name should be same as we are calling the function which is defined by it's file name only.

Also, the function should give the output in terms of our derived equation.

After writing the function definition, we will go to new script and write the code for our main program.

Now, before writing up the main program, try to use comments as much as possible as it is easier for the reader to get what is written in the code, also it is easier for us, if we open after few months and it becomes easier to gasps the code and their meaning.

To create the program initially; clear all, close and clc command should be used.

It is not mandatory but it is a good to clear the command window and workspace and the different windows of the software, so that you program looks clean and clear also closing the different unnessary windows helps to run the software smoothly.

Start with writing the starting points of the links and ending points of the links, so that you will not have a writing all the things at a same time.

It is always a good practice to write the program in small sections and checking for the errors and afterwards combine them together.

Now, after that we have to define the inputs which are fixed and needed for calculation in our program.

After defining inputs, we have to define some initial conditions such as initial velocity and initial displacement of the pendulum.

After that we have to define time for which we have to do analysis.

To define time, I have used the variable t_span.Here, the t_span variable, I have given is from 0 sec to 20 seconds and I wanted 200 points in between these time span, that's why i have used linspace command.

After that, we need to call the function and need to solve our ODE function.

So, to do that, I have created a array of two variables which is time and results i.e. [t, results]

Now, I have used ode45, it is a solver which is used to solve our equation.

The synatx of the ode45 solver is ode45(@(a,b) ode_func(x,y,z) , t_span, initial);

where, @(a,b) is the variable which is changed while solving the equation.

ode_func is our ODE function which is to be solved, and which we have defined earlier.

then comes our time span and after that the initial conditions.

Now, we will run the program and check for the program if it has any problems.

Now, we will use the plot command and plot the displacement and velocity, if it looks good as we wanted then we will move for the animation part of the pendulum.

Now, when we get the plot as required, we also need to animate the pendulum, to check that our solution is 100% correctly.

Note:Here, i am using Octave software so, there are some steps which are different from the matlab programming.

Now, we are using for loop for our animation purpose. Before for loop, i have used 'ct' which is counter variable, it is used for iteration purpose also it is used to check if the program is running perfectly or not, also it is used in different purposes such as plotting different graph in new window etc.

Now, In for loop the range of 'i' is from 1 to 200 (length(results(:,1)).

Now, after that we are required to write the formula for the position of the pendulum.

Now, we know that one end of the pendulum is fixed which is (x0,y0).

The Other end of the pendulum is simply taken as the cos and sin components and in the value of theta, it is taken as the results of the ode function. Here, neguative y direction makes the pendulum to oscillate in downward direction.

And then the plot command is used to plot the positions and after that saveas function is used to save a single image of the position, here point to be noted that matlab uses some other function to plot the image and make a animation from it directly.

Given Data:

Length of pendulum =1 metre

mass =1 kg

Damping co-efficient (b) =0.05

Gravity =9.81 m/s2

Program:

1.Program to define the Ode Function.

%Here, Odinary Differential Equation is defined. 
%ode_func is the variable which is a function of time,angle, damping coefficient....
%...length and mass of the ball at the end of pendulum.
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);
%It returns the output as dtheta1/dt and dtheta2/dt
dtheta_dt = [dtheta1_dt,dtheta2_dt];

end

2.Program to Solve the above Differential Equation.

%Program to solve second order differntial equation

clear all
close all
clc

%Inputs
b = 0.05; %b is the damping co-efficient
g = 9.81; %gravity
l = 1; %length of the pendulum
m = 0.1; %mass 

%Initial Conditions
theta_0 = [0;5]; %Here, it takes initial displacement as zero and intial velocity as 5.

%Time span upto which we need to analysis the pendulum
t_span = linspace(0,20,200);

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

%Saving the image file and making animation.
ct = 1
for i = 1:length(results(:,1));
     %Here, Pendulum position's formula is written
     % 1st point of pendulum fixed at (0,0)
     x0 = 0;
     y0 = 0;
     %Second Position of pendulum changing as per given formula:
     x1 = l*sin(results(i,1));
     y1 = -l*cos(results(i,1));
     
     %Plotting the position of the pendulum's position
     plot([-1.5 1.5],[0 0], 'linewidth',4,'color', 'k');#Pendulum Base.
     hold on
     plot([x0 x1],[y0 y1],'linewidth',4,'color', 'g');#pendulum Position
     hold on
     plot(x1,y1,'o','markersize', 15, 'color','r'); #Pendulum Ball at end
     axis([-1.5 1.5 -1.5 0.5]);
     hold off
     
     %To save the image files of theposition of pendulum 
     Text = sprintf("pendulum%d.png",ct);
     saveas(gca, Text);
     ct = ct + 1;
end

%Plotting Displacement and Velocity in Magnitude Vs Time
figure(101);
plot(results(:,1), 'linewidth', 4);
hold on
plot(results(:,2), 'linewidth', 4);
xlabel('Magnitude', 'color', 'r')
ylabel('Plot', 'color', 'r')
legend('Displacement','Velocity')

Output:

The Displacement and velocity with respect to time is shown as in below plot:

Animation:

The Different Position occupied by the pendulum during it's motion at defined time span is shown by this animation: