How to take the numerical derivative of a function $f(x)$.
Solve for $f^{\prime}(x_i)$: $$ f^{\prime}(x_i) = \frac{f(x_{i+1})-f(x_i)}{\Delta x} - \left[\frac{1}{2}f^{\prime\prime}\Delta x^2 / \Delta x+\ldots\right].$$
$$ f^{\prime}(x_i) = \frac{f(x_{i+1})-f(x_i)}{\Delta x} + \mathcal{O}(\Delta x).$$Subtract the second equation from the first $\rightarrow$ $$ f_{i+1}-f_{i-1} = 2f_i^{\prime}\Delta x + \left[\frac{2}{6}f_i^{\prime\prime\prime}\Delta x^3+\ldots\right].$$
Solve for $f_i^{\prime}$: $$f_i^{\prime} = \frac{f_{i+1}-f_{i-1}}{2\Delta x} - \left[\frac{1}{6}f_i^{\prime\prime\prime}\Delta x^3+\ldots\right]/\Delta x.$$
Or,
$$f_i^{\prime} = \frac{f_{i+1}-f_{i-1}}{2\Delta x} + \mathcal{O}(\Delta x^2).$$Solve for $f_i^{\prime\prime}$
$$f_{i}^{\prime\prime} = \frac{f_{i-1}-2f_i+f_{i+1}}{\Delta x^2} - \mathcal{O}(\Delta x^2).$$For a term like $(-2)^2/2!$, the $(-2)^2$ is because we have $(-2\Delta x)^2 = (-2)^2\Delta x^2$ because we are $-2\Delta x$ away from $j=0$.
To take a linear combination, each column gets multiplied by some coefficient $c_j$ and we add these up.
We want the resulting coefficient of each term $f_0$, $f_0^{\prime}$, etc. to all be zero except the derivative we are approximating.
$\rightarrow$ the $\sum$ of each row, when we have columns times $c_j$ are 0, except the desired derivative, which is $1$.
Form a linear system: $Ac=b$ where $c$ are the coefficients of the linear combinations, $b^{T} = [0,1,0,0,0,0]$ for a first derivative, or $b^T = [0,0,1,0,0,0]$ for a second derivative, etc. The elements of $A$ are in the table above:
This gives $$f_0^{\prime}\approx\frac{1}{\Delta x}\left[-\frac{1}{2}f_{-1} + 0f_0 + \frac{1}{2}f_1\right] = \frac{f_1 - f_{-1}}{2\Delta x}.$$
This gives
$$f_0^{\prime}\approx\frac{1}{\Delta x}\left[-\frac{3}{2}f_0 + 2f_1 - \frac{1}{2}f_2\right] $$ * $\vec{R}_3^T\cdot \vec{c} = -1/3\ne 0$. * So, the order is $k-1 = 3-1 = 2$ or 2$^{nd}$ order. * $[0,\,1/6,\,4/3]\cdot[-3/2,\,2,\,-1/2] = -1/3$.Finite difference approximations to arbitrary order derivatives, using arbitrary stencils.
$$\frac{d^n u}{dx^n} = \frac{1}{\Delta x^n}\sum_k c_ku_k$$
As an example, a second order central second difference:
$$\frac{d^2 u}{dx^2} = \frac{1}{\Delta x^2}(c_{j-1} + c_ju_j + c_{j+1}u_{j+1}),$$
with $c$'s are 1, -2, 1
We assume here that $j=0$ is the base point.
oDer = 1 # order of the derivative to approximate (like 1 for f', or 2 for f", etc.)
npts = 9 # number of points in the stencil. Min is oDer + 1.
jb = 5 # position of the base point in the stencil, starting at 1
# the setup: oDer=1, npts=3, jb=2, will give a second order central difference approx to f'
#----------------------
jp = 1-jb:npts-jb
a_kj = zeros(npts,npts)
for k in 1:npts
a_kj[k,:] = jp.^(k-1) / factorial(k-1)
end
b = zeros(npts)
b[oDer+1] = 1
c_u = a_kj\b
kk = 0
for k ∈ oDer+1:100
aa = jp.^k/factorial(k)
if abs(aa'*c_u) > 1E-10 kk=k; break end
end
println("Grid_position (top), coefficient (bottom) =\n")
println(jp')
println(c_u')
println("Order = ", kk-oDer)
Grid_position (top), coefficient (bottom) = [-4 -3 -2 -1 0 1 2 3 4] [0.0035714285714285744 -0.03809523809523793 0.19999999999999885 -0.7999999999999965 -4.9960036108132076e-15 0.8000000000000043 -0.20000000000000207 0.038095238095238675 -0.003571428571428649] Order = 8
import Symbolics as sym
#----------------------
oDer = 1 # order of the derivative to approximate (like 1 for f', or 2 for f", etc.)
npts = 9 # number of points in the stencil. Min is oDer + 1.
jb = 5 # position of the base point in the stencil, starting at 1
# the setup: oDer=1, npts=3, jb=2, will give a second order central difference approx to f'
#----------------------
jp = 1-jb:npts-jb # index in grid: like [-2 -1 0 1 2 ] for jp=3 for npts=5
sym.@variables A[1:npts,1:npts] b[1:npts]
A = sym.scalarize(A)
b = sym.scalarize(b)
for k in 1:npts
A[k,:] .= jp.^(k-1) // factorial(k-1)
end
b[:] .= 0
b[oDer+1] = 1
c = A \ b
sym.@variables aa[1:npts]
aa = sym.scalarize(aa)
kk = 0
for k ∈ oDer+1:10
aa[:] .= jp.^k // factorial(k)
if aa'*c != 0 kk=k; break end
end
println("Order = ", kk-oDer)
println("Coefficients = ")
display(c)
#---------------------- output results (extra work outside Symbolics) to get the order
using Latexify
ulist = [sym.variable("u_$j") for j in jp]
uc = ulist .* c
str = "\$\$\\frac{du_0}{dx} = "
for (j,u) in enumerate(uc)
if c[j] == 0
str *= "0\\cdot " * latexify(string(ulist[j]))[2:end-1] * " + " # the [2:end-1] strips off the $'s
else
str *= latexify(string(u))[2:end-1] * " + "
end
end
str *= "\\mathcal{O}(\\Delta x^$(kk-oDer)) \$\$"
LaTeXString(str)
Order = 8 Coefficients =