First order derivatives

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Forward Difference¶

The forward difference is so named because it starts at $x$ and looks forward to obtain derivative information.

Consider the Taylor expansion of $f(x)$ from a point $x_i$ , with the step size $h$ :

\begin{aligned} f(x+h) &= f(x) + f'(x) h + f''(x) \frac{h^2}{2} + f'''(x) \frac{h^3}{6} + O(h^4) \end{aligned}

(1)

Truncating the series after the first two terms:

f(x+h) = f(x) + f'(x) h + O(h^2)

(2)

and solving for the first derivative yields the forward difference stencil:

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

(3)

which is to say that:

f'(x) \approx \frac{f(x+h) - f(x)}{h}

(4)

with a linear order of convergence since it depends on $O(h)$ .

Example: Rocket Velocity¶

A rocket has velocity:

v(t) = 2000 \ln \left(\frac{14 \times 10^4}{14 \times 10^4 - 2100t}\right) - 9.8 t

(5)

What is the acceleration as a function of $t$ ? Analytically we can find,

a(t) = v'(t) = \frac{42 \times 10^5}{14 \times 10^4 - 2100 t} - 9.8

(6)

def v(t):
    return 2000 * np.log((14 * 10**4) / (14 * 10**4 - 2100 * t)) - 9.8 * t

def a(t):
    return (42 * 10**5) / (14 * 10**4 - 2100 * t) - 9.8

t = np.linspace(0, 30, 100)
vt, at = v(t), a(t)

fig = make_subplots(rows=1, cols=2, subplot_titles=["Velocity v(t)", "Acceleration a(t)"])
fig.add_trace(go.Scatter(x=t, y=vt, name="v(t)"), row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=at, name="a(t)"), row=1, col=2)
fig.update_xaxes(title_text="Time (s)")
fig.update_yaxes(title_text="Velocity (m/s)", row=1, col=1)
fig.update_yaxes(title_text="Acceleration (m/sÂ²)", row=1, col=2)
fig.update_layout(showlegend=False, height=400, margin=dict(l=40, r=40, t=40, b=40))
fig.show()

Numerical Acceleration with Forward Difference¶

We implement the forward difference stencil and approximate the rocket acceleration across step sizes $h \in \{2, 1, 0.1\}$ :

def forward_difference(f, x, h):
    return (f(x + h) - f(x)) / h

a_fd_h2 = forward_difference(v, t, h=2)
a_fd_h1 = forward_difference(v, t, h=1)
a_fd_h0p1 = forward_difference(v, t, h=0.1)

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=a_fd_h2, mode="lines", name="FD h=2"))
fig.add_trace(go.Scatter(x=t, y=a_fd_h1, mode="lines", name="FD h=1"))
fig.add_trace(go.Scatter(x=t, y=a_fd_h0p1, mode="lines", name="FD h=0.1"))
fig.add_trace(go.Scatter(x=t, y=at, mode="lines", name="Analytic", line=dict(dash="dash", color="black")))
fig.update_layout(title="Forward Difference vs Analytic Acceleration", xaxis_title="Time (s)", yaxis_title="Acceleration (m/sÂ²)", height=420)
fig.show()

Backward Difference¶

By contrast, the backward difference steps backwards. Consider replacing $h$ with $-h$ :

f(x-h) = f(x) - f'(x) h + O(h^2)

(7)

and solving for the derivative:

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

(8)

def backward_difference(f, x, h):
    return (f(x) - f(x - h)) / h

a_bd_h2 = backward_difference(v, t, h=2)
a_bd_h1 = backward_difference(v, t, h=1)
a_bd_h0p1 = backward_difference(v, t, h=0.1)

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=a_bd_h2, mode="lines", name="BD h=2"))
fig.add_trace(go.Scatter(x=t, y=a_bd_h1, mode="lines", name="BD h=1"))
fig.add_trace(go.Scatter(x=t, y=a_bd_h0p1, mode="lines", name="BD h=0.1"))
fig.add_trace(go.Scatter(x=t, y=at, mode="lines", name="Analytic", line=dict(dash="dash", color="black")))
fig.update_layout(title="Backward Difference vs Analytic Acceleration", xaxis_title="Time (s)", yaxis_title="Acceleration (m/sÂ²)", height=420)
fig.show()

Central Difference¶

Note that the approximations are approaching the analytic solution from opposite directions. Can we combine forward and backward differences to get a better result?

Let’s examine the forward and backward expansions:

\begin{aligned} f(x+h) &= f(x) + f'(x) h + f''(x) \frac{h^2}{2} + f'''(x)\frac{h^3}{6} + O(h^4) \\ f(x-h) &= f(x) - f'(x) h + f''(x) \frac{h^2}{2} - f'''(x)\frac{h^3}{6} + O(h^4) \end{aligned}

(9)

Subtracting the two equations cancels the quadratic term:

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

(10)

Rearranging leads to the central difference formula:

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

(11)

which is second order accurate in $h$ since the next correction is $O(h^2)$ !

def central_difference(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

a_cd_h1 = central_difference(v, t, h=1)
a_cd_h2 = central_difference(v, t, h=2)

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=a_bd_h2, mode="lines", name="BD h=2"))
fig.add_trace(go.Scatter(x=t, y=a_fd_h2, mode="lines", name="FD h=2"))
fig.add_trace(go.Scatter(x=t, y=a_cd_h2, mode="lines", name="CD h=2"))
fig.add_trace(go.Scatter(x=t, y=at, mode="lines", name="Analytic", line=dict(dash="dash", color="black")))
fig.update_layout(title="Central Difference vs Forward and Backward", xaxis_title="Time (s)", yaxis_title="Acceleration (m/sÂ²)", height=420)
fig.show()

Error Comparison ( $h=2$ )¶

Let’'s compare the error in forward, backward, and central difference for $h=2$ across the first few time steps:

np.stack([a_fd_h2 - a(t), a_bd_h2 - a(t), a_cd_h2 - a(t)]).T[1:5, :]

Error Scaling: First-Order Schemes¶

Let’'s see how error scales when halving the step size from $h=2$ to $h''=1$ . For forward and backward difference, the error roughly cuts in half: $Error \sim \frac{h}{h''}$ , while the central difference cuts by $Error \sim \left[\frac{h}{h''}\right]^2$

print("Forward difference error ratio (h=2 vs h=1):")
print(np.stack([a_fd_h2 - a(t), a_fd_h1 - a(t)]).T[1:5, :])
print("\nBackward difference error ratio (h=2 vs h=1):")
print(np.stack([a_bd_h2 - a(t), a_bd_h1 - a(t)]).T[1:5, :])
print("\nCentral difference error ratio (h=2 vs h=1):")
np.stack([a_cd_h2 - a(t), a_cd_h1 - a(t)]).T[1:5, :]

Summary of First-Order Derivatives¶

The central difference algorithm:

Finds a more accurate solution in the same number of function calls.
Becomes more accurate with step size quadratically ( $O(h^2)$ ).
Requires information both before and after the evaluation point, meaning it cannot be used at boundary endpoints without one-sided stencils.
Requires equal step size!