MATLAB Examples


% Author:  Darren Rowland
% Version: 1.0
% Date:    8 November 2010

Points with equal spacing

When the x-values are equally spaced, the derivatives at the interior points as calculated by DIFFXY are approximated in the same way as the native Matlab function GRADIENT, i.e. using the central difference method. The endpoints are different because DIFFXY uses a higher-order approximation.

x = linspace(-1,2,20);
y = exp(x);
dy = diffxy(x,y);
dy2 = diffxy(x,dy);  % Or, could use >> dy2 = diffxy(x,y,[],2);

Dy = gradient(y)./gradient(x);
Dy2 = gradient(Dy)./gradient(x);
hold on
title('Relative error in derivative approximation')
legend('diffxy: dy/dx','diffxy: d^2y/dx^2',...
       'gradient: dy/dx','gradient: d^2y/dx^2')
hold off

Points with unequal spacing

DIFFXY was developed specifically to cope with unevenly spaced data. Consider the result when GRADIENT attempts to approximate the derivative of y=exp(x) when the x-values are randomly distributed.

x = 3*sort(rand(20,1))-1;
y = exp(x);

title('Approximation of d(exp(x))/dx by GRADIENT')

Compare this to the approximation obtained by DIFFXY

title('Approximation of d(exp(x))/dx by DIFFXY')

The distinction is even more pronounced as higher-order derivatives are approximated.

dy2 = diffxy(x,y,[],2);
Dy = gradient(y)./gradient(x);
Dy2 = gradient(Dy)./gradient(x);
title('Approximation of d^2(exp(x))/dx^2 for unequally spaced x')
legend('analytical','diffxy: d^2y/dx^2','gradient: d^2y/dx^2',...

N-dimensional data

As well as operating on univariate data as above, DIFFXY can be used in general multidimensional space, obtaining the derivative in any single dimension one-at-a-time. This is different to GRADIENT, which returns the approximate derivatives in each coordinate direction simultaneously.

x = sort(rand(10,1));
y = sort(rand(15,1));
[X,Y] = ndgrid(x,y);

f = sqrt(X).*Y.*cos(2*Y);
% Exact partial derivatives
f_x = f./(2*X);
f_y = sqrt(X).*(cos(2*Y)-2*Y.*sin(2*Y));

dfdx = diffxy(x,f,1);
dfdy = diffxy(y,f,2);

hold on
hold off
view(gca,[-68.5 40]);
hold on
hold off
view(gca,[-68.5 40]);