Week 3 - Solving second order ODEs

Aim: A program that solves second order ODE.

Theory:

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time. This ODE represents the equation of motion of a simple pendulum with damping.

In the above equation,

g = gravity in m/s2, L = length of the pendulum in meter, m = mass of the ball in kg,b=damping coefficient.

Solution:

We solve above second order ODE by using Numerical Solution technique,2nd Order ODE is converted to two 1st order ODE which are as follows.Intial values are given as `theta(t=0)=0; (d/dt(theta)) at(t=0) =3`

`[ d/dt(theta1) ; d/dt(theta2)]=[ theta1 ; -g/L*sin(theta1)-b/m*d/dt*(theta1)`

Above Vector containing two 1st order ODE along with vector containing intial values with time span values are pased to lsode() function.lsode stands for Livermore Solver for Ordinary Differential Equation which gives us Final result.Final result contains a matrix with two column with 1st column as angular displacment and 2nd column as angular velocity.

1)Function "odefunction" is defined.

function theta_dt = odefunction(theta,t),
         L=1;        %length of pendulum in m.
         m=1;        %mass of bob in Kg
         b=0.05;     %Damping cofficent in Ns/m^2.
         g=9.81;     %acceleration due to gravity in m/s^2.
 
          theta_dt=zeros(2,1);
          theta_dt(1)=theta(2);
          theta_dt(2)=-g/L*sin(theta(1))-b/m*theta_dt(1);
 
 end

2)Main Programme (Code is in OCTAVE).

clear all,close all,clc;

[result]=lsode("odefunction",[0;3],t=linspace(0,20,100));
%Result section 1st column will contain angular displacment and 2nd will contain 
%angular velocity

    
%Creation of plot in windows,two figures are created,1st Figure will contain two subplots 
% 1)Angular Displacment,Angular Velocity Vs time 2)Velocity vs Displacment.
% 2nd Figure contains animation

figure(1),
ct=1;
for i=1:length(result(:,1))
     L=1;
     x0=0;                      
     y0=0;
     
     x1=L*sin(result(i,1));
     y1=-L*cos(result(i,1));
     
     %1st plot() gives a black color line for reference.
     plot([-1 1],[0 0],'color','k','linewidth',3),
     
     
     %2nd plot() function gives blue color line of pendulum length.
     %2nd and 3rd plot() needs to be saved so "hold on" is used
     hold on, 
     plot([x0 x1],[y0 y1],'color','b','linewidth',4),
     
     %3rd plot() functions produces red color bob.
     %Arguments passed are marker,markersize and markerfacecolor.
     % marker is a circle'o' and markersize is 40 and with facecolor red "r".
     plot(x1,y1,'marker','o','markersize',40,'markerfacecolor','r'),
     axis([-1 1 -1.25  0]);       
     hold off,
     pause(0.03);
     
     %Code to save the plot as png.
     %These newly created png files are stiched together to form animation.
     Text=sprintf("pendulum%d.png",ct);
     saveas(1,Text); 
     ct+=1; 
end  

figure(2),
%SubPlot of Angular displacment vs time.
subplot(1,2,1);                                        
plot(t,result);
legend("Angular Displacment","Angular Velocity");
title("Angular Displacment,Angular Velocity Vs time");
xlabel("Time");
%SubPlot of Angular velocity vs time.
subplot(1,2,2);
plot(result(:,1),result(:,2),'color','c');
xlabel("Angular Displacment radians");
ylabel("Angular Velocity radians/sec^2");

print -dpng "PlotPendulum.png"

Explanation:

a)First a Function named 'odefunction' is created.Two arguments are passed in this case theta and t.Function returns a vector of dimension 2 which is given below

`[ d/dt(theta1) ; d/dt(theta2)]=[ theta1 ; -g/L*sin(theta1)-b/m*d/dt*(theta1)`

b) 'lsode' function is used to solve this 2nd order ODE,here 3 arguments are passed, function "odefunction" and intial values and time span.Initial values are `theta(t=0)=0; (d/dt(theta)) at(t=0) =3`.Time span is a 1D array containing time from 0 to 20`

c)Output from "lsode" function is a matrix with two columns,first column contain angular displacment and second column contains angular velocity.

d)Creation of plot in windows,two figures are created,2nd Figure will contain two subplots
1)Angular Displacment,Angular Velocity Vs time 2)Velocity vs Displacment. and subplot(1,2,1) contain Angular Displacment,Angular Velocity Vs time and subplot(1,2,2) contain Velocity vs Displacment

e)1nd Figure contains animation.figure(1) function is used to create 1st figure which has pendulum animation. 'For' loop is used for this purpose.Initial points of pendulum x0,y0 are fixed and x1 and y1 will vary its position.x1 and y1 and `x1=L*sin(theta) and y1=-L.cos(theta)`this value of`theta` will be taken from 1st column of matrix result.

f)At the position x1,y1 we will place a mass,so for this purpose plot() function is used.Arguments passed are marker,markersize,markerfacecolor x1,y1 positions are passed down as argument with a marker "o" and marker size is 40 and markerfacecolor is red.

h) Following Code is used to save the figures

%Code to save the plot as png.
%These newly created png files are stiched together to form animation.
Text=sprintf("pendulum%d.png",ct);
saveas(1,Text);
ct+=1;

Output:

Plot of Angular Displacment,Angular Velocity vs time.

Youtube Video link:

https://youtu.be/OQaY7QWgL4g

Error and its Explanation:

a)To get the proper animation hold on and hold off should be used and they should be positioned carefully.'hold off' command was not used hence it was showing this type of plot,But after positiong of 'hold off' it gave correct results

solution