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.
import numpy as np
#----------------------
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 = np.arange(1-jb, npts-jb+1) # index in grid: like [-2 -1 0 1 2 ] for jp=3 for npts=5
a_kj = np.zeros((npts,npts))
j = np.arange(npts)
for k in range(npts) :
a_kj[k,j] = jp[j]**k / np.math.factorial(k)
b = np.zeros(npts)
b[oDer] = 1
c_u = np.linalg.solve(a_kj,b)
# check the order
aa = np.zeros(npts)
for k in range(oDer+1,100) :
aa[j] = (jp[j]**k)/np.math.factorial(k)
if np.abs(np.sum(aa*c_u)) > 1.0E-10 :
break
#----------------------
# output the results
print('Grid_position (top), coefficient (bottom) =\n')
print(np.vstack([jp, c_u]))
print('\nOrder = ', k-oDer)
Grid_position (top), coefficient (bottom) = [[-4.00000000e+00 -3.00000000e+00 -2.00000000e+00 -1.00000000e+00 0.00000000e+00 1.00000000e+00 2.00000000e+00 3.00000000e+00 4.00000000e+00] [ 3.57142857e-03 -3.80952381e-02 2.00000000e-01 -8.00000000e-01 -4.99600361e-15 8.00000000e-01 -2.00000000e-01 3.80952381e-02 -3.57142857e-03]] Order = 8
import sympy as sp
sp.init_printing(long_frac_ratio=1)
import numpy as np
from sympy.printing.latex import LatexPrinter
from IPython.display import Markdown as md # hack to pretty output the latex...
#----------------------
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 = np.arange(1-jb, npts-jb+1, dtype=np.int32) # index in grid: like [-2 -1 0 1 2 ] for jp=3 for npts=5
A = sp.zeros(npts,npts)
for k in range(npts):
for j in range(npts):
A[k,j] = jp[j]**k / sp.factorial(k)
B = sp.zeros(npts,1)
B[oDer] = 1
C = A**-1 * B
c = np.array(C.T).astype(np.float64)[0] # numpy array of coefficients
#---------------------- check the order
j = np.arange(npts)
for k in range(oDer+1,100) :
aa = (jp[j]**k)/np.math.factorial(k)
if np.abs(np.sum(aa*c)) > 1.0E-10 :
break
#---------------------- output results
points = [sp.symbols(f'u_{j}') for j in jp]
Δx, u, x, = sp.symbols('Δx, u, x')
u0 = sp.symbols('u_0')
rr = []
for i in range(npts):
rr.append(C[i]*points[i])
rhs = sp.Add(*rr,evaluate=False) # keep the order of the additive terms...
lhs = LatexPrinter(dict(order='none'))._print(sp.Derivative(u0,x,oDer))
rhs = LatexPrinter(dict(order='none'))._print_Add(rhs)
expr = lhs + "=" + rhs + f"+O(\Delta x^{k-oDer})"
md(f"$${expr}$$")
C