Week 3 - Solving second order ODEs

Solving second order ODE , Animating the motion of the pendulum.

Objective:

To write a program that solves the ODE representing the equation of motion of a simple pendulum and simulate the pendulum motion.

given, g = gravity in m/s2 = 9.81

L = length of the pendulum in m = 1

m = mass of the ball in kg = 1

b=damping coefficient = 0.05

Procedure:

1. The ODE is solved using ode45 function, which has the nomenclature of

[t, y] = ode45 (fun, trange, init)

where, fun is a function handle or string containing the name of the function that defines the ODE.

trange specifies the time interval over which the ODE will be evaluated.

init contains the initial value for the unknowns.

2. The ode function constructed contains two unknowns,

theta1 = angular displacement,

theta2 = rate of change of angular displacement = angular velocity.

3. It is given that at time t = 0sec, angular velocity = 3rad/sec.

4. Total time of 20 sec is considered for simulating pendulum motion

5. After initializing th given variables, the initial condition for equation is set as theta_0 = [0,5] (displacement and rate of displacement)

6. A vector variable tspan is created to store time interval values using iinspace command.

7. ode45 function is used to solve ODE.

Code:

Main Function

  close all
  clear all
  clc

% Input variables

   b = 0.05;
   g = 9.81;
   l = 1;
   m = 0.1;
   
% Initial condition

  theta_0 = [0; 3];
  
% Time points

  time_span = linspace(0, 20, 500);
 
% Solve ODE
 
  [t, result] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m), time_span, theta_0);
  
% plotting
    
   figure(1);
   
   plot(t, result(:, 1)) % Displacement vs Time graph
   hold on
   plot(t, result(:, 2)) % Velocity vs Time graph
   xlabel("plot")
   ylabel("time")
  
% Animating

  dispacement = result(:, 1);
  ct = 1;
  
  
  
  for i= 1: length(time_span)
    
    x1 =0;
    y1=0;
    
    x2 = l*cos(dispacement(i)-(pi/2)); % claculating X coordinate of line Using Trignometry 
    y2 = l*sin(dispacement(i)-(pi/2)); % claculating Y coordinate of line Using Trignometry  
    
    figure(2)
    
    plot ([x1,x2],[y1,y2],"linewidth", 2)                     % Plotting Pendulum line
    hold on
    plot(x2,y2, '-o', "markersize", 25,'markerfacecolor','r') % Ploting Pendulum bob
    plot([-0.5,0.5],[0,0],"linewidth", 2, 'color', 'b')       % Plotting Support line
    hold off
    axis([-1.5 1.5 -1.5 0.5])                                 % Setting limits
    pause(0.002)
    
   Text = sprintf("pendulum_%d.png",ct);
   saveas(2,Text)
   ct = ct+1;
  
  end

ODE function ode_fun

  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

Output:

1. Ang Velocity vs Time and Displacement Vs Time Graph

Angular velovity at time, 0sec is 3rad/sec

2. Simulation of Pendulum

Link for Youtube video of pendulum simulation : https://youtu.be/7BIFAx2ksc8

Errors faced:

1. Errors faced due to mistakes in function syntax.

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

Error solved by follwing instruction give on command window.

2. Errors faced due to mistakes commited while typing.

thet_0 = [0; 3];

Error solved by proper inspection of code.