## Welcome to the World of Modelling and Simulation

### What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ... ### An Example of a Finite Difference Method in MATLAB to Find the Derivatives

In this tutorial, I am going to apply the finite difference approach to solve an interesting problem using MATLAB. This example is based on the position data of two squash players - Ramy Ashour and Cameron Pilley - which was held in the North American Open in February 2013. The positions (in meters) of the left and right feet of the players (Ashour and Pilley, respectively) during an 8 second segment of their match are provided in this link. With the available position data for the two players, we are interested to see the velocity and acceleration of each player.

Since velocity and acceleration are the first and second derivative of the position respectively, using the finite difference approach, we can determine the derivatives numerically. Here, we will use centered finite difference approach for both derivatives, which has an accuracy of second order. A program is written in MATLAB, which evaluates the derivatives numerically using the centered finite difference.

The data shows us the x and y positions of the feet for the two players at different time steps. As we know that the rate of change of position gives us the sense of velocity and rate of change of velocity provides the notion of acceleration accordingly. So, to numerically evaluate the first and second derivatives of the position data, we can easily apply the fundamentals of the centered finite difference using the direct formulae for velocity and acceleration. The following MATLAB codes do the same thing by implementing several 'for loops' to organize the data and calculate the derivatives as well.

## MATLAB Program for Centered Finite Difference:

n = length(Dataset(:,1)); % For for-loop limit
Time = Dataset(:,1);   % First column of dataset is time

% X and Y average of position for Ashour
PL1_Xavg = (Dataset(:,2) + Dataset(:,4))./2;
PL1_Yavg = (Dataset(:,3) + Dataset(:,5))./2;
% X and Y average of position for Pilley
PL2_Xavg = (Dataset(:,6) + Dataset(:,8))./2;
PL2_Yavg = (Dataset(:,7) + Dataset(:,9))./2;

% Resultant Position for Ashour and Pilley
PL1 = sqrt(PL1_Xavg.^2 + PL1_Yavg.^2);
PL2 = sqrt(PL2_Xavg.^2 + PL2_Yavg.^2);

% Calculating velocities for Ashour
% Centered finite difference approach
i=2; % Initiating the loop for velocity
for i = 2:n-1
V_1X(i-1) = (PL1_Xavg(i+1)-PL1_Xavg(i-1))/(Time(i+1)-Time(i-1));
V_1Y(i-1) = (PL1_Yavg(i+1)-PL1_Yavg(i-1))/(Time(i+1)-Time(i-1));
i = i+1;
end
% Calculating velocities for Pilley
for i = 2:n-1
V_2X(i-1) = (PL2_Xavg(i+1)-PL2_Xavg(i-1))/(Time(i+1)-Time(i-1));
V_2Y(i-1) = (PL2_Yavg(i+1)-PL2_Yavg(i-1))/(Time(i+1)-Time(i-1));
i = i+1;
end
% Calculating accelerations for Ashour
% Centered finite difference approach
i=2; % Initiating the loop for acceleration
for i=2:n-2
A_1X(i) = (PL1_Xavg(i+1)-(2*PL1_Xavg(i))+PL1_Xavg(i-1))/((Time(i+1)-Time(i)).^2);
A_1Y(i) = (PL1_Yavg(i+1)-(2*PL1_Yavg(i))+PL1_Yavg(i-1))/((Time(i+1)-Time(i)).^2);
i = i+1;
end
% Calculating accelerations for Pilley
for i=2:n-2
A_2X(i) = (PL2_Xavg(i+1)-(2*PL2_Xavg(i))+PL2_Xavg(i-1))/((Time(i+1)-Time(i)).^2);
A_2Y(i) = (PL2_Yavg(i+1)-(2*PL2_Yavg(i))+PL2_Yavg(i-1))/((Time(i+1)-Time(i)).^2);
i = i+1;
end
% Plotting X position, velocity & acceleration of Ashour
figure(1)
subplot(311), plot(Time,PL1_Xavg,'b'), grid
title('X Position of Ashour with Time')
xlabel('Time'), ylabel('Position')
subplot(312), plot(Time(2:n-1),V_1X(1:n-2),'r'), grid
title('X Velocity of Ashour with Time')
xlabel('Time'), ylabel('Velocity')
subplot(313), plot(Time(2:n-1),A_1X(1:n-2),'g'), grid
title('X Acceleration of Ashour with Time')
xlabel('Time'), ylabel('Acceleration')

% Plotting Y position, velocity & acceleration of Ashour
figure(2)
subplot(311), plot(Time,PL1_Xavg,'b'), grid
title('Y Position of Ashour with Time')
xlabel('Time'), ylabel('Position')
subplot(312), plot(Time(2:n-1),V_1Y(1:n-2),'r'), grid
title('Y Velocity of Ashour with Time')
xlabel('Time'), ylabel('Velocity')
subplot(313), plot(Time(2:n-1),A_1Y(1:n-2),'g'), grid
title('Y Acceleration of Ashour with Time')
xlabel('Time'), ylabel('Acceleration')

% Plotting X position, velocity & acceleration of Pilley
figure(3)
subplot(311), plot(Time,PL2_Xavg,'b'), grid
title('X Position of Pilley with Time')
xlabel('Time'), ylabel('Position')
subplot(312), plot(Time(2:n-1),V_2X(1:n-2),'r'), grid
title('X Velocity of Pilley with Time')
xlabel('Time'), ylabel('Velocity')
subplot(313), plot(Time(2:n-1),A_2X(1:n-2),'g'), grid
title('X Acceleration of Pilley with Time')
xlabel('Time'), ylabel('Acceleration')

% Plotting Y position, velocity & acceleration of Pilley
figure(4)
subplot(311), plot(Time,PL2_Yavg,'b'), grid
title('Y Position of Pilley with Time')
xlabel('Time'), ylabel('Position')
subplot(312), plot(Time(2:n-1),V_2Y(1:n-2),'r'), grid
title('Y Velocity of Pilley with Time')
xlabel('Time'), ylabel('Velocity')
subplot(313), plot(Time(2:n-1),A_2Y(1:n-2),'g'), grid
title('Y Acceleration of Pilley with Time')
xlabel('Time'), ylabel('Acceleration')

The following plot shows the horizontal position, velocity, and acceleration components for Ashour with time.

The following plot shows the vertical position, velocity, and acceleration components for Ashour with time.

The following plot shows the horizontal position, velocity, and acceleration components for Pilley with time.

The following plot shows the vertical position, velocity, and acceleration components for Pilley with time.

Reference: