x values when taking a numerical derivative

I want to calculate the first derivative f(x) = dy/dx for data that is irregularly spaced in x. I think (correct me if I'm wrong) this can be done by diff(y)./diff(x)

But the resulting vector is one less than the length of the x and y vectors, so what are the corresponding x values? I want to then calculate x dy/dx, so it would be helpful to know what to use.

Answer 1

[Star Strider]

‘I think (correct me if I'm wrong) this can be done by diff(y)./diff(x)’

Not actually wrong, simply mistaken with respect to expecting that the result of diff will be what you want.

Use the gradient function (or Torsten’s approach, that does essentially the same operation) instead:

dydx = gradient(y) ./ gradient(x);

.

[L'O.G.]

Thank you!

[L'O.G.]

Just to make sure I understand, would you be able to tell me why you recommend using gradient rather than diff here? They give quite different results.

[Star Strider]

As always, my pleasure!

The gradient funciton returns an estimate of the derivative with the length (and more generally, size) of the output being the same as the input. (It can also be used to calculate spatial derivatives.) It is more accurate than diff, that simply calculates sequential differences. It computes a numerical derivative by approximating it with central differences, as described in the documentation section on Algorithms, and also described by Torstein’s code, that describes it in more detail.

.

[L'O.G.]

Thank you again. Sorry, just one more question (I promise this will be it for this thread): the support page for gradient says that FX = gradient(F) corresponds to /, so why do you recommend dividing by the gradient in x in this case?

[Torsten]

the support page for gradient says that FX = gradient(F) corresponds to /, so why do you recommend dividing by the gradient in x in this case?

Because in the case FX = gradient(F), the distance between the x values is assumed to be 1 for all of them.

Or do you think MATLAB knows how your x vector looks like if you don't supply it ?

[Star Strider]

@Torsten — Thank you!

Answer 2

[Sulaymon Eshkabilov]

Yes, you are right, e.g.:

x = [0   .2 .3 .45 .65 .75 .96]

x = 1×7

0 0.2000 0.3000 0.4500 0.6500 0.7500 0.9600

y = [-3  10 11  12  13 12  9]

y = 1×7

-3 10 11 12 13 12 9

dydx = diff(y)./diff(x)

dydx = 1×6

65.0000 10.0000 6.6667 5.0000 -10.0000 -14.2857

yyaxis left
plot(x, y, 'b-o', 'MarkerFaceColor', 'y', 'DisplayName', 'y(x)')
ylabel('y(x)')
yyaxis right
plot(x(1:end-1), dydx, 'r--p', 'MarkerFaceColor','c','DisplayName', 'dy/dx')
ylabel('dy/dx')
xlabel('x')
grid on
legend show

[L'O.G.]

Thanks. So just to be clear, why does dydx(end-1) correspond to x(end)?

[Sulaymon Eshkabilov]

difference1 is between 1 and 2 and then 2 and 3, etc., and thus end-1

Answer 3

[Tushar Behera]

Hi L'O.G,

I believe you want to calculate derivate for two separate datasets. In order to do that you can use 'diff(y)./diff(x)' for example:

clc;
clear;
close all;
x=linspace(0,2*pi,100);
y=sin(x);
yp=cos(x);
dx=diff(x);
dy=diff(y);
yp_hat=dy./dx;
err=yp(1:end-1)-yp_hat;
figure;
subplot(1,2,1);
plot(x,y);
hold on;
plot(x(1:end-1),yp_hat)
xlabel('x');
ylabel('y');
legend('original function','Approx derivative');
grid on;

subplot(1,2,2);
plot(x(1:end-1),err);
xlabel('x');
ylabel('error');

here 'x' and 'y' are two vectors and by using 'diff(y)./diff(x)' you can calculate the first order derivate which is 'cos(x)'. To answer the question which values of x corresponds to which 'yp_hat' . you can get that by using,

x(1:end-1)

I hope this solves your query.

Regards,

Tushar

Answer 4

[Torsten]

Use

n = numel(y);
dydx(1)     = (y(2) - y(1))/(x(2) - x(1));
dydx(2:n-1) = (y(3:n) - y(1:n-2))./(x(3:n) - x(1:n-2));
dydx(n)     = (y(n) - y(n-1))/(x(n) - x(n-1));   

Answer 5

[Sulaymon Eshkabilov]

There will be some significantly different results from diff() and gradient() if the increment of x varies. See this simulation:

x = [0   .2 .3 .45 .65 .75 .96];
y = [-3  10 11  12  13 12  9];
dy1 = diff(y)./diff(x);
dy2 = gradient(y)./gradient(x);
for ii = 1:length(x)-1
dy3(ii) = (y(ii+1)-y(ii))/(x(ii+1)-x(ii));
end
N = numel(y);
dy4(1) = (y(2)-y(1))/(x(2)-x(1));
dy4(2:N-1) = (y(3:N)-y(2:N-1))./(x(3:end)-x(2:end-1));
dy4(N) = (y(end)-y(end-1))/(x(end)-x(end-1));
plot(x(1:end-1), dy1, 'b-o', 'MarkerFaceColor', 'y', 'DisplayName', 'diff', 'markersize', 13)
hold on
plot(x(1:end), dy2, 'rs--', 'MarkerFaceColor', 'c', 'DisplayName', 'gradient')
plot(x(1:end-1), dy3, 'g--p', 'MarkerFaceColor','y','DisplayName', 'Loop computed difference', 'MarkerSize', 10)
hold on
plot(x(1:end), dy4, 'k--h', 'MarkerFaceColor','c','DisplayName', 'vectorized: gradient')
ylabel('dy/dx')
xlabel('x')
grid on
legend show

% Note that as the increment of x gets smaller the error (offset) will also
% diminish. See this example: dx = 0.063467 vs. dx = 0.006289:

x=linspace(0,2*pi,100);
dx = x(2);
y=sin(x);
ANS=cos(x);
dY1=diff(y)./diff(x);
dY2 = gradient(y)./gradient(x);
for ii = 1:length(x)-1
dY3(ii) = (y(ii+1)-y(ii))/(x(ii+1)-x(ii));  % The same as diff()
end
n = numel(y);
dY4(1) = (y(2)-y(1))/(x(2)-x(1));
dY4(2:n-1) = (y(3:end)-y(1:end-2))./(x(3:end)-x(1:end-2));
dY4(n) = (y(end)-y(end-1))/(x(end)-x(end-1));

E1=ANS(1:end-1)-dY1;
E2 = ANS-dY2;
E3 = ANS(1:end-1)-dY3;
E4 = ANS-dY4;
fprintf(['Norm of errors @ dx = %f: ' ...
    'E_diff = %f; E_gradient = %f; ' ...
    'E_loop = %f; E_grad_vect = %f \n'], [dx, norm(E1) norm(E2) norm(E3) norm(E4)])

Norm of errors @ dx = 0.063467: E_diff = 0.223238; E_gradient = 0.004770; E_loop = 0.223238; E_grad_vect = 0.004770 

% 10 times smaller incremental step of x than the previous example leads to
% the reduction of error norm to more than 3 times
x=linspace(0,2*pi,1000);
dx = x(2);
y=sin(x);
ANS=cos(x);
dY1=diff(y)./diff(x);
dY2 = gradient(y)./gradient(x);
for ii = 1:length(x)-1
dY3(ii) = (y(ii+1)-y(ii))/(x(ii+1)-x(ii));  % The same as diff()
end
n = numel(y);
dY4(1) = (y(2)-y(1))/(x(2)-x(1));
dY4(2:n-1) = (y(3:end)-y(1:end-2))./(x(3:end)-x(1:end-2));
dY4(n) = (y(end)-y(end-1))/(x(end)-x(end-1));

E1=ANS(1:end-1)-dY1;
E2 = ANS-dY2;
E3 = ANS(1:end-1)-dY3;
E4 = ANS-dY4;
fprintf(['Norm of errors @ dx = %f: ' ...
    'E_diff = %f; E_gradient = %f; ' ...
    'E_loop = %f; E_grad_vect = %f \n'], [dx, norm(E1) norm(E2) norm(E3) norm(E4)])

Norm of errors @ dx = 0.006289: E_diff = 0.070283; E_gradient = 0.000147; E_loop = 0.070283; E_grad_vect = 0.000147