Week 3 - Solving second order ODEs

Solving 2^nd order ODE

Aim: To solve the equation regarding the motion of pendulum and generate animation regarding simple pendulum.

Governing Equation:

In this we taken pendulum as an example.

Movement of pendulum depends on the Newtons second law. When this law is written down we get second order ODE and by using this ODE we can compute the position of the ball w.r.t to time.

This ODE represents the equation of motion of a simple pendulum with damping

In the given equation

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

In this problem we are taking below values for solving the ODE

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

We are simulating the motion of ball between 0-20 sec, taking initial values of ODE as angular displacement=0, angular velocity=3 rad/sec.

We have breaked this 2nd order ODE into 2 single order differential equation and solving the ODE.

Coding Process:

Code1: This code is regarding the function of ODE which we are going to solve

function [dtheta_dt] = pendulum_ball(t,theta,b,m,g,l)
  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

In this we have converted 2nd order ODE into single ODE's.

As we assumed theta as array of theta1&2 .We are assuming theta1 as theta of 1 and theta2 as theta of 2.

We have solved and made 2nd orde ODE into 2 single ODE's like

dtheta1_dt = theta2

dtheta2_dt = -(b/m)*theta2 -(g/l)*sin(theta1)

and assuming the dtheta1_dt and dtheta2_dt as dtheta_dt.

Afer that we will get

dtheta_dt = [dtheta1_dt; dtheta2_dt].

This programme is used to solve the ODE in the main pendulum function.

Code2: This code is used to solve the ODE and create an animation of the pendulum

% Simple pendulum Differential Equation Solver
  clear all
  close all
  clc

% Input Values
  b = 0.05   % damping coefficient
  g = 9.81   % gravity in m/s2
  l = 1      % lenght of pendulum in m 
  m = 1      % mass of ball in Kg
  
% Initial values 
  initial_a = [0; 3]   
% Time range
  t_span = linspace (0,20,250)
  
% Solving Differential Equation:

  [t, result_data] = ode45(@(t,theta) pendulum_ball(t,theta,b,m,g,l),t_span , initial_a)

% After getting our values we can simulate the pendulum
  ct = 1
  for i = 1:length(result_data(:,1))
  % Defining the coordinate points of pendulum
  x0 = 0
  y0 = 0
  y1 = -l*cos(result_data(i,1))
  x1 = l*sin(result_data(i,1))
  % Plot and fill commands to create a fixed support in animation  
  plot([ -1.5 1.5 ],[ 0 0 ],'color','k','linewidth',3)
  hold on  
  plot([-1.5 -1.5 1.5 1.5 ],[ 0 0.15 0.15 0 ],'color','k','linewidth',3) 
  fill([ -1.4 -1.4 1.4 1.4 ], [ 0 0.15 0.15 0 ],'k')
  % title command to name the plot
  title('Pendulum Animation','fontsize',25,'verticalalignment','middle')
  % Plot command to generate  centre axis in plot
  plot([ 0 0 ],[ 0.3 -1.25 ],'linestyle','-.','color','b','linewidth',2)
  % To display the pendulum ball and string
  plot(x0,y0,'o','markersize',15)
  plot([x0 x1],[y0 y1],'color','g','linewidth',3)
  plot(x1,y1,'.','markersize',50,'color','r')
  % Defining the axis limits
  axis([-2 2 -2 1.5 ])
  % naming the labels of plot
  xlabel('X-axis')
  ylabel('Y-axis')
  hold off
  Imagea = sprintf('Ball%d.png',ct)
  saveas(gca,Imagea)
  ct = ct + 1
end

First we have defined the inputs for the proble like damping coefficient,length of pendulum,mass,gravity.

After that we defined the initial values of the differential equation like displacement as 0 and velocity as 3m/s2

Getting an array of 250 time values between 0 to 20 seconds by using linspace command.

After that we have used ode45 function to get the results

first we have to give the result array name like in this programme we used 'result_data' and t represents the time values and used ode45 function to solve.

In ode function we have first defined main function variables time ,theta after that we called main ode function 'pendulum_ball' here and after that we give time span and initial values.

After getting the results we will be getting displacement in first column of result_data vector.

After that we have opened 'for loop' and defined the range is from 1 to length of result_data column 1

After that I defined the coordinate points of pendulum like x0,y0,x1,y1.

x1 and y1 are calculated by general trignometric relation by assuming the length of pendulum as hypotenuse.

I have used negative to correct the position of pendulum as I got error like plot coming on wrong axis So I used negative sign.

After that I used plot and fill command to create a fixed support for pendulum and title command to name the plot including hold command to hold all plots.

Next also I used plot command to create an centre axis by changing line style

After that I used marker command to display the end of pendulum as ball and increased size and color.

After that I increased the line width of plot between x0,y0 and x1,y1

Defined axis limits and labels x&y axis.

Used hold off to stop plots saving for next iteration

Used sprintf to save the images in array with different names by using the varible ct

lastly we used saveas command to save the plots with each iteration in the programme directory.

Code3: To generate plots for the same ODE.

% Simple pendulum Differential Equation Solver
  clear all
  close all
  clc

% Input Values
  b = 0.05   % damping coefficient
  g = 9.81   % acceleration
  l = 1      % lenght of pendulum
  m = 1      % mass of pendulum
  
% Initial values
  
  initial_a = [0; 3]
  
% Time range

  t_span = linspace (0,20,500)
  
% Solving Differential Equation:

  [t, result_data] = ode45(@(t,theta) pendulum_ball(t,theta,b,m,g,l),t_span , initial_a)
  
% plotting the values
  % Plotting between time & Velocity,Displacement
  figure(1)
  hold on
  grid on
  box on
  xlabel('Time')
  ylabel('Velocity & Displacement')
  plot(t,result_data)
  % Plotting between Velocity & Displacement
  figure(2) 
  hold on
  grid on
  box on
  xlabel('Displacement')
  ylabel('Velocity')
  plot(result_data(:,1),result_data(:,2))
   

In this code we used the same code we used for pendulum animation for getting results.

After that we used plot command to plot between results(2 columns) and time array.

After that we used plot command to plot between results column 1 and 2.

We used grid on and box on command to switch it on.

Errors and Corrections:

Firstly I am unable to plot the the pendulum but after seeing the video I got some information and after that I plotted the pendulum.

After creating the code , I am not getting plot on bottom but after the correction in support session I resolved the issue by changing the x1 and y1 formulas and using the negative sign.

After that I have got problem regarding save as function and I corrected by keeping the function inside for loop but I kept outside for loop before.

Results:

Youtube Link : https://youtu.be/QlE-GEtGvkk