Week 3 - Solving second order ODEs

Aim :- Write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

Equation of motion of simple pendulum –

Aim :- Write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

Equation of motion of simple pendulum –

d²(theta)/dt² + (b/m)*(d(theta)/dt) + ((g/l)*(sin(theta)))

Where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

TERMS :-

• ODE :- Ordinary differential equation, in mathematics, an equationrelating a functionf of one variable to its derivatives. It is used in various application to determine the transient behaviour.

• ODE45 solver :- MATLAB’s standard solver for ordinary differential equations (ODEs) is the function ode45. This function implements a Runge-Kutta method with a variable time step for efficient computation.

dx /dt = f(t, x) , x(t0) = x0

where, t is the independent variable, x is a vector of dependent variables to be found and f(t, x) is a function of t and x.

ODE Simplification -

In order to solve an ODE of nth order, first we have to reduce it into n first order ODEs.

Since, our equation is quadratic, therefore it will be reduced into 2 first order ODEs.

Let, theta = theta1

d(theta)/d(t) = d(theta1)/d(t) = theta2

d²(theta)/dt² = d²(theta1)/dt² = (d/dt)*(d(theta1)/dt) = (d/dt)*(theta2)

1^st single order ODE –

d²(theta)/dt² = d(theta2)/d(t) – eqn.1

2^nd single order ODE –

D(theta2)/dt = (-b/m)*(d(theta1)/dt) –(g/l)*(sin(theta1)) = 0 - eqn.2

Code explaination :-

1. Define the required input values.
2. Define the initial condition for both the ODEs.
3. Theta_0 = [0,3], here 0 shows ‘initial displacement’ & 3 shows ‘initial velocity’.
4. Define time array as t_span and enter the time for pendulum motion. In our case it is 0 - 20sec. with 200 values in between them.
5. In order to solve both first order ODEs for all 200 values of time and to show relation among variables, create a function. Function name – ‘ode_func’.
6. To solve ODE use MATLAB inbuilt solver ODE45 as a function of ‘t &theta’.
7. Extract the values of angular-displacement into another variable ‘a_d’.
8. Plot the angular-displacement & velocity in figure1.
9. Create a for loop to get all the positions (200) of pendulum.

• Plot the pendulum for all the positions.
• Create & save the animation as ‘.avi’ file.

Main program code :-

clear all
close all
clc
%Inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;

% Initial condition
theta_0 = [0;3];

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

% solve ODE

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

%plotting of velocity- displacement relation
figure(1)
plot (t, results(:,1))
hold on 
plot (t, results(:,2))
xlabel ('Time')
ylabel ('plot')
legend('angular-displacement','velocity')


for i = 1:length(a_d)
   
    x0 = 0; 
    y0 = 0;
    x1 = l*sin(a_d(i));
    y1 = -l*cos(a_d(i));
    % Plotting for animation 
    figure(2)
    plot([-1 1], [0 0],'linewidth',3,'color','b')
     hold on
    plot(0,0,'-o','markersize',5,'markeredgecolor','black','markerfacecolor','y')
    plot(x1,y1,'-o','markersize',15,'markeredgecolor','black','markerfacecolor','r')
    hold on
    plot([x0,x1], [y0,y1],'linewidth',1,'color','g')
   hold off 

    axis ([-1.5 1.5 -1.5 1])
    
    pause(0.0001)
    P(ct) = getframe(gcf);
    ct = ct+1;
       
    end

     % To_create & save movie
      movie(P)
      Pendulum = VideoWriter ('simple_Pendulum.avi','uncompressed AVI')
      open (Pendulum)
      writeVideo (Pendulum,P)
      close (Pendulum)

Function code :-

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

Results & discussion :-

1.) Angular-displacement velocity plot (figure 1)

figure - 1

As it is a damped motion, therefore both angular-displacement & velocity is decreasing in each cycle and after sometime it comes to zero.

Animation :- youtube link - https://youtu.be/HODYQcF28IE