Ch5.1: Finite Differences & Derivatives

• R code for finddif

findiff <- function(f,x,
               h = x*sqrt(.Machine$double.eps))
  {return((f(x + h)-f(x))/h)}

\[ f'(x) \cong \frac{f(x+h)-f(x)}{h} \]

• R code for symdiff

symdiff <- function(f,x,
             h = x*(.Machine$double.eps)^(1/3))
{return((f(x+h)-f(x-h))/(2*h))}

\[ f'(x) \cong \frac{f(x+h)-f(x-h)}{2h} \]

• How do we find the right step size h?
  • Small h leads to better accuracy
  • Too small leads to cancellation error
  • Too small leads to division by zero
• Trial & error with h values and see what is best?
• Use Calculus somehow to minimize error?

\[ \small{ \begin{aligned} f'(x) & \cong \frac{f(x+h)-f(x)}{h} \\ f'(x) &\cong \frac{f(x+h)-f(x-h)}{2h} \end{aligned} } \]

f <- function(x){x^2 + 1}; x <- 2

h <- c(3e-4,3e-8,3e-12,3e-16)
findiff(f,2,h)

[1] 4.000300 4.000000 3.999763 5.921189

(h <- x*sqrt(.Machine$double.eps))

[1] 2.980232e-08

findiff(f,x)

[1] 4

f <- function(x){x^2 + 1}; x <- 2

h <- c(1e-1,1e-5,1e-12,1e-16)
symdiff(f,x,h)

[1] 4.000000 4.000000 4.000356 0.000000

(h <- x*(.Machine$double.eps)^(1/3))

[1] 1.211091e-05

symdiff(f,x)

[1] 4

f <- function(x){sin(x)}; x <- 0.9

h <- c(5e-2,5e-6,5e-10,5e-14)
symdiff(f,x,h)

[1] 0.6213510 0.6216100 0.6216101 0.6217249

(h <- x*(.Machine$double.eps)^(1/3))

[1] 5.449909e-06

symdiff(f,x)

[1] 0.62161

• Find error bound for forward (backward) difference.

\[ \small{ f(x+h) = f(x) + f'(x)\left((x+h)-x)\right) + \frac{f''(x)}{2}\left((x+h)-x)\right)^2 + \ldots } \]

• It follows that

\[ \small{ \begin{aligned} f(x+h) - f(x) & = f'(x)h + \frac{f''(x)}{2}h^2 + \ldots \\ f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{f''(x)}{2}h + \ldots \\ f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{f''(c(x))}{2}h, \,\, \,\, c(x) \in [x,x+h] \\ \end{aligned} } \]

\[ \small{ \begin{aligned} f(x+h) & = f(x) + f'(x)\left((x+h)-x)\right) + \frac{f''(x)}{2}\left((x+h)-x)\right)^2 + \ldots \\ &= f(x) + f'(x)h + \frac{f''(x)}{2}h^2 + \frac{f'''(x)}{6}h^3 + \ldots \\ f(x-h) & = f(x) - f'(x)h + \frac{f''(x)}{2}h^2 - \frac{f'''(x)}{6}h^3 + \ldots \\ \end{aligned} } \]

• It follows that

\[ \small{ \begin{aligned} f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{f'''(x)}{6}h^3 + \ldots \\ f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c(x)), \,\, \,\, c(x) \in [x-h,x+h] \\ \end{aligned} } \]

• Using Taylor polynomials, we found error terms that go with the forward and centered difference formulas:

\[ \small{ \begin{aligned} f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f''(c), \,\,\, \, c(x) \in [x,x+h]\\ \\ f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c), \, \,\,\, c(x) \in [x-h,x+h] \end{aligned} } \]

• Truncation error for each is defined as

\[ \begin{aligned} E_T & = \frac{h}{2}f''(c(x)), \,\,\,\, c(x) \in [x,x+h] \\ E_T & = \frac{h^2}{6}f'''(c(x)), \,\,\,\, c(x) \in [x-h,x+h] \end{aligned} \]

• Error bounds are found as follows:

\[ \begin{aligned} |E_T| & \leq \frac{h}{2}M, \,\, M = \max_{u \in [x,x+h]}|f''(u)|, \\ |E_T| & \leq \frac{h^2}{6}M, \,\, M = \max_{u \in [x-h,x+h]}|f'''(u)| \end{aligned} \]

• Recall: How do we find the right step size h?
  • Small h leads to better accuracy
  • Too small leads to cancellation error
  • Too small leads to division by zero
• Trial & error with h values and see what is best?
• Use Calculus somehow to minimize error?

\[ \small{ \begin{aligned} f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f''(c(x)), \\ f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c(x)) \end{aligned} } \]

• For the forward difference results,

\[ \small{ f'(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f''(c), \,\, c(x) \in [x,x+h] } \]

• Residual (error) term from truncating infinite Taylor polynomial:

\[ \small{ E_T = \frac{h}{2}f''(c(x)) } \]

• Next, we look at roundoff error (see next slide).

• For roundoff error, denote the computed values by

\[ \begin{aligned} \tilde{f}(x+h) &= f(x+h) + e_1(x) \\ \tilde{f}(x) &= f(x) + e_2(x) \end{aligned} \]

• Then the computed forward difference formula is

\[ \small{ \begin{aligned} \tilde{f}'(x) &= \frac{\tilde{f}(x+h) - \tilde{f}(x)}{h} \\ &= \frac{f(x+h) - f(x) }{h} + \frac{e_1(x) - e_2(x)}{h} \\ &= \frac{f(x+h) - f(x) }{h} + E_R \end{aligned} } \]

• From previous slide,

\[ \begin{aligned} \tilde{f}'(x) = \frac{f(x+h) - f(x) }{h} + E_R \end{aligned} \]

• Here, we have denoted roundoff error by \( E_R \), with

\[ E_R = \frac{e_1(x) - e_2(x)}{h} \]

• Thus the total error \( E \) is given by

\[ E = E_R + E_T = \frac{e_1(x) - e_2(x)}{h} + \frac{h}{2}f''(c(x)) \]

• Our total error \( E \) is given by

\[ E = \frac{e_1(x) - e_2(x)}{h} + \frac{h}{2}f''(c) \]

• Observe that round off increases as \( h \) gets smaller, while truncation error decreases as \( h \) get smaller.

• We will use calculus to determine the optimal value of \( h \).

• The total error \( E \) is given by

\[ E = \frac{e_1(x) - e_2(x)}{h} + \frac{h}{2}f''(c) \]

• First, we observe that

\[ \left|e_1 - e_2\right| \leq |e_1| + |e_2| \leq 2 \max \left\{ |e_1|, |e_2| \right\} = 2e \]

• Also, let

\[ M = \max_{x \leq u \leq x+h} |f''(u)| \]

\[ \small{ \begin{aligned} f(x+h) = f(x) + f'(x)h & + \frac{f''(x)}{2}h^2 + \frac{f'''(x)}{6}h^3 + \frac{f^{(4)}(x)}{24}h^4 + \ldots \\ f(x-h) = f(x) - f'(x)h & + \frac{f''(x)}{2}h^2 - \frac{f'''(x)}{6}h^3 + \frac{f^{(4)}(x)}{24}h^4 + \ldots \end{aligned} } \]

• Add together and solve for \( f''(x) \):

\[ \small{ \begin{aligned} f(x+h) + f(x-h) & = 2f(x) + f''(x)h^2 + \frac{f^{(4)}(x)}{12}h^4 + \ldots \\ \\ f''(x) & = \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} - \frac{h^2}{12}f^{(4)}(c(x)) \end{aligned} } \]