using CalculusWithJulia
using Plots
plotly()
using ForwardDiff
using SymPy
using Roots
23 Numeric derivatives
This section uses these add-on packages:
SymPy
returns symbolic derivatives. Up to choices of simplification, these answers match those that would be derived by hand. This is useful when comparing with known answers and for seeing the structure of the answer. However, there are times we just want to work with the answer numerically. For that we have other options within Julia
. We discuss approximate derivatives and automatic derivatives. The latter will find wide usage in these notes.
23.0.1 Approximate derivatives
By approximating the limit of the secant line with a value for a small, but positive, \(h\), we get an approximation to the derivative. That is
\[ f'(x) \approx \frac{f(x+h) - f(x)}{h}. \]
This is the forward-difference approximation. The central difference approximation looks both ways:
\[ f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}. \]
Though in general they are different, they are both approximations. The central difference is usually more accurate for the same size \(h\). However, both are susceptible to round-off errors. The numerator is a subtraction of like-size numbers - a perfect opportunity to lose precision.
As such there is a balancing act:
- if \(h\) is too small the round-off errors are problematic,
- if \(h\) is too big the approximation to the limit is not good.
For the forward difference \(h\) values around \(10^{-8}\) are typically good, for the central difference, values around \(10^{-6}\) are typically good.
Example
Let’s verify that the forward difference isn’t too far off.
f(x) = exp(-x^2/2)
= 1
c = 1e-8
h = (f(c+h) - f(c)) / h fapprox
-0.6065306479285937
We can compare to the actual with:
@syms x
= diff(f(x), x)
df = convert(Float64, df(c))
factual abs(factual - fapprox)
1.1784039744000552e-8
The error is about \(1\) part in \(100\) million.
The central difference is better here:
= 1e-6
h = (f(c+h) - f(c-h)) / (2h)
cdapprox abs(factual - cdapprox)
1.5675682973892435e-11
The FiniteDifferences and FiniteDiff packages provide performant interfaces for differentiation based on finite differences.
23.0.2 Automatic derivatives
There are some other ways to compute derivatives numerically that give much more accuracy at the expense of slightly increased computing time. Automatic differentiation is the general name for a few different approaches. These approaches promise less complexity - in some cases - than symbolic derivatives and more accuracy than approximate derivatives; the accuracy is on the order of machine precision.
The ForwardDiff
package provides one of several ways for Julia
to compute automatic derivatives. ForwardDiff
is well suited for functions encountered in these notes, which depend on at most a few variables and output no more than a few values at once.
The ForwardDiff
package was loaded in this section; in general its features are available when the CalculusWithJulia
package is loaded, as that package provides a more convenient interface. The derivative
function is not exported by ForwardDiff
, so its usage requires qualification. To illustrate, to find the derivative of \(f(x)\) at a point we have this syntax:
derivative(f, c) # derivative is qualified by a module name ForwardDiff.
-0.6065306597126334
The CalculusWithJulia
package defines an operator D
which goes from finding a derivative at a point with ForwardDiff.derivative
to defining a function which evaluates the derivative at each point. It is defined along the lines of D(f) = x -> ForwardDiff.derivative(f,x)
in parallel to how the derivative operation for a function is defined mathematically from the definition for its value at a point.
Here we see the error in estimating \(f'(1)\):
= D(f)(c) # D(f) is a function, D(f)(c) is the function called on c
fauto abs(factual - fauto)
0.0
In this case, it is exact.
The D
operator is defined for most all functions in Julia
, though, like the diff
operator in SymPy
there are some for which it won’t work.
Example
For \(f(x) = \sqrt{1 + \sin(\cos(x))}\) compare the difference between the forward derivative with \(h=1e-8\) and that computed by D
at \(x=\pi/4\).
The forward derivative is found with:
f(x) = sqrt(1 + sin(cos(x)))
= pi/4, 1e-8
c, h = (f(c+h) - f(c))/h fwd
-0.20927346522370271
That given by D
is:
= D(f)(c)
ds_value - fwd ds_value, fwd, ds_value
(-0.20927346371432803, -0.20927346522370271, 1.5093746807970376e-9)
Finally, SymPy
gives an exact value we use to compare:
= diff(f(x), x) fp
\(- \frac{\sin{\left(x \right)} \cos{\left(\cos{\left(x \right)} \right)}}{2 \sqrt{\sin{\left(\cos{\left(x \right)} \right)} + 1}}\)
= convert(Float64, fp(PI/4))
actual - ds_value, actual - fwd actual
(-5.551115123125783e-17, 1.5093746252858864e-9)
Convenient notation
Julia
allows the possibility of extending functions to different types. Out of the box, the '
notation is not employed for functions, but is used for matrices. It is used in postfix position, as with A'
. We can define it to do the same thing as D
for functions and then, we can evaluate derivatives with the familiar f'(x)
. This is done in CalculusWithJulia
along the lines of Base.adjoint(f::Function) = D(f)
.
Then, we have, for example:
f(x) = sin(x)
'(pi), f''(pi) f
(-1.0, -1.2246467991473532e-16)
Example
Suppose our task is to find a zero of the second derivative of \(k(x) = e^{-x^2/2}\) in \([0, 10]\), a known bracket. The D
function takes a second argument to indicate the order of the derivative (e.g., D(f,2)
), but we use the more familiar notation:
k(x) = exp(-x^2/2)
find_zero(k'', 0..10)
1.0
We pass in the function object, k''
, and not the evaluated function.
23.1 Recap on derivatives in Julia
A quick summary of the \(3\) different ways for finding derivatives in Julia
presented in these notes:
- Symbolic derivatives are found using
diff
fromSymPy
- Automatic derivatives are found using the notation
f'
which utilizesForwardDiff.derivative
- approximate derivatives at a point,
c
, for a givenh
are found with(f(c+h)-f(c))/h
.
For example, here all three are computed and compared:
f(x) = exp(-x)*sin(x)
= pi
c = 1e-8
h
= diff(f(x),x)
fp
fp(c), f'(c), (f(c+h) - f(c))/h fp,
(-exp(-x)*sin(x) + exp(-x)*cos(x), -exp(-pi), -0.043213918263772265, -0.04321391703978389)
The use of '
to find derivatives provided by CalculusWithJulia
is convenient, and used extensively in these notes, but it needs to be noted that it does not conform with the generic meaning of '
within Julia
’s wider package ecosystem and may cause issue with linear algebra operations; the symbol is meant for the adjoint of a matrix.
23.2 Questions
Question
Find the derivative using a forward difference approximation of \(f(x) = x^x\) at the point \(x=2\) using h=0.1
:
Using D
or f'
find the value using automatic differentiation
Question
Mathematically, as the value of h
in the forward difference gets smaller the forward difference approximation gets better. On the computer, this is thwarted by floating point representation issues (in particular the error in subtracting two like-sized numbers in forming \(f(x+h)-f(x)\).)
For 1e-16
what is the error (in absolute value) in finding the forward difference approximation for the derivative of \(\sin(x)\) at \(x=0\)?
Repeat for \(x=\pi/4\):
Question
Let \(f(x) = x^x\). Using D
, find \(f'(3)\).
Question
Let \(f(x) = \lvert 1 - \sqrt{1 + x}\rvert\). Using D
, find \(f'(3)\).
Question
Let \(f(x) = e^{\sin(x)}\). Using D
, find \(f'(3)\).
Question
For Julia
’s airyai
function find a numeric derivative using the forward difference. For \(c=3\) and \(h=10^{-8}\) find the forward difference approximation to \(f'(3)\) for the airyai
function.
Question
Find the rate of change with respect to time of the function \(f(t)= 64 - 16t^2\) at \(t=1\).
Question
Find the rate of change with respect to height, \(h\), of \(f(h) = 32h^3 - 62 h + 12\) at \(h=2\).