How to calculate a derivative of function numerically?

I have a continuous function of two variables:

f = @(x,y) x.^2 + (y-1).^2

How do I calculate a derivative of this function numerically? I tried two different ways and got errors:

% First try
df = @(x,y) diff(f(x,y),x)
df(1,2)
ans = []

% Second try
df = @(x,y) diff(f,x)
df(1,2)
Undefined function 'diff' for input arguments of type 'function_handle'

I must not set the derivative separately.

[Denis Perotto]

Like this:

[Denis Perotto]

I found out that this solves my problem:

syms x y

X = [2 3 4];
Y = [0 1 4];
f = @(x,y) x.^2 + (y-1).^2;

g1 = eval(['@(x,y)' char(diff(f,x))]);
g1(1,2)

But it never works if function 'f' contains interpolation (for example, 'interp1'). So sad =(

[Bruno Luong]

What sad about? If you want numerical differential then program it as you like, for the starter:

f = @(x,y) x.^2 + (y-1).^2;
dxdy=1e-6;
df = @(x,y) [f(x+dxdy,y)-f(x,y); f(x,y+dxdy)-f(x,y)]/dxdy;

df(5,5)

[Denis Perotto]

Well, that's okay, but I hoped there is an inner function to evaluate the procedure without extra code

[onur karakurt]

Thanks so much,

But I have another question that showing in picture:

[madhan ravi]

ask a separate question @onur

[John D'Errico]

You can use higher order finite differences.

x = 2;
fun = @(x) sin(x);

-sin(x)
ans =
        -0.909297426825682

dx = 1e-6;
(fun(x - dx) - 2*fun(x) + fun(x + dx))/dx^2
ans =
        -0.909494701772928

dx = 1e-7;
(fun(x - dx) - 2*fun(x) + fun(x + dx))/dx^2
ans =
        -0.899280649946377

But be careful, as if you use too small a delta, the scheme will go to hell quickly. Even derivest has difficulties for higher order derivatives, as things get unstable.

[df,errest] = derivest(fun,2,'deriv',2)
df =
        -0.909297426825812
errest =
      2.98009666234719e-12

Answer 1

[Jan]

A numerical differentiation means a variation of the argument and a quotient of the differences:

df = (f(x + h) - f(x)) / h           % one sided
df = (f(x + h) - f(x - h)) / (2*h)   % two sided

Answer 2

[John D'Errico]

You can use symbolic tools to compute the derivative, if you have that toolbox. But if you wish todo so numerically, then you need to use numerical tools. They will never give a perfectly correct answer of course, since they at best use approximations to the result.

You can also use my derivest tools on the File Exchange. They use numerical approximations to compute the derivative, as well as estimate the error in that estimate.

For example, consider the function sin(x). We know the derivative of sin(x) is cos(x).

fun = @(x) sin(x);

We know that the derivative is defined in terms of a limit. Thus,

x = 1;
dx = 1e-12;
df = (fun(x + dx) - fun(x))/dx
df =
         0.540345546085064

At x==1, we know the derivative is cos(1).

cos(1)
ans =
          0.54030230586814

So the simple estimate, using dx == 1e-12 was not too terribly poor. Be careful though. Make dx too small, and things go to hell.

dx = eps
dx =
      2.22044604925031e-16

df = (fun(x + dx) - fun(x))/dx
df =
                       0.5

which is clearly incorrect, even though MATLAB can compute the limit using symbolic tools.

syms x dx
limit((sin(x + dx) - sin(x))/dx,dx,0,'right')
ans =
cos(x)

Derivest is able to get itright, since it is designed to do just that.

[df,errest] = derivest(fun,1)
df =
          0.54030230586814
errest =
      1.47246694029184e-15

derivest nails it, and it tells you that it is confident of the result to within roughly 1.5e-15.

A better approximation comes from a higher order estimate.

dx = 1e-12;
df = (fun(x + dx) - fun(x - dx))/(2*dx)
df =
         0.540290034933832

Again, in the limit, as x gets too small, things will get bad here too.

dx = eps
dx =
      2.22044604925031e-16

df = (fun(x + dx) - fun(x - dx))/(2*dx)
df =
                       0.5

Limits don't always work that well when you push them in term of numerical computations. Derivest succeeds because it uses an extrapolatory scheme.