9 - Numerical Differentiation

Differentiation

The Problem

• Given $f : \mathbb{R} \to \mathbb{R}$, evaluate $f'(x)$ for a given input $x$
• Similar to integration, finding the derivative analytically is expensive for some functions

เราอาจจะไม่ได้ต้องการว่า $f’(x)$ เป็นอะไร แต่อยากรู้ว่า $f’(2)$ คือเท่าไหร่ ประมาณนี้แหละ!

Goal

• Approximate $f'(x)$ for a given input $x$ numerically without computing the formula for $f'$

Analogy: Think of this like estimating the slope of a hill at a specific point by measuring how much you climb over a small distance, rather than using calculus formulas.

Finite Difference

Forward Difference

Let $f : \mathbb{R} \to \mathbb{R}$. By Taylor series:
$f(x + h) = f(x) + f'(x)h + \frac{f''(x)}{2!}h^2 + \frac{f'''(x)}{3!}h^3 + \ldots \quad (1)$
Rearranging:
$f'(x)h = f(x + h) - f(x) - \frac{f''(x)}{2!}h^2 - \frac{f'''(x)}{3!}h^3 - \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} + O(h)$

เป็น $O(h)$ ก็เพราะว่าหลังจากก็จะเป็น Degree of $h$ ไปเรื่อย ๆ

Forward difference formula:
$\boxed{f'(x) \approx \frac{f(x + h) - f(x)}{h}}$

• This has first-order accuracy (error is $O(h)$)

$h$ should be very small and finite number!

Backward Difference

By Taylor series:

$f(x - h) = f(x) - f'(x)h + \frac{f''(x)}{2!}h^2 - \frac{f'''(x)}{3!}h^3 + \ldots \quad (2)$
Rearranging:

$f'(x)h = f(x) - f(x - h) + \frac{f''(x)}{2!}h^2 - \frac{f'''(x)}{3!}h^3 + \ldots$

$f'(x) = \frac{f(x) - f(x - h)}{h} + \frac{f''(x)}{2}h + \ldots$

$f'(x) = \frac{f(x) - f(x - h)}{h} + O(h)$

Backward difference formula:
$\boxed{f'(x) \approx \frac{f(x) - f(x - h)}{h}}$

• Also first-order accurate (error is $O(h)$)

Centered Difference

Subtract equation (2) from equation (1): จากด้านบนนะ

$f(x + h) - f(x - h) = 2f'(x)h + 2\frac{f'''(x)}{3!}h^3 + \ldots$

$2f'(x)h = f(x + h) - f(x - h) - \frac{f'''(x)}{3}h^3 - \ldots$

$f'(x) = \frac{f(x + h) - f(x - h)}{2h} - \frac{f'''(x)h^2}{6} + \ldots$

$f'(x) = \frac{f(x + h) - f(x - h)}{2h} + O(h^2)$

Centered difference formula:
$\boxed{f'(x) \approx \frac{f(x + h) - f(x - h)}{2h}}$

• Second-order accurate (error is $O(h^2)$)

Analogy: This is like measuring the slope by looking both forward and backward and taking the average - more balanced and accurate!

Centered Difference for Second Derivative

Add equations (1) and (2):
$f(x + h) + f(x - h) = 2f(x) + f''(x)h^2 + 2\frac{f^{(4)}(x)}{4!}h^4 + \ldots$
$f''(x)h^2 = f(x + h) + f(x - h) - 2f(x) - \frac{f^{(4)}(x)}{12}h^4 - \ldots$
$f''(x) = \frac{f(x + h) - 2f(x) + f(x - h)}{h^2} - \frac{f^{(4)}(x)}{12}h^2 - \ldots$
$f''(x) = \frac{f(x + h) - 2f(x) + f(x - h)}{h^2} + O(h^2)$
Centered difference for $f''$:

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

• Second-order accurate

Example 1

Problem: Approximate $f'(2)$ where $f(x) = 2x^4 - 3x^3$ with the three methods. Use $h = 0.1$.

Solution:

Forward difference:
$f'(2) \approx \frac{f(2 + 0.1) - f(2)}{0.1} = \frac{(2(2.1)^4 - 3(2.1)^3) - (2(2^4) - 3(2^3))}{0.1} \approx 31.1320$

Backward difference:
$f'(2) \approx \frac{f(2) - f(2 - 0.1)}{0.1} = \frac{(2(2^4) - 3(2^3)) - (2(1.9)^4 - 3(1.9)^3)}{0.1} \approx 25.1280$

Centered difference:
$f'(2) \approx \frac{f(2 + 0.1) - f(2 - 0.1)}{2(0.1)} = \frac{(2(2.1)^4 - 3(2.1)^3) - (2(1.9)^4 - 3(1.9)^3)}{0.2} \approx 28.1300$

Note:

• $f'(x) = 8x^3 - 9x^2$
• True value: $f'(2) = 8(2^3) - 9(2^2) = 28$
• Centered difference is closest!

Second derivative approximation: Approximate $f''(1)$. Use $h = 0.1$.

$f''(1) \approx \frac{f(1.1) - 2f(1) + f(0.9)}{(0.1)^2}$

$= \frac{(2(1.1)^4 - 3(1.1)^3) - 2(2(1)^4 - 3(1)^3) + (2(.9)^5 - 3(.9)^3)}{.01} \approx 6.0400$

Note:

• $f''(x) = 24x^2 - 18x$
• True value: $f''(1) = 24 - 18 = 6$

Exercise 2

1. Approximate $f'(-1)$ where $f(x) = x - 2x^2$ with the three methods. Use $h = 0.2$
2. Approximate $f''(-1)$. Use $h = 0.2$

Solution

Calculate function values:

• $f(-1 + 0.2) = f(-0.8) = (-0.8) - 2(-0.8)^2 = -2.08$
• $f(-1) = (-1) - 2(-1)^2 = -3$
• $f(-1 - 0.2) = f(-1.2) = (-1.2) - 2(-1.2)^2 = -4.08$

Forward difference:
$f'(-1) \approx \frac{f(-0.8) - f(-1)}{0.2} = 4.6$

Backward difference:
$f'(-1) \approx \frac{f(-1) - f(-1.2)}{0.2} = 5.4$

Centered difference:
$f'(-1) \approx \frac{f(-0.8) - f(-1.2)}{0.4} = 5$

Centered difference for $f''$:
$f''(-1) \approx \frac{f(-0.8) - 2f(-1) + f(-1.2)}{(0.2)^2} = -4$

Example 3

Problem: Use Taylor series to derive a second-order accurate approximation to $f'(x)$ in terms of the values of $f(x)$, $f(x - 2h)$, and $f(x + 4h)$.

Given Taylor series:
$f(x - 2h) = f(x) - 2f'(x)h + 2f''(x)h^2 - \frac{4f'''(x)h^3}{3} + O(h^4) \quad (3)$
$f(x + 4h) = f(x) + 4f'(x)h + 8f''(x)h^2 + \frac{32f'''(x)h^3}{3} + O(h^4) \quad (4)$

Solution

To get a second-order accurate approximation, we need to eliminate the $h^2$ term.

Multiply (3) by 4:

$4f(x - 2h) = 4f(x) - 8f'(x)h + 8f''(x)h^2 - \frac{16f'''(x)h^3}{3} + O(h^4) \quad (5)$

Subtract (5) from (4):

$f(x + 4h) - 4f(x - 2h) = -3f(x) + 12f'(x)h + 16f'''(x)h^3 + O(h^4)$

$12f'(x)h = f(x + 4h) - 4f(x - 2h) + 3f(x) - 16f'''(x)h^3 + O(h^4)$

$f'(x) = \frac{f(x + 4h) - 4f(x - 2h) + 3f(x)}{12h} - \frac{4}{3}f'''(x)h^2 + O(h^3)$

Final formula:

$\boxed{f'(x) \approx \frac{f(x + 4h) - 4f(x - 2h) + 3f(x)}{12h}}$

• This is second-order accurate as the error term is $O(h^2)$

Alternative Method for Deriving Finite Difference Formulas

Key Idea

• The previous Taylor-series method is cumbersome for higher-order or higher-derivative methods
• Alternative approach: Interpolate $f$ and find the derivative of the interpolant

Setup

• Let $f : \mathbb{R} \to \mathbb{R}$ be the function whose derivatives we want to find
• Let $t_i$, $i = 1, \ldots, n$ be equally spaced points in $\mathbb{R}$, with step size $h = t_{i+1} - t_i$
• Let $y_i = f(t_i)$

Two-Point Interpolation (Forward Difference)

Let's interpolate two points $(t_i, y_i)$ and $(t_{i+1}, y_{i+1})$.

Lagrange interpolant:
$p(t) = y_i\frac{t - t_{i+1}}{t_i - t_{i+1}} + y_{i+1}\frac{t - t_i}{t_{i+1} - t_i}$
$= y_i\frac{t - t_{i+1}}{-h} + y_{i+1}\frac{t - t_i}{h}$
$= y_{i+1}\frac{t - t_i}{h} - y_i\frac{t - t_{i+1}}{h}$

Take the derivative:

$p'(t) = \frac{y_{i+1}}{h} - \frac{y_i}{h} = \frac{y_{i+1} - y_i}{h} = \frac{f(t_{i+1}) - f(t_i)}{h} = \frac{f(t_i + h) - f(t_i)}{h}$

This gives the forward difference formula for $f'$.

Two-Point Interpolation (Backward Difference)

Now interpolate $(t_{i-1}, y_{i-1})$ and $(t_i, y_i)$.

Lagrange interpolant:

$p(t) = y_{i-1}\frac{t - t_i}{t_{i-1} - t_i} + y_i\frac{t - t_{i-1}}{t_i - t_{i-1}}$

$= y_{i-1}\frac{t - t_i}{-h} + y_i\frac{t - t_{i-1}}{h}$

$= y_i\frac{t - t_{i-1}}{h} - y_{i-1}\frac{t - t_i}{h}$

Take the derivative:

$p'(t) = \frac{y_i}{h} - \frac{y_{i-1}}{h} = \frac{y_i - y_{i-1}}{h} = \frac{f(t_i) - f(t_{i-1})}{h} = \frac{f(t_i) - f(t_i - h)}{h}$

This gives the backward difference formula for $f'$.

Three-Point Interpolation (Centered Difference)

Interpolate three data points: $(t_{i-1}, y_{i-1})$, $(t_i, y_i)$, $(t_{i+1}, y_{i+1})$.

Lagrange interpolant:

$p(t) = y_{i-1}\frac{(t - t_i)(t - t_{i+1})}{(t_{i-1} - t_i)(t_{i-1} - t_{i+1})} + y_i\frac{(t - t_{i-1})(t - t_{i+1})}{(t_i - t_{i-1})(t_i - t_{i+1})}$

$+ y_{i+1}\frac{(t - t_{i-1})(t - t_i)}{(t_{i+1} - t_{i-1})(t_{i+1} - t_i)}$

Substituting $h = t_{i+1} - t_i$:

$p(t) = y_{i-1}\frac{(t - t_i)(t - t_{i+1})}{(-h)(-2h)} + y_i\frac{(t - t_{i-1})(t - t_{i+1})}{h(-h)}$

$+ y_{i+1}\frac{(t - t_{i-1})(t - t_i)}{2h(h)}$

$= y_{i-1}\frac{(t - t_i)(t - t_{i+1})}{2h^2} - y_i\frac{(t - t_{i-1})(t - t_{i+1})}{h^2} + y_{i+1}\frac{(t - t_{i-1})(t - t_i)}{2h^2}$

Using the product rule ($\frac{d}{dt}(uv) = \frac{du}{dt}v + u\frac{dv}{dt}$):

$p'(t) = y_{i-1}\frac{(t - t_i) + (t - t_{i+1})}{2h^2} - y_i\frac{(t - t_{i-1}) + (t - t_{i+1})}{h^2}$

$+ y_{i+1}\frac{(t - t_{i-1}) + (t - t_i)}{2h^2} \quad (6)$

Evaluate at $t = t_i$:

$p'(t_i) = y_{i-1}\frac{(t_i - t_i) + (t_i - t_{i+1})}{2h^2} - y_i\frac{(t_i - t_{i-1}) + (t_i - t_{i+1})}{h^2}$

$+ y_{i+1}\frac{(t_i - t_{i-1}) + (t_i - t_i)}{2h^2}$

$= y_{i-1}\frac{-h}{2h^2} - y_i\frac{h - h}{h^2} + y_{i+1}\frac{h}{2h^2}$

$= -\frac{y_{i-1}}{2h} + \frac{y_{i+1}}{2h} = \frac{y_{i+1} - y_{i-1}}{2h}$

$= \frac{f(t_{i+1}) - f(t_{i-1})}{2h} = \frac{f(t_i + h) - f(t_i - h)}{2h}$

This gives the centered difference formula for $f'$.

Second Derivative

Take the derivative of equation (6) with respect to $t$:
$p''(t) = y_{i-1}\frac{2}{2h^2} - y_i\frac{2}{h^2} + y_{i+1}\frac{2}{2h^2}$
$= \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}$
$= \frac{f(t_{i+1}) - 2f(t_i) + f(t_{i-1})}{h^2}$
$= \frac{f(t_i + h) - 2f(t_i) + f(t_i - h)}{h^2}$

This gives the centered difference formula for $f''$.

General Approach for Higher-Order Methods

To derive higher-order methods or methods for higher derivatives:

1. Add more data points, e.g., $(t_{i-2}, y_{i-2})$, $(t_{i+2}, y_{i+2})$, $(t_{i-3}, y_{i-3})$, $(t_{i+3}, y_{i+3})$, etc.
2. Find the interpolant
3. Take the derivative of the interpolant

Analogy: It's like fitting a smoother curve through more points to get a better estimate of the slope.

Step Size Considerations

Key Points

• In many problems, including numerical integration or differentiation, smaller step sizes give more accurate results
  • For quadrature: think of the widths of the subintervals for composite rule as step sizes
• But step sizes cannot be too small

Issues with Very Small Step Sizes

If step size is too small:

• Too expensive (e.g., too many function evaluations for quadrature)
• Too high error (e.g., catastrophic cancellation for finite difference)

Analogy: Using a microscope with too much magnification - you see so much detail that you lose perspective, and tiny measurement errors become huge problems.

Richardson Extrapolation

Another idea to improve accuracy for finite difference.

Setup

• Let $F(h)$ denote the value obtained from a given method with step size $h$

  • Example: the value obtained from a finite difference method using step size $h$

• Suppose:
  $F(h) = a_0 + a_1h^p + O(h^r)$
  as $h \to 0$ for some $p$ and $r$, with $r > p$

• We know the values of $p$ (and $r$), but not $a_0$ or $a_1$

  • ตอนทำจริงไม่จำเป็นต้องรู้ค่าแท้จริงของ $r$ ก็ได้ แค่รู้ว่ามัน $r>p$ ก็พอ เพราะสุดท้าย ใน Richardson Extrapolation Formula ไม่มี $r$ อยู่ซักตัวในสมการ5555

• $a_0$ is the true value we want to find (the true derivative, for instance)

• $p$ is the order of the formula $F$

มันก็เหมือน Forward Difference นั่นแหละ ที่ $p=1,r=2$

Derivation

Suppose we compute $F$ for two step sizes: $h$ and $h/q$ for some positive integer $q$.

We have:
$F(h) = a_0 + a_1h^p + O(h^r) \quad (7)$

and
$F(h/q) = a_0 + a_1\left(\frac{h}{q}\right)^p + O(h^r)$

$F(h/q) = a_0 + a_1q^{-p}h^p + O(h^r) \quad (8)$

This is a system of two linear equations in two unknowns $a_0$ and $a_1$.

Solving for $a_0$:

Multiply (7) by $q^{-p}$:
$q^{-p}F(h) = q^{-p}a_0 + a_1q^{-p}h^p + O(h^r)$
Subtract (8):
$q^{-p}F(h) - F(h/q) = (q^{-p} - 1)a_0 + O(h^r)$

$(q^{-p} - 1)a_0 = q^{-p}F(h) - F(h/q) - O(h^r)$

$(q^{-p} - 1)a_0 = q^{-p}F(h) - F(h/q) + O(h^r)$

$a_0 = \frac{q^{-p}F(h) - F(h/q)}{q^{-p} - 1} + O(h^r)$
Rearranging:
$a_0 = \frac{q^{-p}F(h) - F(h) + F(h) - F(h/q)}{q^{-p} - 1} + O(h^r)$

$a_0 = \frac{(q^{-p} - 1)F(h) + F(h) - F(h/q)}{q^{-p} - 1} + O(h^r)$

$a_0 = F(h) + \frac{F(h) - F(h/q)}{q^{-p} - 1} + O(h^r)$

Richardson Extrapolation Formula

$\boxed{a_0 = F(h) + \frac{F(h) - F(h/q)}{q^{-p} - 1}}$
Key observations:

• The above has error $O(h^r)$, which is smaller than the errors of $F(h)$ and $F(h/q)$
• If $F(h)$ is known for several values of $h$, can repeat the extrapolation process to get even more accurate approximations (Not going through that detail na)

Analogy: It's like taking two rough measurements and using the pattern in their errors to calculate a much more accurate estimate.

Example 4

Problem: Use Richardson extrapolation to improve the accuracy of the first-order forward difference formula to approximate the derivative of $\sin(x)$ at $x = 1$. Use $h = 0.5$ and $h/2 = 0.25$.

Solution

Note: for forward difference, $p = 1$ and $r = 2$.

$F(h) = F(0.5) = \frac{\sin(1.5) - \sin(1)}{0.5} = 0.312048$
$F(h/2) = F(0.25) = \frac{\sin(1.25) - \sin(1)}{0.25} = 0.430055$

The extrapolated value is:

$a_0 = F(h) + \frac{F(h) - F(h/q)}{q^{-p} - 1}$

$= 0.312048 + \frac{0.312048 - 0.430055}{2^{-1} - 1}$

$= 0.312048 + \frac{0.312048 - 0.430055}{(1/2) - 1}$

$= 0.548061$
For comparison:

• The correctly rounded result: $\cos(1) = 0.540302$
• Much closer than either original approximation!

Automatic Differentiation (AD) Preview

จะสอนแบบไม่ได้เจาะลึกมากเพราะยาก

What is AD?

• A tool/software that takes a program code that computes $f$ and outputs a program code that computes $f'$ (and $f$)
• The output code computes the true derivatives (subject to roundoff error)
  • In contrast, finite difference formulas have truncation error too
• Works for any code regardless of how complicated $f$ is
  • (Although the output derivatives code may be very slow or very expensive spacewise)

Analogy: It's like having a robot that can watch you bake a cake and automatically write down the recipe for making the cake rise twice as fast.

The Idea Behind AD

Core concept:

• Any computer program is just a sequence of basic arithmetic operations (addition, subtraction, multiplication, and division), whose derivatives are easy to compute
• Use chain rules (properly) on each statement of the code for $f$ to get the code for $f'$

Two modes:

• Forward mode
• Reverse mode
• They differ in how the chain rule is used

How Forward-mode AD Works

Original code:

function y = fun(x)
    u = 2*x+3;
    v = u*(x+2);
    y = u+v;
end

Forward mode approach:

• For each statement $w = \ldots$, add a statement to compute $dw/dx$ immediately after it
• For simplicity, let $\dot{z}$ denote $dz/dx$
• Note: $\dot{z}$ is a new variable in the program

AD-generated code (forward mode):

function [y, ẏ] = fun(x)
    u = 2*x+3;
    u̇ = 2;
    v = u*(x+2);
    v̇ = u̇*(x+2) + u;
    y = u+v;
    ẏ = u̇+v̇;
end

Analogy: Forward mode is like watching a domino chain fall and calculating how fast each domino hits the next one as you go forward.

Reverse-mode AD [Optional]

Approach:

1. Run the original code once first
2. From bottom to top, for each statement $w = g(t_1, \ldots, t_n)$, add a statement that computes $\frac{\partial y}{\partial w}\frac{dg}{dt_i}$ and adds it to $dy/dt_i$ for all $t_i$

Analogy: Reverse mode is like tracing backward from the final result to see how much each earlier step contributed to it - like forensic analysis of a chain reaction.

Forward Mode vs Reverse Mode

When to use each:

• Suppose $f : \mathbb{R}^n \to \mathbb{R}^m$
• Forward mode is faster when $m$ is much larger than $n$
• Reverse mode is faster when $n$ is much larger than $m$
• Reverse mode takes a lot more space (except for some special cases)

Practical note:

• For machine learning (many inputs, one output loss function): reverse mode (backpropagation) is typically used
• For simulations (few inputs, many outputs): forward mode is typically better

Summary

Key Formulas

First Derivative:

• Forward difference (1st order): $f'(x) \approx \frac{f(x+h) - f(x)}{h}$
• Backward difference (1st order): $f'(x) \approx \frac{f(x) - f(x-h)}{h}$
• Centered difference (2nd order): $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$

Second Derivative:

• Centered difference (2nd order): $f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}$

Richardson Extrapolation:
$a_0 = F(h) + \frac{F(h) - F(h/q)}{q^{-p} - 1}$

Key Concepts

• Higher order = more accurate (but requires more function evaluations)
• Centered difference is generally more accurate than forward/backward
• Step size trade-off: smaller is more accurate, but too small causes numerical issues
• Richardson extrapolation can boost accuracy without reducing step size
• Automatic differentiation computes exact derivatives by applying calculus rules to code