Approximate derivatives and automatic differentiation

With Automatic differentiation, you get results very similar to what you would get by symbolically manipulating the values, then plugging in.

For example, we have for the function $f(x) = e^x \sin(x)$, a derivative computed by automatic differentiation:

using Roots # calls in ForwardDiff package
f(x) = exp(x) * sin(x)
D(f, 3)(1)     # third derivative of f at 1

-1.6373226945259143

Whereas, Symbolically we have:

using SymPy
x = symbols("x")
f_xxx = diff(f(x), x, 3)

$$2 \left(- \sin{\left (x \right )} + \cos{\left (x \right )}\right) e^{x}$$

And evaluating at 1 gives:

f_xxx(1) |> N

-1.6373226945259145

However numeric derivatives are a different matter. There are formulas along the line of this formula to take a second derivative:

$$~ f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} ~$$

But they get complicated. It gets tempting to instead compute the answer along the lines of

$$~ \begin{align} f'(x) &\approx \frac{f(x+h) - f(x)}{h} \\ f'(x+h) &\approx \frac{f(x+2h) - f(x+h)}{h}\\ f''(x) &\approx \frac{f'(x+h) - f'(x)}{h} \end{align} ~$$

And then using the approximations for $f'(x)$ and $f'(x+h)$ to do a further approximation.

Repeat this to find the approximate third derivative of this $f(x)$ at $x=1$ using $h=10^{-6}$. What value do you get?

Based on your answer, would you recommend this?

(If you were using Julia, you can make a function to find the approximate derivative with:

h = 1e-6
app_x(f) = x -> (f(x+h) -f(x))/h

app_x (generic function with 1 method)

Then finding the above is just:

app_xx(f) = x -> (app_x(f)(x + h) - app_x(f)(x))/h

app_xx (generic function with 1 method)

So, if $f(x) = \exp(x)\sin(x)$ and $x=1$ you have:

f(x) = exp(x) * sin(x)
app_xx(f)(1)

2.936761944738464

But –- you say – isn't this just doing the same thing twice? if so, can't it be done with composition? Yes, of course:

app_xx(f) = app_x(app_x(f))

app_xx (generic function with 1 method)

and we get the same answer as above:

app_xx(f)(1)

2.936761944738464

Except for notation that is satisfying. To address that, we can overload the postfix operator ' as follows:

Base.ctranspose(f::Function) = app_x(f)

ctranspose (generic function with 36 methods)

And then, the second derivative computed this way can be done with

f''(1)

2.936761944738464

To compare this with the exact answer at 1, we would have

app_xx(f)(1), D(f, 2)(1)

(2.936761944738464,2.9373878798317703)

That is the error is in the third decimal point. Maybe that isn't too bad.

What happens if you try this once more? Compare to the actual value.