How do I write code for solving partial derivatives numerically?

As stated in the title. I am trying to write a function which evaluates the partial derivative at two points (a,b) for f.

However, the output of the partial derivative evaluated at (0,0) is way too large.

My supposition is that my algorithm for calculating the partial derivative is wrong. But I don't see how.

It has been a long time since I've last used MATLAB, so I do apologise if I've made some errors or used a inefficent way of writing my code.

My code is below:

function derivative = PartialDeriv(f, a, b, i)
h = 0.0001;
fn=zeros(1,2);
if i == 1
fn(i) = (f(a+h,b)-f(a,b)/h);
elseif i==2 fn(i) = (f(a,b+h)-f(a,b)/h);
end
derivative = fn(i);
end

Calling my function I get:

PartialDeriv(f, a, b, i)

where f is

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

I get:

f = -1.9998e+04

Doing it by hand I should get -2.

The i which is seen among the parameters for:

 PartialDeriv(f,a,b,i)

denotes my index, inorder to distinguish the partial derivative with respect to x and y.

Meaning that fn(1) is the partial derivative with respect to x and fn(2) is the partial derivative with respect to y.

1 Answer

Instead of having $f$ accept two arguments, $a, b$, let it accept an argument which is a vector: a_vec = [a1, a2, ..., an]. You seem to need only n=2, but this format will enable you to have the code work for any dimension.

Now:

===========

function derivative = PartialDeriv(f, a_vec, i)

h = 0.0001;

a_dim = length(a_vec)

zero_vector = zeros(1, a_dim)

fn = zero_vector;

for i == 1:a_dim, increment_vec = zero_vector increment_vec(i) = h

fn(i) = ( f( a_vec + increment_vec ) - f( a_vec ) ) / h; end

derivative = fn(i);

end

===========

You might want do use double-sided finite differences instead of the above one-sided one:

fn(i) = ( f( a_vec + increment_vec ) - f( a_vec - increment_vec) ) / (2*h);