Week 3 - Solving second order ODEs

Solving 2^nd order ODE using MATLAB

AIM: To solve a 2^nd order Ordinary Differential Equation that represents the motion of pendulum with damping.

OBJECTIVES:

• To develop an understanding of discretization of higher order ODE’s
• To write a MATLAB code to solve a 2^nd order ODE representing pendulum motion with damping.
• To create an animation in using MATLAB of the pendulum motion for given set of parameters.

THEORY:

In Engineering, ODE is used to describe the transient behaviour 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 these systems.

The following equation represents motion of a simple pendulum with damping, and can be solved to find its angular displacement and angular velocity. It also represents a spring, mass and damper system, solving it would give moment of mass. Thus, it has application in vehicle dynamics to design suspension systems.

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

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b = damping coefficient

Discretizing above equation into a 1^st order ODE-

Let, `ϴ= ϴ_1` and `(dϴ)/dt= ϴ_2`

Therefore, `(d^2 ϴ)/dt= d/dt ((dϴ)/dt)= (dϴ_2)/dt`

Substituting these values in above equation: `(dϴ_2 )/dt + b/m ϴ_2 + g/l sin⁡ϴ_1 =0` which implies:

`(dϴ_2 )/dt= - b/m ϴ_2- g/l sin⁡ϴ_1`

Now we can save the values of ϴ1 and ϴ2 in matrix form:

`ϴ= d/dt (ϴ_1¦ϴ_2 )= (ϴ_2¦(- b/m ϴ_2- g/l sin⁡ϴ_1 ))`

MATLAB CODE:

1) Function containing discretized eqation:

% ODE Equation 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

2) Main Code:

% Solving 2nd order ODE
clc

%Inputs
b = 0.05;    % damping coefficient.
g = 9.81;    % gravity in m/s2 
l = 1;       % length of the pendulum in m
m = 1;       % mass of the ball in kg

%Initial condition
theta_0 = [0;3];

%Time points
t_span = linspace(0,20,200);

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

%Ploting
figure(1);                          %damping plot of displacement & velocity curves
plot(t,results(:,1), 'color','r');
hold on
plot(t,results(:,2), 'color','b');
xlabel('Time [s]')
legend('Displacement', 'Velocity')
title('Damping of Ang. Displacement & Velocity')
hold off

figure(2);                          %Variation of angular velocity with angular displacement 
plot(results(:,1),results(:,2), 'color','b');
xlabel('Displacement [radian]')
ylabel('Angular velocity [rad/s]')

%Iterating to create animation
n = 1;
for i = 1:length(t_span)
    %fixing start point
    x0 = 0;
    y0 = 0;
   
    %defining variation of end point
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    
    %Ploting
    figure(3);                                      %Pendulum position
    plot([-1 1],[0 0],'linewidth',4,'color','r');   %fixing base line
    axis([-2 2 -5 2]);                              %fixing axes
    line([x0 x1],[y0 y1],'linewidth',2,'color','g');%ploting string
    hold on
    plot(x1,y1,'o','markers',15,'markerfacecolor','b','markeredgecolor','b'); %ploting the bob
    hold off
    title('Pendulum Motion');
    legend('Base','String','Bob');
    
    %Storing frames one by one
    z(n) = getframe(gcf);
    n = n + 1;
end

%Creating video file
movie(z)
videofile = VideoWriter('pendulum','Uncompressed avi');
open(videofile)
writeVideo(videofile,z)
close(videofile)
    

CODE EXPLANATION:

1. First, we define a function ‘ode_func’ that will give an array containing values of`(dϴ_1)/dt`and`(dϴ_2)/dt` as the output.
2. Now to star the main code by initializing the input values of damping coefficient, gravitational acceleration, length of string and mass of the bob.
3. Next we store the initial values of angular displacement and angular velocity of 0 and 3 respectively as an array in the variable ‘theta_0’.
4. We also need to need to create an array of time values (‘t_span’) between 0-20 seconds.
5. The most important part of the code is using ‘ode45’ method to solve the required ODE.

Its syntax is as follows:

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

Where, [t, results] stores the output of ode45 in two variables ‘t’ & ‘results’

@(t,theta) independent and dependent variables

ode_func(t,theta,b,g,l,m) function which is needed to be solved

t_span values of independent variable ‘t’

theta_0 values of dependent variable ‘ϴ’

6. Plot the graphs showing damping of angular displacement and velocity with respect to time. And variation of angular velocity with respect to angular displacement.

7. Now to create an animation we need to save various positions of bob to generate frames. For that purpose, we use a for loop to iterate through the process.
8. To start with the loop, we fix the start point (x0, y0) and then define the motion of end point (x1, y1) which will be an array.
9. To make the plot look like a pendulum we must plot a base line and the string of the pendulum separately along with the bob at the end point of the string.
10. After finishing up with plot to store the frames we use the ‘getframe’ method paired with a variable to store them. This is done in the loop itself. Which also finishes up the loop.
11. Lastly to create a video out of the stored frames we use the ‘movie’ and ‘’VideoWriter’’ methods. Which completes our code.

CODE OUTPUT:

1) pendulum motion:

2)Plot of angular displacement vs anguar velocity:

3) Damping plot:

4) Animation:

https://youtu.be/B7djDp9Y5XQ

ERRORS OCCURED:

1) In the function cntaining the equation 'ode_func' I made a typo in'sin(theta1)' accidently writting it as sind(theta1), which gave a completely different damping plot than expected since th 'sind()' method takes the input in degrees than in radians. I had to correct it, but it was quite a bug.

The incorrect plot:

RESULTS:

As we can see the plots are not verry smooth, we can increase the smoothness by inicreasenig the no. of time points. This can be done by increasing the no. of splits in the 'linspace' method while declaring the t_span array.

CONCLUSION:

Hence we have successfully learned how to use discretization and implement it to write a code in MATLAB to solve the ODE more efficiently. Also we created an animation of the pendulum motion which proves the correctness of our code.