Backwards finite difference for mixed partials at higher order

Question

I'm trying to understand what is the general method for calculating a backwards difference for a mixed partial of $n$ variables. Let's start with one variable:

The forward and backward finite differences are

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

and

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

However, when we go to mixed partials for $n=2$ variables discretized on a 2D grid with $h = x_{i+1}-x_{i}$, we can use (a central difference).

$$\left(\frac{\partial^2 u}{\partial x\partial y} \right)_{i,j} = \frac{u_{i+1,j+1} - u_{i+1,j-1} - u_{i-1,j+1} + u_{i-1,j-1}}{4 \Delta x \Delta y} + \mathscr{O}\left((\Delta x)^2, (\Delta y)^2 \right) \qquad (1)$$

How would one go about finding $\mathscr{O}\left((\Delta x)^2, (\Delta y)^2 \right)$? Furthermore, what about backwards difference? Notice the above central difference needs information of the $i+1,j+1$ point which is problematic near or approaching a singularity.

EDIT: Is this approached with a matrix method of the Taylor expansion of the mixed partials?

Answer 1

You can write the Taylor series for the function at the points used in (1), e.g.:

$$ f(x+\Delta x, y+ \Delta y) = \sum_{i=0}^{\infty} \dfrac{(\Delta x)^i}{i!} \dfrac{\partial^{i}}{\partial x^i} f(x, y+ \Delta y) = \sum_{i=0}^{\infty} \dfrac{(\Delta x)^i}{i!} \sum_{j=0}^{\infty} \dfrac{(\Delta y)^j}{j!} \dfrac{\partial^{i+j}}{\partial y^j \partial x^i} f(x, y) $$

, and use them to find the formal accuracy of the scheme.

Answer 2

There are multiple ways to find finite difference coefficients for any arbitrary order of accuracy. Let us try to achieve second order accuracy to maintain simplicity. I have listed down three ways which could be used:

1. Method I: Use Lagrange polynomials. Assume that the solution u can be approximated by a Lagrange polynomial within both the x and y directions. Use three points to generate the interpolant in both x and y directions to achieve second order accuracy.
2. Method II: Assume that the derivative can be approximated by a weighted sum of points where the weights are undetermined. To achieve second order accuracy, you might need 3 or more points in both directions. Verify truncation error is O($\Delta x^2, \Delta y^2$) in the remaining Taylor expansion.
3. Method III: Use discrete differential operators. Apply these operators one after the other but note that this method does not guarantee second order accuracy even if the respective operators are individually second order accurate. As a result, one needs to verify the order of accuracy by performing Taylor expansion of the terms generated after performing the differencing.

I am going to use Method II first and then show how the same result could be obtained more easily via a simple application of Method III. I am using central difference in the x direction and backward difference in the y direction.

Let us assume $(u_{x,y})_{i,j} = a_1 u_{i+1, j} + a_2 u_{i-1, j} + a_3 u_{i+1, j-1} + a_4u_{i-1, j-1} + a_5 u_{i+1, j-2} + a_6 u_{i-1, j-2}$ where $a_1, a_2, \ldots, a_6 \in R$ are to be determined.

Now consider a Taylor expansion of each one of the $u$ terms in the equation above assuming constant spacing $h$ in both $x$ and $y$ directions.

$u_{x,y} = a_1( u_{i,j} + h(u_x)_{i,j} + (h^2/2)(u_{xx})_{i,j} + (h^3/6)(u_{xxx})_{i,j} + (h^4/24)(u_{xxxx})_{i,j} + \ldots) + a_2( u_{i,j} - h(u_x)_{i,j} + (h^2/2)(u_{xx})_{i,j} - (h^3/6)(u_{xxx})_{i,j} + (h^4/24)(u_{xxxx})_{i,j} + \ldots) + a_3( u_{i,j} + h(u_x)_{i,j} -h(u_y)_{i,j} + (h^2/2)(u_{xx})_{i,j} - (h^2)(u_{xy})_{i,j}+ (h^2/2)(u_{yy})_{i,j} + (h^3/6)(u_{xxx})_{i,j} - (h^3/6)(u_{yyy})_{i,j} - (h^3/2)(u_{xxy})_{i,j} + (h^3/2)(u_{xyy})_{i,j} + (h^4/24)(u_{xxxx})_{i,j} - (h^4/6)(u_{xxxy})_{i,j} + (h^4/4)(u_{xxyy})_{i,j} - (h^4/6)(u_{xyyy})_{i,j} + (h^4/24)(u_{yyyy})_{i,j} \ldots ) +a_4( u_{i,j} - h(u_x)_{i,j} - h(u_y)_{i,j} + (h^2/2)(u_{xx})_{i,j} + (h^2)(u_{xy})_{i,j} + (h^2/2)(u_{yy})_{i,j} - (h^3/6)(u_{xxx})_{i,j} - (h^3/6)(u_{yyy})_{i,j} - (h^3/2)(u_{xxy})_{i,j} - (h^3/2)(u_{xyy})_{i,j} + (h^4/24)(u_{xxxx})_{i,j} + (h^4/6)(u_{xxxy})_{i,j} + (h^4/4)(u_{xxyy})_{i,j} + (h^4/6)(u_{xyyy})_{i,j} + (h^4/24)(u_{yyyy})_{i,j} \ldots ) +a_5( u_{i,j} + h(u_x)_{i,j} - 2h(u_y)_{i,j} + (h^2/2)(u_{xx})_{i,j} - (2h^2)(u_{xy})_{i,j} + (2h^2)(u_{yy})_{i,j} + (h^3/6)(u_{xxx})_{i,j} - (4h^3/3)(u_{yyy})_{i,j} - (h^3)(u_{xxy})_{i,j} + (2h^3)(u_{xyy})_{i,j} + (h^4/24)(u_{xxxx})_{i,j} - (h^4/3)(u_{xxxy})_{i,j} + (h^4)(u_{xxyy})_{i,j} - (4h^4/3)(u_{xyyy})_{i,j} + (2h^4/3)(u_{yyyy})_{i,j} \ldots ) +a_6( u_{i,j} - h(u_x)_{i,j} - 2h(u_y)_{i,j} + (h^2/2)(u_{xx})_{i,j} + (2h^2)(u_{xy})_{i,j} + (2h^2)(u_{yy})_{i,j} - (h^3/6)(u_{xxx})_{i,j} - (4h^3/3)(u_{yyy})_{i,j} - (h^3)(u_{xxy})_{i,j} - (2h^3)(u_{xyy})_{i,j} + (h^4/24)(u_{xxxx})_{i,j} + (h^4/3)(u_{xxxy})_{i,j} + (h^4)(u_{xxyy})_{i,j} + (4h^4/3)(u_{xyyy})_{i,j} + (2h^4/3)(u_{yyyy})_{i,j} \ldots )$

I have used approximations until fourth order for the above expansion to make sure that the error terms remains second order accurate after division by $h^2$. Now, equating both sides of the equation for the terms $u, u_{x}, u_{y}, u_{xx}, \ldots$, we get a total of six equations for six variables:

$a_1 + a_2 + a_3 + a_4 + a_5 + a_6 = 0$
$a_1 - a_2 + a_3 - a_4 + a_5 - a_6 = 0$
$- a_3 - a_4 - 2a_5 - 2a_6 = 0$
$- a_3 + a_4 - 2a_5 + 2a_6 = 1/h^2$
$a_3/2 + a_4/2 + 2a_5 + 2a_6 = 0$
$-a_3/6 - a_4/6 - 4a_5/6 - 4a_6/6 = 0$

One can solve for all the a's and make some simplifications to obtain the below approximation for $(u_{xy})_{i,j}$

$ (u_{xy})_{i,j} = ( (3u_{i+1,j} -4u_{i+1, j-1} + u_{i+1, j-2}) - (3u_{i-1,j} -4u_{i-1, j-1} + u_{i-1, j-2}) )/(4h^2) $

One can verify from the above expansion that the error term is $(h^2/6) u_{xxxy} - (h^2/3)u_{xyyy}$ which is $O(h^2)$.

Note that this approximation could also be obtained after successive application of the discrete differential operators. Let $\Delta_x$ and $\Delta_y$ denote the discrete operators in the x and y directions respectively.

Then $u_{xy} \approx \Delta_y \Delta_x u = \Delta_y(u_{i + 1,j} - u_{i-1,j})/(2h) = (\Delta_y (u_{i + 1,j}) - \Delta_y (u_{i-1,j}) )/(2h) = ( (3u_{i+1,j} -4u_{i+1, j-1} + u_{i+1, j-2}) - (3u_{i-1,j} -4u_{i-1, j-1} + u_{i-1, j-2}) )/(4h^2)$

In this case, it turns out that the discrete operators result in the same equation as the one above and hence, achieve second order accuracy. Any other type of differencing in y or x can also be done in the same manner as described above.