Week 3 - Solving second order ODEs

Aim: Simulating a Pendulum Motion by using 2nd Order ODE.

Governing Formula:

d^2theta/2 + b/m(dtheta/dt) + g/l = 0

Main Objective of Project:

To creat a animation of swinging pendulum by using 2nd order ODE in Matlab.

Procedure:

• First we will divide 2nd order formula in two arrays
  1. Let's assume: Theta = Theta1
  2. dtheta/dt = dtheta1/dt = theta2 ------ eq-1
  3. d^2theta1/dt = d/dt(dtheta1/dt) = dtheta2/dt --------- eq-2
  Applying above equations on ODE formula
  Therefore, dtheta2/dt +b/m(theta2) + g/l(sin(theta1)) = 0
  
  Dividing above eq into to arrays:
  dtheta2/dt = -b/m(theta2) - g/l(sine(theta1))
• We wil start writing program by typing 'clear All', 'close All', 'clc' to reset command window and workspace
• Now we will type in the given values i.e.
  Length = 1 meter
  Mass = 1kg
  Dampning Coeff. = 0.05
  Gravity = 9.81 m/s^2
• Then we will type in initial condition i.e.
  Initial displacement = 0
  Initial velocity = 3 rad/sec
• Now we will set the duration of the animation for 20sec and we will divide the time span in 500 division using 'linspace' commande and save it in variable 't_span'
• To solve ODE we will use inbuilt solver called 'ode45' which will difine time and theta as 't' and 'theta' to our function 'ode_func'. This wil give as two results i.e. 't' and 'result'.
• Ode function definition:
  1. First we will define name of the function as 'ode_func' and type all the variables required in the formula as function input.
  2. First we will define 'theta1'(displacement) and theta2 (velocity).
  3. Now we will use the equation whice we derived earlier.
  dtheta1_dt = theta2
  dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1)
  4. then we will merge them in single output i.e.
  dtheat_dt = [dtheta1_dt; dtheta2_dt]
  5. We will end this function by typing 'end'
• For Loop:
  1. result of the above equation has to columns i.e. displacement and velocity. To loop from 1st column of result first we need to creat a variable where we can save the whole column.
  2. 'theta_1' will have displacement only.
  3. we will loop through 'theta_1' and on each loop we wil collect data as 'THETA'.
  4. Then we will give the starting point (hinge) of pendulum.
  5. Then we will use sin and cos theta to plot the final point af the other side of pendulum.
  6. Plot the line on thegraph by using 'plot' command and put both x and y walue for both condition as arraies.
  7. pause time by 0.03sec to slow down the animation.
  8. Set the graph limit by using 'axis' command.
  9. We will capture frames by using command 'M(ct) = getframe(gcg)' and count the frame in 'ct' variable which we created earlier.
  10. End the loop by typing 'end'.
• Now we will create a movie by difining 'movie(M)'
• 'VideWriter' will create a new folder in avi formate named as 'pendulum_animation.avi' and save it in variable 'videofile'.
• open file by using command 'open'
• write movie by using command 'writeVideo'
• Then close the file by command 'close'

Conclusion:

• We will get the animation showing the swinging of pendulum.
• Problem faced during coding:
  1. Was not able to create a for-loop for 1st column of 'result' variable.
  solution: Created a variable to collect all the content of column 1 in 'result' variable.
  2. Could not animat pendulum properly.

Youtube link for animation:

https://youtu.be/4PPqOpCA7BY

clear all
close all
clc

%given values
b=0.05;
g=9.81;
l=1;
m=1;

% initial condition
theta_0 = [0;3];

%time point
t_span = linspace(0,20,500);

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

ct = 1;
theta_1 = result(:,1);
for i = l:length(theta_1)
    THETA = theta_1(i);
    
    % Initial point
    x0 = 0;
    y0 = 0;
    
    % Final point
    
    x1 = l*cos(THETA)
    y1 = l*sin(THETA)
    
    %Plotting
    plot([x0 x1],[y0 y1],'linewidth',3)
    pause(0.03)
    axis([-2 2 -2 2])
    M(ct) = getframe(gcf);
    ct=ct+1;
    
end

movie(M)
videofile = VideoWriter('pendulam_animation', 'uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

%------------------ Funcition Definition -------------------%

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