Fourier characteristics of repeated numerical derivative

Question

Background

I am trying to analyse fourier characteristics of a derivative. For example if I have a first order derivative approximated as following: $$\frac{\partial \Psi(x)}{\partial x} = \frac{\Psi_{i+1}-\Psi_{i-1}}{\Delta x}$$

to get its fourier characteristics, substitute $\Psi(x)= e^{j\, kx}$ in LHS, and its discrete counter part $\Psi(x=i\Delta x) = \Psi_{i}= e^{j \,ki\Delta x}$ RHS. Here $j=\sqrt{-1}$ in and k represents the wavenumber. It gives

$$F \left[\frac{\partial \Psi(x)}{\partial x}\right] = \frac{e^{jk(i+1)\Delta x}-e^{jk(i-1)\Delta x}}{\Delta x}$$ Comparing both sides gives $$k= \sin(j \, k\Delta x)$$

Problem statement

I like to carry out similar analysis for implicit derivatives of lele (1992), given as following $$f_i' + \alpha (f_{i+1}' + f_{i-1}') + \beta (f_{i+2}' + f_{i-2}') \nonumber \\ =a\frac{f_{i+1}- f_{i-1}}{(2h)^2} + b\frac{f_{i+2}- f_{i-2}}{(4h)^2} + c\frac{f_{i+3}- f_{i-3}}{(6h)^2} $$ The coefficients are determined after Taylor series expansion of both sides and then comparing coefficients.

My problem contains a derivative term of type $\frac{\partial}{\partial x} \left[a(x) \frac{\partial}{\partial x} \right]\psi $. It can be simplified if I assume constant media for Fourier analysis:

$$\frac{\partial}{\partial x} \left[a(x) \frac{\partial}{\partial x} \right]\psi = a(x) \frac{\partial}{\partial x} \frac{\partial\psi}{\partial x} \psi= a(x) \frac{\partial^2\psi}{\partial x^2} $$

In simulation I am using first derivative two times for computing the second derivatives as following:

1. step computes : $\zeta = a(x) \frac{\partial}{\partial x} $
2. step compute : $\frac{\partial \zeta}{\partial x}$

For explicit derivatives (example given in background part) the equality holds i.e. $\partial_x.\partial_x ==\partial_x^2$ is valid in numerical form.

So my question is

1. whether $\partial_x.\partial_x$ is equal to $ \partial_x^2$ for implicit derivative.
2. If it doesn't then how can I carry out its analysis.
3. How to handle the cross derivative terms e.g. $\partial_x.\partial_y$.

Any help or suggestion is appreciated!

Answer 1

Regarding point 1: generally, applying a discretized derivative (stencil) twice is not equivalent to applying an equivalent (i.e. same number of points) stencil for the second derivative.

Example in matrix representation: A three point stencil for the first derivative would correspond to a tridiagonal matrix. The same is true for the three point stencil for the second derivative. Taking the square of the matrix representing the first derivative leads to a matrix where also the second sub-diagonal is non-zero (the corresponding five-point form is different from the five-point stencil obtained via Taylor expansion and thus, unlike the latter, does not have the (generally optimal) $O\!\left({\rm gridspacing}^4\right)$ error for the second derivative).

The reason for squaring these discretized representations of the first derivative not being equivalent to the direct discretization of the second derivative is that the corresponding basis set is incomplete. Only in the limit of infinitely small grid spacing between the stencil points will the square of the approximated first derivative equal the second derivative approximation. Infinitely small grid spacing also corresponds to a complete basis for (the considered interval of) $\mathbb{R}$.

An exception to this incomplete basis set problem are Chebyshev polynomial basis sets: here the second derivative is equivalent to squaring the first derivative operator in the Chebyshev basis.

Regarding point 3: Since $a$ depends on $x$, I don't follow the steps presented (assuming constant media would mean $a$ is a constant and not $a(x)$?). I would transform according to the product rule for differentiation: $$ \frac{\partial}{\partial x}\left[a(x)\frac{\partial}{\partial x}\right] \psi = \frac{\partial a(x)}{\partial x} \frac{\partial \psi}{\partial x} + a(x) \frac{\partial^2 \psi}{\partial x^2} $$

The representations of $\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$ can be applied in any order (the corresponding matrices commute). Care only has to be taken in the case of $x$ and $y$ being the same variable as discussed above.

Here is a specific example for a differential operator in real space representation:

In real space, the three point stencil would correspond to the first derivative operator $$D_{ij} = -\frac{1}{2} \delta_{i,j+1} + \frac{1}{2} \delta_{i,j-1}.$$ Taking the square $$D_{ij}^2 = \frac{1}{4} \delta_{i,j+2} - \frac{1}{2} \delta_{i,j} + \frac{1}{4} \delta_{i,j-2}.$$ This is different from (and a worse approximation than) starting from Taylor expansion and using the three point stencil for the second derivative directly: $$D^{\rm second}_{ij} = \delta_{i,j+1} - 2 \delta_{i,j} + \delta_{i,j-1}.$$ Was the grid spacing infinitely small (meaning making the basis set complete), $D$ would go over to the analytical differential operator, and applying $D$ twice and second derivative would be identical.

Thus, generally, if the basis set in which $D$ is represented is finite (or really incomplete), repeated application of $D$ is not equivalent to creating higher differential operators with the same technique as used for $D$ in this finite basis set. In a complete basis set (infinitely small grid spacing limit as an example), repeated application of $D$ is equivalent to higher order derivative operators.

As a final example, there can also be finite basis sets, in which powers of $D$ are identical to higher derivatives. In Fourier space, e.g., derivatives simply become multiplication with the frequency variable $\omega$ and thus repeated multiplication by $\omega$ is identical to higher derivatives, no matter how many values of $\omega$ are included in the finite basis.