Second Order Differential

AIM: TO write a program that will simulate the motion of a simple pendulum

GOVERNING EQUATIONS:   

     A simple pendulum ,in its simple form consists of a heavy bob suspended at the end of a light inextensible, flexible string and the other end of the string is fixed at O. the pendulum is in equilibrium when the bob is at A. If the bob is moved from B to C and released, it will start vibrating between the positions B and C with A as mean position.

Consider the equilibrium of the system at C,the weight mg of the bob can be resolved into two components i.e., mgcos`theta` and mg sin`theta` at right angles to each other. The mgcos`theta` will act along the string and balances the tension in the string. The component mg sin`theta` being unbalanced will give rise to an acceleration in the direction of CA.

              In Engineering, ODE Order differential equations are used to describe the transient behavior of a system. Here we have taken the system as simple pendulum. The way the pendulum moves depends on Newton’s second law of motion. When this law is written down, second order ordinary differential equations describe the position of the ball w.r.t time.

`frac{d^2theta}{dt^2}=(frac{b}{m})frac{d theta}{d t}+(frac{g}{L})sin theta=0`

Where b=damping coefficient

            m=mass of the ball, kg

              g=gravity in m/s2

               L=length of the pendulum in m

 To solve this ODE, we have to integrate and solve for the variable. we use numeric solver by discretizing the equation that is we convert the ODE into a linear equation and then solve this equation to approximate the solution.

                       `theta=theta1`

                        `frac{d theta}{d t}=frac{d theta1}{d t}=theta2` ------------1

                          `frac{d^2theta}{d t^2}=frac{d^2 theta1}{d t^2}=frac{d (frac{d theta1}{d t})}{d t}=frac{d (theta2)}{d t}`--------2

Now substitute equation 1 and 2 in ODE of simple pendulum, `frac{d^2theta}{dt^2}=(frac{b}{m})frac{d theta}{d t}+(frac{g}{L})sin theta=0`

`frac{d^2theta1}{dt^2}=(frac{b}{m})frac{d theta1}{d t}+(frac{g}{L})sin theta1=0`

`frac{d theta1}{d t}=theta2`----3

`frac{d theta2}{d t}=-(frac{b}{m})theta2-(frac{g}{L})sintheta1`----4

when we solve the above linear equations 3 and 4, we get output results as below

`frac{d[[theta1],[theta2]]}{d t}=[[theta2],[-(frac{b}{m})theta2-(frac{g}{L})sintheta1]]`                

  OBJECTIVE:

1. To write a program that solve the ODE of the pendulum
2. To do simulation the motion of a pendulum between 0- 20 seconds
3. Upload the animation file.

Inputs:

L=1m

g=9.81m/s2

b=0.05

m=1kg

CODING:

The ODE45 function we use to  solve ordinary differential equation is ode_func(). The variable we have taken here are time,t and angular displacement, theta. This function takes three inputs ,they are function which is to integrate , time span and the initial conditions.

The function is ode_func()

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

The time span is 0 to 20 sec and in between 500 discrete values are taken so that at each value the angular displacement and velocity is calculated using ODE. The time span would be the integration limits which will be represented as a vector of 2 values the initial time and final time with a spacing using the linspace function. The output assigned to a variable t_span.

     t_span=linspace(0,20,500)

Initial conditions is taken as column vector of two values ie.,first one is angular displacement  value taken is zero and the second one is angular  velocity i.e., the pendulum moves with an angular velocity of 3m/s at zero angular displacement.

      theta_0=[0;3]

The function ode45 is in the form of an anonymous function,

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

 The output of this ode45 is a vector of two columns and 500 rows. The first column is the angular displacement and the second column is the corresponding angular velocity at each discrete value of time between the time span of 20 seconds. The obtained matrix is assigned to the variable 'results' .

  The plot we get by taking the output values  of ODE45 function, 'results' with respect to time t. The blue curve represents angular displacement aginst time stating at inital value theta=0. The orange curve represents angular velocity against time at initial value ,v=3m/s.

                        plot(t,results(:,1)) ---values of first column of results matrix against time are taken
                        hold on                ---- hold on command is avoid overwrite the output results.
                        plot(t,results(:,2))-------values of second column of results matrix against time are taken
                        xlabel('Time')
                        ylabel('plot')

SIMULATION: 

         To do the simulation of the simple pendulum, we take for loop with a condition of i=1:results(:,1). here the variable results is the output of ODE45 function, which is a matrix of 2 columns and 500 rows.The first column values, angular displacement of the bob are taken for iteration of the loop. So at each value of angular displacement, x1 and y1  are calculated while x0 and y0 are taken as (0,0) and is fixed. For each iteration plot is captured and displayed in a single frame by using command getframe(gcf).

Marker command is used to create the bob of size 15 diameter and color red.

ct command is used to for each frame iteration.

clear all
close all
clc

b=0.05;
g=9.81;
l=1;
m=0.1;

theta_0=[0;3];
t_span=linspace(0,20,500);

[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('plot')
 
figure(2)
ct=1
for i=1:length(results(:,1))
     x0=0;
     y0=0;
     
     x1=l*sin(results(i));
     y1=-l*cos(results(i));
  
     plot([x0 x1],[y0 y1],'linewidth',3)
     hold on
     plot(x1,y1,'marker', 'o','markersize',15,'markerfacecolor','r');
     hold off
     axis([-1 1 -1 1]);
     M(ct)=getframe(gcf)
     pause(0.03)
     
     ct=ct+1
end
movie(M)
videofile=VideoWriter('Simple Pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Errors: ODE function is given inside the for loop,got this error.

youtube link: https://youtu.be/VyheJuAXeW8

CONCLUSION: Using  ODE functions we can solve many such system. These methods are of greater importance when we realise that computing machines are now available which reduces numerical work considerably. ODE45 function used in this program uses Ranga-Kutta Method.Many differential equations arising from physical problems are linear but have variable coefficients and do not permit a general solution in terms of known functions. such equation can be be solved by Numerical methods.