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:

load Dataset.mat  % Loading the data for two players

          

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: https://sameradeeb-new.srv.ualberta.ca/introduction-to-numerical-analysis/numerical-differentiation/