Computers, Waves, Simulations
Finite-Difference Method - High-Order Taylor Operators

29. High-Order Taylor Operators#

This exercise covers the following aspects

  • Learn how to define high-order central finite-difference operators

  • Investigate the behaviour of the operators with increasing length

29.1. Basic Equations#

The Taylor expansion of f(x+dx) around x is defined as

f(x+dx)=n=0f(n)(x)n!(dx)n

Finite-difference operators can be calculated by seeking weights (here: a, b, c) with which function values have to be multiplied to obtain an interpolation or a derivative. Example:

af(x+dx)=a[f(x)+f(x)dx+12!f(x)(dx)2+]bf(x)=b[f(x)]cf(xdx)=c[f(x)f(x)dx+12!f(x)(dx)2]

This can be expressed in matrix form by comparing coefficients, here seeking a 2nd derivative

a+b+c=0ac=0a+c=2!dx2

which leads to

(111101101)(abc)=(002!(dx)2)

and using matrix inversion we obtain

(abc)=(12(dx)222(dx)212(dx)2)

This is the the well known 3-point operator for the 2nd derivative. This can easily be generalized to higher point operators and higher order derivatives. Below you will find a routine that initializes the system matrix and solves for the Taylor operator.

29.2. Calculating the Taylor operator#

The subroutine central_difference_coefficients() initializes the system matrix and solves for the difference weights assuming dx=1. It calculates the centered differences using an arbitrary number of coefficients, also for higher derivatives. The weights are defined at x±idx and i=0,..,(nop1)/2, where nop is the length of the operator. Careful! Because it is centered nop has to be an odd number (3,5,…)!

It returns a central finite difference stencil (a vector of length nop) for the nth derivative.

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import factorial
# Define function to calculate Taylor operators


def central_difference_coefficients(nop, n):
    """
    Calculate the central finite difference stencil for an arbitrary number
    of points and an arbitrary order derivative.

    :param nop: The number of points for the stencil. Must be
        an odd number.
    :param n: The derivative order. Must be a positive number.
    """
    m = np.zeros((nop, nop))
    for i in range(nop):
        for j in range(nop):
            dx = j - nop // 2
            m[i, j] = dx**i

    s = np.zeros(nop)
    s[n] = factorial(n)

    # The following statement return oper = inv(m) s
    oper = np.linalg.solve(m, s)
    # Calculate operator
    return oper

29.3. Plot Taylor operators#

Investigate graphically the Taylor operators. Increase nop for the first n=1 or higher order derivatives. Discuss the results and try to understand the interpolation operator (derivative order n=0).

# Calculate and plot Taylor operator


# Give length of operator (odd)
nop = 25
# Give order of derivative (0 - interpolation, 1 - first derivative, 2 - second derivative)
n = 1

# Get operator from routine 'central_difference_coefficients'
oper = central_difference_coefficients(nop, n)
# Plot operator
x = np.linspace(-(nop - 1) / 2, (nop - 1) / 2, nop)

# Simple plot with operator
plt.figure(figsize=(10, 4))
plt.plot(x, oper, lw=2, color="blue")
plt.plot(x, oper, lw=2, marker="o", color="blue")
plt.plot(0, 0, lw=2, marker="o", color="red")
# plt.plot (x, nder5-ader, label="Difference", lw=2, ls=":")
plt.title("Taylor Operator with nop =  %i " % nop)
plt.xlabel("x")
plt.ylabel("Operator")
plt.grid()
../_images/aa34e16efc1213a7f88dedf5b11fdb9c3d83b72636714b65ed878acf1e8e26ed.png

29.4. Conclusions#

  • The Taylor operator weights decrease rapidly with distance from the central point

  • In practice often 4th order operators are used to calculate space derivatives