Week 3 - Solving second order ODEs

OBJECTIVE :

To write a program to solve the ODE equation which is required to simulate the motion of a simple pendulum with damping.

MATLAB CODE :

l = 3;
% Initial Condition 
% Angular Displacement - 0
% Angular Velocity - 3 rad/sec
ct = 1;
Angular_Displacement_and_Velocity = [0;3];

% Simulation time span

Motion = linspace(0,20,80)

% ODE EQUATION

[t,results] = ode45(@(t,theta) pendulum_func(t,theta),Motion,Angular_Displacement_and_Velocity)


theta = results(:,1);

for i = 1 : length(Motion)
    X0 = 0;
    Y0 = 0;
    X1 = l*sin(theta(i))
    Y1 = -l*cos(theta(i))
    plot([-4,4],[0,0])
    ylabel('Plot')
    xlabel('time')
    hold on
    plot([X0 X1],[Y0 Y1],'linewidth',3,'Color','k')
    plot(X1,Y1,'ko','MarkerSize',20,'MarkerFaceColor','m','MarkerEdgeColor','m')
    pause(0.003)
    axis([-4 4 -4 4])
    M(ct) = getframe(gcf);
    ct = ct + 1;
    hold off
    
end
%%matlab creates a movie for this code by using "Movie" Command
movie(M)  
Simulation_pendu = VideoWriter('Simulation_Pendulum.avi','Uncompressed AVI');
open(Simulation_pendu)
writeVideo(Simulation_pendu,M) % you need to write the video file in M
close(Simulation_pendu)

• The function Line for ODE 45 solver is by the syntax [t,y] = ode45(@(t0,y0) func_name(t0,y0), tspan , y0)

Where,

t,y = Output arguments ,where ode 4 solver gives you the value of Y by derivating for different set of t values. t is the scalar quantiy which is an independent variable and Y - is the vector quantity which is a state / dependent variable

t0 , y0 - Input arguments / Variables to make understand matlab that we are giving inputs for t0,y0 so that these variables are outputed after the function code gets executed.

Function_name - This is the function template we call part of the solver

tspan - This is the time steps, time span, or the steps the matlab considers part of integration.

y0 - Intial conditions which is angular displacement and angular velocity here in this condition

• So primarily we provide the initial conditions for angular displacement and angular velocity.
• Now we provide the time span under the motion variable from o to 20 seconds considering 80 steps / values
• [t,results] = ode45(@(t,theta) pendulum_func(t,theta),Motion,Angular_Displacement_and_Velocity) Now i've given the ode solver function code line
• Under the pendulum_func, i've attached the function code line wherein i've provided values for b (damping coefficient) , g (graviy) , L (Length of the pendulum rod ) , m (Mass of the ball)
• Now we primarily do a function code line to solve a polynomial function it can be first order , second order ot nth order polynomial.
• So we give dheta_dt as the output argument and we call the function template and provide input values for the variables specified.

  function[dtheta_dt] = pendulum_func(t,theta,b,g,L,m)
  
  % Parameters for calculating the second order equation for simple pendulum
  b = 0.05; % Damping Coefficient
  m = 1; % Mass of the ball in kg
  g = 9.81; % Acceleration due to Gravity in m/s^2
  L = 1; % Length of the pendulum in 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
  

• theta1 = theta(1) means theta(1) from derivation is represented as theta1 as a variable similarly for theta2
• We have a second order polynomial, so inorder to solve the equation we need to derivate or solve it to a series of first order polynomial. Here we need two first order polynomials which is dtheta1/dt and detheta2 / dt
• So,`(d theta)/dt = [ (D theta_1) / (dt) ; (D theta_2) / (dt)]`
• So we are going to obtain a column vector for dheta_dt
• Upon running the program we get a 80 rows X 1 column for t ; and 2 outputs of the same order for results
• Now to get the pendulum oscillation with a fixed end, i need to know the position of the pendulum for a sequence of 80 steps and i'm using a loop counter inside a for loop as 'i' and im counting to the length of "Motion"
• In order to find the respective theta value, i'm using the coordinates formula , (l cos`theta`, l sin `theta`)
• Here , when i find the values using the right angled triangle logic for sin theta and i get the values for x & y.
• For y = - l cos `theta` because, the pendulum is below the fixed plane and not above, hence plotting accordingly.
• Sintheta = opposite / hypotenuse ; costheta = adjacent/hypotenuse.
• we get the values using right angled triangle logic.
• Now we plot the fixed plane, X0 ,Y0 & X1 Y1 , and give the respective labels.
• Nowe we hold on to fix a marker icon at X1 & Y1 and pause for 0.003 seconds and we get the frames using the counter and initialise it to variable M.
• Later we write the video file and obtain the animation.

ERRORS :

Some of the errors obtained were inconsistent matrices which i initialised by using double loop counters, later i realised a single loop counter is enough to plot the position and orientation of pendulum.

I mistook that i need to plot time vs position part of the plot and thats how i get the output simulation. I later understood, position vector which is results(:,1) is the only input needed to get the simulation of the pendulum.

I made small errors when initialising the solver , i tried doing it without the '@' , later i realised those are input arguments whose order is followed exactly inside the function template code.

I ploted the fixed plan axis wrong, i gave values which are not part of the scale. Later rectified it by giving appropriate coordinates.

OUTPUT :

https://youtu.be/rJmE4P3E758