1

I'm doing a project required to compute the precise(up to a certain precession N) numerical derivatives of a function. The usual approach was to use the finite difference types of algorithms, i.e.

(f(x+h)-f(x))/h 

However, there's not many document to mention how small the h should be, or how to check the convergence or the stability of the results, as h=0 had a different meaning in the numerical computation than symbolic computation.

How to tell how small h is in the finite difference for the precession N? Also, is there any numerical differentiation method that's stable under recurrent application, since the aim was to apply the differentiation up to 1000 times?

1 Answers1

1

In IEEE754 double precision the optimal value of h depends on f(x) and the value of f''(x), the second derivative of f at x. If you have an upper bound M of its absolute value, then the optimal h is (Scilab code below)

h_opt = 2*sqrt(%eps*f(x)/M)

For example, if f(x)=cos(x) and try to approximate f'(1), you have M=1 and

--> h_opt = 2*sqrt(%eps*cos(1))

h_opt =

   2.191D-08

If you try values of h between 10^-1 and 10^-16 and compute the difference between the forward approximation and the true value -sin(1)

format e
h=10.^-(1:16)';
[h -sin(1)-(cos(1+h)-cos(1))./h]

ans =

   1.000D-01   2.559D-02
   1.000D-02   2.687D-03
   1.000D-03   2.700D-04
   1.000D-04   2.701D-05
   1.000D-05   2.702D-06
   1.000D-06   2.701D-07
   1.000D-07   2.806D-08
   1.000D-08  -3.025D-09
   1.000D-09   1.302D-07
   1.000D-10   3.522D-07
   1.000D-11   3.522D-07
   1.000D-12   7.807D-05
   1.000D-13  -1.032D-03
   1.000D-14   2.299D-03
   1.000D-15   4.671D-02
   1.000D-16  -8.415D-01

you can see that the best precision is attained for h=10^-8, which is well predicted by the magnitude of h_opt. This video https://youtu.be/yeFJWtSnV7Y?t=2013 (in french) explains the computation of the optimal h.

Stéphane Mottelet
  • 2,853
  • 1
  • 9
  • 26
  • But how to estimate the lowest upper bound of `f''(x)` if `f('x)` was not known? Because for h smaller than 1e-8 where M to the function was not known and set to a large value, it's clearly not accurate. – ShoutOutAndCalculate Jul 05 '23 at 17:38
  • 1
    There is no alternative, unless you can use a complex step (see https://nhigham.com/2020/10/06/what-is-the-complex-step-approximation/) – Stéphane Mottelet Jul 06 '23 at 18:31
  • Is there any paper with the Cauchy's integral formula? Does the integration better control the rounding error better than the finite differences? – ShoutOutAndCalculate Jul 06 '23 at 18:58
  • 1
    Using complex Cauchy integrals allow to compute derivatives of order > 1, but for first derivative the complex step finite difference gives exact (up to machine precision) results.See https://www.hedonisticlearning.com/posts/complex-step-differentiation.html – Stéphane Mottelet Jul 07 '23 at 13:40