Numerical differentiation for function including vector

I have implemented a numerical differentiation with central differences. I differentiate the function with respect to a. This appears to be relatively simple for functions containing scalars, e.g., . However, once the function also contains a vector, e.g., , where c is a 3 by 1 vector, I am struggling to implement this.

Below I share my code:

rng default
b = randn(1);
c = randn(3, 1);

f1 = @(a) a^4 * b^3;
%f2 = @(a) (a^4 * b^3 * c.').';

a = 1; 
k = 3; % no. derivative
h = 0.1; % step size
cent_diff = zeros(k+1,1);

for i = 0:k
    cent_diff(i+1, 1) = sum((-1)^i*nchoosek(k,i)*f1(a+(k/2-i)*h));
end

cent_diff = sum(cent_diff)/h^k;
cent_diff = 3.7304

% check result with symbolic differentiation
syms a
f1 = a^4 * b^3;
diff_3 = diff(f1, a, 3);
a = 1;
vpa(subs(d3)) = 3.7304

The differentiation of f2 can be easily achieved using the symbolic approach. How can I get the corresponding result using finite differences?

[Walter Roberson]

    cent_diff(i+1, 1) = sum((-1)^i*nchoosek(k,i)*f1(a+(k/2-i)*h));

What is the sum over? k i h are scalar. f1 returns a scalar for scalar input.

If you were looking ahead to the case where f1 might return a vector, then you would not want to sum over that vector.

[M R]

I refer to this wiki article under point 3 higher order derivatives with the finite difference method. For f1 it is fine because, as you mentioned, it returns a scalar for a scalar input. However, for f2 it is different because c is a 3 by 1 vector. I am able to replicate the result of the symbolic approach by using the explicit formula for the 3rd derivative of the finite difference method. However, the difficulties arise when using the generalized formula for the k-th derivative as described in the wiki article.

The code to see that it works with the explicit formula is as follows:

f2 = @(a) (a^4 * b^3 * c.').';
D = (f2(a+(2*h)) - 2*f2(a+h) + 2*f2(a-h) - f2(a-(2*h))) / (2*h^3)
D =
    6.8411
   -8.4263
    3.2162

Using the symbolic approach:

f2 = (a^4 * b^3 * c.').';
diff_3 = diff(f2, a, 3);
f2_diff = vpa(subs(d3))
f2_diff =
 6.8411
-8.4263
 3.2162

Thus, the question now is how to use the generalized version described above when the inputs and outputs are not scalars.

Answer 1

[Walter Roberson]

The wikipedia article shows .

for i = 0:k
    cent_diff(i+1, 1) = sum((-1)^i*nchoosek(k,i)*f1(a+(k/2-i)*h));
end

inside that loop, i is scalar. There is no sum at that level. You are effectively recording the terms there

cent_diff = sum(cent_diff)/h^k;

and that is where you are doing the summation.

Now if you use a function such as f2 that returns a vector, then your stray sum() is adding together the individual terms for all of the vector elements. Which is not what you want.

[Walter Roberson]

rng default
b = randn(1);
Nc = 3;
c = randn(Nc, 1);

f1 = @(a) a^4 * b^3;
f2 = @(a) (a^4 * b^3 * c.').';

a = 1; 
k = 3; % no. derivative
h = 0.1; % step size
cent_diff = zeros(Nc, k+1);

for i = 0:k
    cent_diff(:, i+1) = (-1).^i .* nchoosek(k,i) .* f2(a+(k/2-i)*h);
end

cent_diff

cent_diff = 3×4

0.4985 -1.0394 0.6965 -0.1488 -0.6141 1.2803 -0.8579 0.1833 0.2344 -0.4887 0.3275 -0.0700

cent_diff = sum(cent_diff,2)./h^k;

cent_diff

cent_diff = 3×1

6.8411 -8.4263 3.2162

Answer 2

[Sulaymon Eshkabilov]

These answer discussions provide appropriate points how to apply the finite difference method:

https://www.mathworks.com/matlabcentral/answers/165748-how-do-you-solve-a-nonlinear-ode-with-matlab-using-the-finite-difference-approach?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/answers/366599-problem-with-finite-difference

https://www.mathworks.com/matlabcentral/answers/328722-matlab-simple-loop-for-different-function-variables-finite-difference?s_tid=srchtitle

[M R]

Thanks for the links. Unfortunately, they are not really helpful. I know how to implement the 1st and 2nd order derivatives using the finite difference method. However, in this question I am interested in the general solution for the nth order difference.