Week 3 - Solving second order ODEs

Solving Second Order ODE in Matlab/Octave

An ordinary differential equation (ODE) is a differential equation which has one or more functions of one independant variable and the derivatives of those functions. ODEs are widely used in multiple fields, where description of change uses derivatives and differentials. Differential equation is a result that describes dynamically changing phenomena, where quantities are defined as the rate of change of other quantities.

This program is written to solve a second order ODE in Octave. The second order ODE considered is the motion of a simple pendulum, as shown above.

To solve the equation, the first step is to convert the second order ODE into 2 first order ODEs. This is done by introducing new variables θ1 and θ2, as shown below:

Here, new variables are introduced, such that the second degree ODE is converted to two First degree ODE. The equations are written in matrix form, which is accepted by the software.

1) Program to solve the ODEs to animate the damped motion of a pendulum:

https://drive.google.com/file/d/1UX_IZqW0Uhf-dL4QfDMipaeHruoCMr8W/view?usp=sharing

%program for solving a first order differential equation 

clear all
close all
clc

b = 0.05;          %damping coefficient
g= 9.81;           %acceleration due to gravity
l = 1;             %pendulum length
m = 1;             %mass of bob

theta_0 = [0;3];                    %specifying the position and angular velocity
t_span = linspace(0,20,1000);       %time span divided into intervals

%ode45:inbuilt differential equation solver to solve for position and velocity using the inputs given
[t,Results] = ode45(@(t,theta) pendulum_ode_damped(t,theta,b,g,l,m),t_span,theta_0);

subplot(5,4,[13,14,17,18]);                    %subplot command to create multiple plots within one figure
plot(t,Results(:,1),'color','r')               %plotting the graph for displacement
hold on;
plot(t,Results(:,2),'color','g')               %plotting the graph for velocity
xlabel('Time');
ylabel('Position,Velocity');
title('Velocity/Displacement vs Time');                                 %title of the graph
h = legend('Position','Velocity');
legend(h,'location','southoutside','orientation','horizontal');         %legend to distinguish the data sets
hold off;

subplot(5,4,[15,16,19,20]);                    %subplot command to create velocity versus position graph
plot(Results(:,1),Results(:,2))                %plotting the obtained values of postions and velocities
xlabel('Position');                        
ylabel('Angular velocity');
title('Angular Velocity vs Position');

ct = 1;                                         
for i=1:length(Results(:,1))                   %for loop for creating animation
  x0=0;                                        %co-ordinates for fixed point of pendulum
  y0=0;
  x1=l*sin(Results(i,1));                      %co-ordinates for centre of the bob of the pendulum
  y1=-l*cos(Results(i,1));
  subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);              %graph for pendulum animation
  plot([-1 1],[0 0],'linewidth',4,'color','k');           %line to represent the ceiling for suspending the pendulum
  hold on;
  line([x0 x1],[y0 y1],'linewidth',2,'color','b');        %line to represent the string of the pendulum
  hold on;
  plot(x1,y1,'-o','markersize',20,'markerfacecolor','r')  %marker to represent the bob of the pendulum
  grid on;
  axis([-1.5 1.5 -1.5 0.5])
  pause(0.01);
  hold off;
  title('Pendulum Animation');
  
  file_text=sprintf('Pendulum%04d.png',ct)           %Name of the file to store the images
  saveas(gca, file_text)                             %Save command for file_text
  ct=ct+1;
end

First the input parameters, namely damping coefficient(b), acceleration due to gravity(g), length of pendulum(l) and mass(m) are specified. Then the initial position and angular velocities are specified in vector form, assigning theta_0. Then the timespan for the whole process is entered (0-20 seconds), with the span divided into 1000 intervals. The ode45 is an inbuilt command in both matlab and octave, used for solving first order differential equations. It utilises the Runge Kutta Method to solve the equations. The function is used by first specifying the variables (t,theta), then calling the function and its inputs [pendulum_ode_damped(t,theta,b,g,l,m)], then specifying the limits for the integration(t_span) and finally the initial conditions(theta_0). The results are plotted in two columns(t,results) for displacement and velocities respectively.

Next the graphs are plotted with the outputs obtained. The subplot command is useful for plotting multiple graphs in one figure window. The subplot works by specifying the dimensions of the entire figure window in the form of a virtual matrix and then specifying the elements of the matrix where the specific graph is to be entered. In this case the matrix is 5X4. The first plot plots the outputs of ode45 (position and velocity) with respect to time. The legend is also prepared for this plot, differentiating between the two graphs, by specifying the values and position of the legend. Then the graph for angular velocity versus position graph was plotted by specifying the position within the subplot command. The last graph plotted after this was the animation of the pendulum. This was done by executing a for loop with a counter. The co ordinates for the fixed and moving ends of the pendulum was written to simulate the string within the loop, along with the line representing the ceiling. The axes of the graph are also specified along with the bob of the pendulum, which is represented by a marker. All the graphs are given a title. Then the saveas and sprintf commands are used to save the animation in the form of a series of images, which are joined together to create a gif file in the imagemagick software, in the same way done within the robotic arm challenge. The gif is converted to mp4 format and uploaded to youtube.

Youtube link: https://youtu.be/TfN37DjcIe0

2) Function called within the program

https://drive.google.com/file/d/1lldxI9eLxAGM_5bkb1BcBMfGpJhKV1Nz/view?usp=sharing

function [dtheta_dt] = pendulum_ode_damped(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];
  
endfunction

The function is used to take the inputs from the original program to calculate the matrix shown in the first image (Equation 4). The value 'dtheta_dt' is considered by ode45 to calculate the required parameters.

Error: An error occurred during the execution of the program using lsode command. This was solved by using ode45 command.

Conclusion: The ordinary differential equation was solved within the program using the ode45 command. New features like legend, marker, subplot were also learnt.