Discretized scalar fields¶

A time-dependent scalar field $c({\bf x, t})$ with ${\bf x}\in \mathbb{R}^d$, $t \in (0, \infty)$ in $d$-dimensional space, can be approximated by a discretized scalar field of the form

(1)\[\begin{split} c({\bf x}) & \approx & \sum_{i=1}^{N} c_{i}(t) \psi_i({\bf x}),\end{split}\]

where $C_{i}$ is the value of the scalar field at the $i$ th point, $N = \prod_{i=1}^{d} n_i$, $\,n_i$ is the number of lattice points in the $i$ th dimension, $d$ is the dimension of the space, and $\psi_i$ are interpolation functions. We assume that the node points are uniformly spaced with grid-step $h$ the same in all dimensions, and use the notation ${\bf x^{(n)}}=(x_{1}^{(n)},x_{2}^{(n)},\ldots, x_{d}^{(n)})$ to designate the $n$ th lattice point.

Three-dimensional space: derivatives at lattice points¶

We will make use of the Taylor expansion for functions with three-dimensional spatial dependence

(2)\[\begin{split} f( {\bf x} + {\bf h} ) & \approx & f( {\bf x} ) + \sum_{i=1}^3 \frac{ \partial f({\bf x})}{\partial x_i} h + \frac{1}{2!} \sum_{i=1}^3 \sum_{j=1}^3 \frac{ \partial^2 f({\bf x}) }{\partial x_i \partial x_j} h^2 + \mathcal{O}(h^3).\end{split}\]

We will also make use of the notation ${\bf h}_i \equiv ( \delta_{i1}, \delta_{i2}, \delta_{i3})h$, where $\delta_{ij}$ is the Kronecker delta

(3)\[\begin{split} \delta_{ij} & = & \left\{ \begin{array}{rr} 0 & i \ne j \\ 1 & i = j \end{array} \right. .\end{split}\]

Gradient¶

To approximate first derivative in the $x_1$-direction at the lattice points, we use ${\bf h}=(\pm h,0,0)$ with the Taylor expansions (2). Then

(4)\[\begin{split} f({\bf x}^{(n)} + {\bf h}) & \approx & f({\bf x}^{(n)}) +\frac{\partial f({\bf x}^{(n)})}{\partial x_1} h + \frac{1}{2!} \frac{\partial^2 f({\bf x}^{(n)})}{\partial x_1^2 } h^2 + \mathcal{O}(h_1^3) \\ f({\bf x}^{(n)} - {\bf h}) & \approx & f({\bf x}^{(n)}) - \frac{\partial f({\bf x}^{(n)})}{\partial x_1} h + \frac{1}{2!} \frac{\partial^2 f({\bf x}^{(n)})}{\partial x_1^2 } h^2 + \mathcal{O}(h_1^3) .\end{split}\]

Subtracting the two equations above gives the $x_1$-derivative to first order at the $n\text{th}$ lattice point

(5)\[\begin{split} \frac{\partial f({\bf x}^{(n)})}{\partial x_1} & \approx & \frac{ f({\bf x}^{(n)}+h) - f({\bf x}^{(n)}-h_1) }{2 h} + \mathcal{O}(h^2) .\end{split}\]

Generalizing to other dimensions and using (5) and (3), the gradient operator to first order (in a rectilinear coordinate-system) is

(6)\[ \nabla f({\bf x}^{(n)}) = \sum_{i=1}^3 \frac{\partial f({\bf x}^{(n)})}{\partial x_i} {\bf\hat{x}_i} \approx \sum_{i=1}^3 \frac{1}{2 h} \left[ f( {\bf x}^{(n)} + {\bf h}_i ) - f( {\bf x}^{(n)} - {\bf h}_i ) \right] {\bf\hat{x}_i},\]

where ${\bf\hat{x}_i}$ is the unit vector in the $i\text{th}$ dimension.

Gradient at boundary nodes¶

We only consider boundary surfaces whose unit normal are coincident with the rectilinear axes of the space.

The central differences in (6) cannot be evaluated at the boundaries. To approximate the gradient at the boundary, we use a polynomial approximation of the function to get a first order estimate of the derivative in the direction perpendicular to the boundary surface. For example, in the vicinity of the boundary point $(x_1=0, x_2, x_3)$, we approximate the function as $f(x_1,x_2,x_3) \approx c_0 + c_1 x_1 + c_2 x_1^2$. Then,

(7)\[\begin{split} \renewcommand{\arraystretch}{2} \begin{array}{rcl} f(0,x_2,x_3) & \approx & c_0 \\ f(h,x_2,x_3) & \approx & c_0 + c_1 h + c_2 h^2 \\ f(2h,x_2,x_3) & \approx & c_0 +2 c_1 h + 4 c_2 h^2 \end{array}\end{split}\]

Solving for the coefficients $c_0$ and $c_1$ gives

(8)\[\begin{split} \renewcommand{\arraystretch}{2} \begin{array}{rcl} c_0 & = & f(0,x_2,x_3) \\ c_1 & = & (2h)^{-1} \left[ - 3f(0,x_2,x_3)+ 4f(h,x_2,x_3) - f(2h,x_2,x_3) \right] \end{array}\end{split}\]

Differentiating the polynomial and evaluating at $x=0$ gives

(9)\[ \left. \frac{df(x_1,x_2,x_3)}{dx_1} \right|_{(0,x_2,x_3)} = c_1 = (2h)^{-1} \left[ - 3f(0)+ 4f(h,x_2,x_3) - f(2h,x_2,x_3) \right],\]

which approximates the derivative in the direction perpendicular to the boundary in terms of the values at the boundary node and neighbor nodes.

Similarly, in the vicinity of the boundary point $(x_1=L, x_2, x_3)$ we approximate the function as $f(x_1,x_2,x_3) \approx c_0 + c_1 x_1 + c_2 x_1^2$. Then,

(10)\[\begin{split} \renewcommand{\arraystretch}{2} \begin{array}{rcl} f(L,x_2,x_3) & \approx & c_0 + c_1 L + c_2 L^2 \\ f(L-h,x_2,x_3) & \approx & c_0 + c_1 (L-h) + c_2 (L-h)^2 \\ f(L-2h,x_2,x_3) & \approx & c_0 + c_1 (L-2h) + c_2 (L-2h)^2 \end{array}\end{split}\]

From the first equation we obtain $c_0 = f(L,x_2,x_3) - c_1 L - c_2 L^2$. Multiplying the second equation by 4, subtracting the third equation and using the above result, gives

(11)\[\begin{split} c_1 + 2 L c_2 & = & (2h)^{-1} \left[ 3 f(L,x_2,x_3) - 4 f(L-h,x_2,x_3) + f(L-2h,x_2,x_3) \right]\end{split}\]

Differentiating the polynomial at $x=L$ and using the result above gives

(12)\[ \left. \frac{df(x_1,x_2,x_3)}{dx_1} \right|_{(L,x_2,x_3)} = c_1 + 2Lc_2 = (2h)^{-1} \left[ 3 f(L,x_2,x_3) - 4 f(L-h,x_2,x_3) + f(L-2h,x_2,x_3) \right].\]

Approximation of the derivatives at other boundaries follow in a similar manner.

Laplacian¶

To approximate second derivatives, we add the two equations in (4). This gives the second derivative in the $x_1$-direction (to first order)

(13)\[\begin{split} \frac{ \partial f^2( {\bf x}^{(n)} )}{\partial x_1^2} & \approx & \frac{1}{ h^2} \left( f( {\bf x}^{(n)} + {\bf h_1} ) + f( {\bf x}^{(n)} - {\bf h_1} ) - 2f( {\bf x}^{(n)} ) \right) + \mathcal{O}(h^2)\end{split}\]

where again we are using the notation ${\bf h}_i \equiv ( \delta_{i1}, \delta_{i2}, \delta_{i3})h$.

Generalizing to other dimensions and using (5) and (3), the Laplacian operator to first order (in a rectilinear coordinate-system) is

(14)\[ \nabla^2 f( {\bf x}^{(n)} ) = \sum_{i=1}^3 \frac{\partial^2 f( {\bf x}^{(n)} )}{\partial x_i^2} \approx \sum_{i=1}^3 \frac{1}{h_i^2} \left[ f( {\bf x}^{(n)} + {\bf h}_i ) + f( {\bf x}^{(n)} - {\bf h}_i ) - 2f( {\bf x}^{(n)} )\right] .\]

Diffusion¶

Consider a scalar field $f({\bf x},t)$ describing the time-dependent concentration of molecules in volume $V$ with surface $S$. The time evolution of $f({\bf x},t)$ due to diffusion is given by

(15)\[\begin{split} \frac{\partial f({\bf x},t)}{\partial t} & = & \zeta \nabla^2 f({\bf x},t),\end{split}\]

subject to boundary conditions

(16)\[ {\bf \nabla} f({\bf x}_s) \cdot {\bf \hat{n}} = - \zeta^{-1}\; \phi({\bf x}_s) \;\; \mbox{ and/or } \;\; f({\bf x}_s) = \mbox{constant}\]

where $\zeta$ is the diffusion coefficient, $\phi({\bf x})$ is the flux (amount of material per unit area per time) through a surface element at position ${\bf x}_s$ on $S$, and ${\bf \hat{n}}$ is the unit normal to the surface, which points out of the volume.

Using (14), the Laplacian applied to lattice points on the interior is, to first order, given by

(17)\[\begin{split} \frac{\partial f({\bf x^{(n)}},t)}{\partial t} & \approx & \zeta \sum_{i=1}^3 \frac{1}{h_i^2} \left[ f( {\bf x^{(n)}} + {\bf h}_i ) + f( {\bf x^{(n)}} - {\bf h}_i ) - 2f({\bf x^{(n)}})\right] ,\end{split}\]

where again ${\bf h}_i \equiv ( \delta_{i1}, \delta_{i2}, \delta_{i3}) h$ and $\delta_{ij}$ is the Kronecker delta.

Boundary conditions¶

Flux (Neumann)¶

When we impose flux boundary conditions, we specify the flux ${\bf \phi}({\bf x}_s)$ at the surface point ${\bf x}_s$. Then from (17), this imposes a gradient at the surface

(18)\[ \nabla f({\bf x}_s,t) \cdot {\bf \hat{n}} = \sum_{i=1}^{3} \frac{d f({\bf x}_s, t) }{d x_i} n_i = - \frac{\phi({\bf x}_s)}{\zeta} .\]

Flux at natural boundaries

In general, the discretized points of the volume may not be coincident with the discretized points of the surface. Such cases will occur for curved surfaces or for surfaces moving with respect to the volume, such as with cells. We define natural boundaries to be those boundaries for which the lattice points of the surface are coincident with lattice points of the volume.

Consider a lattice point ${\bf x_s}^{(n)}$ on the surface of a volume. Using (17) and (6), the flux condition on the gradient at ${\bf x_s}^{(n)}$ can be written as

(19)\[ \nabla f({\bf x_s}^{(n)},t) \cdot {\bf \hat{n}} = \sum_{i=1}^3 \frac{ n_i}{2 h_i} \left[ f( {\bf x_s}^{(n)} + {\bf h}_i , t) - f( {\bf x_s}^{(n)} - {\bf h}_i, t ) \right] = - \frac{\phi({\bf x_s}^{(n)})}{\zeta} ,\]

where again ${\bf h}_i \equiv ( \delta_{i1}, \delta_{i2}, \delta_{i3}) h$. For ${\bf x_s}^{(n)}$ on the surface, 1 to 3 of the evaluation points ${\bf x}_s^{(n)} \pm {\bf h}_i$ lie outside the volume. We impose the flux by calculating the value of the field at these {it virtual (image)} lattice points as function of the applied flux (using (19)), then substitute these resulting expressions for the virtual points in the application of the Laplacian at the surface point (14).

Let us illustrate with some simple cases.

Consider the scalar field $f({\bf x}, t)$ defined in the volume $V$, which is a rectangular prism with extents $0 \le x_1 \le L_1$, $0 \le x_2 \le L_2$, and $0 \le x_3 \le L_3$, with flux $\phi$ specified at the natural boundaries.

Example 1: Consider the point $x_1=L_1$ on the boundary , with ${\bf x_{1L}} \equiv (L_1,x_2,x_3)$ and ${\bf \hat{n}} = (1, 0,0)$. Then using (19), we can compute the value of the virtual lattice point $f(L_1+h_1,x_2,x_3)$

(20)\[\begin{split} \frac{ 1}{2 h_1} \left[ f(L_1+h_1,x_2,x_3) - f(L_1-h_1,x_2,x_3) \right] & = & - \frac{\phi( L_1,x_2,x_3 )}{\zeta} \\ f(L_1+h_1,x_2,x_3) & = & - \frac{2 h_1 \phi( L_1,x_2,x_3 ) }{\zeta} + f(L_1-h_1,x_2,x_3) .\end{split}\]

Substituting this in the equation for the Laplacian (14) gives the second derivative in the $x_1$ direction at the boundary point as a function of the applied flux

(21)\[\begin{split} \frac{ \partial f^2( {\bf x_{1L} } )}{\partial x_1^2} & = & \frac{1 }{h_1^2} \left[ f(L_1+h_1,x_2,x_3) + f(L_1-h_1,x_2,x_3) - 2 f(L_1,x_2,x_3) \right] \\ & = & \frac{1}{h_1^2} \left[ - \frac{2 h_1 \phi(L_1,x_2,x_3 ) }{\zeta} + 2 f(L_1-h_1,x_2,x_3) - 2 f(L_1,x_2,x_3) \right] \\ & = & \frac{2}{h_1^2} \left[ f(L_1-h_1,x_2,x_3) - f(L_1,x_2,x_3) \right] - \frac{2 \phi(L_1,x_2,x_3) }{h_1 \zeta} .\end{split}\]

Example 2: Consider the point on the boundary $x_1=0$, with ${\bf x_{10}} \equiv (0,x_2,x_3)$ and ${\bf \hat{n}} = (-1, 0,0)$. Using (19), we can compute the value of the virtual lattice point $f(-h_1,x_2,x_3)$

(22)\[\begin{split} - \frac{ 1}{2 h_1} \left[ f(h_1,x_2,x_3) - f(-h_1,x_2,x_3) \right] & = & - \frac{\phi( 0,x_2,x_3 )}{\zeta} \\ f(-h_1,x_2,x_3) & = & - \frac{2 h_1 \phi(0,x_2,x_3) }{\zeta} + f(h_1,x_2,x_3) .\end{split}\]

Substituting this in the equation for the Laplacian (14) gives the second derivative in the $x_1$ direction at the boundary point as a function of the applied flux

(23)\[\begin{split} \frac{ \partial f^2( {\bf x_{10} } )}{\partial x_1^2} & = & \frac{1 }{h_1^2} \left[ f(h_1,x_2,x_3) + f(-h_1,x_2,x_3) - 2 f(0,x_2,x_3) \right] \\ & = & \frac{1}{h_1^2} \left[ - \frac{2 h_1 \phi(0,x_2,x_3 ) }{\zeta} + 2 f(h_1,x_2,x_3) - 2 f(0,x_2,x_3) \right] \\ & = & \frac{2}{h_1^2} \left[ f(h_1,x_2,x_3) - f(0,x_2,x_3) \right] - \frac{2 \phi(0,x_2,x_3) }{h_1 \zeta} .\end{split}\]

Flux from a tiny sphere¶

Consider a spherical surface $S$ with flux $\phi$ at an arbitrary location ${\bf x_s}$ in the volume $V$. (In this case, the surface is not a natural boundary.) If the radius $\rho$ of the sphere is much less than the distance $h$ between lattice points $\rho \ll h$, then the flux acts like a point source and the rate of change of the concentration at ${\bf x}_s$ is give by

(24)\[ \frac{ \partial f({\bf x_s})}{\partial t} = \frac{4\pi\rho^2}{h^3} \phi .\]

We use interpolation to distribute the change in concentration to the lattice points of the voxel that contains the point:math:{bf x}_s.

Dirichlet (fixed)¶

A Dirichlet boundary conditions specifies a fixed value $\psi({\bf x}_s)$ for the function at the boundary point ${\bf x}_s$. Therefore, the value of the function at the boundary does not change in time and

(25)\[ \frac{\partial f({\bf x}_s,t)}{\partial t} = \zeta \nabla^2 f({\bf x}_s,t) = 0 .\]

Toroidal¶

Essentially, there are no boundaries for the toroidal case; material that diffuses through one face enters through the opposite face.

Example 1: Consider the point $x_1=L_1$ on the boundary, with ${\bf x_{1L}} \equiv (L_1,x_2,x_3)$. Applying the toroidal boundary condition $f(0,x_2,x_3) = f(L_1,x_2,x_3)$ to the second derivative in the $x_1$-direction (13) gives

(26)\[\begin{split} \frac{ \partial f^2( {\bf x_{1L} } )}{\partial x_1^2} & = & \frac{1 }{h_1^2} \left[ f(L_1+h_1,x_2,x_3) + f(L_1-h_1,x_2,x_3) - 2 f(L_1,x_2,x_3) \right] \\ & = & \frac{1 }{h_1^2} \left[ f(h_1,x_2,x_3) + f(L_1-h_1,x_2,x_3) - 2 f(L_1,x_2,x_3) \right]\end{split}\]

Example 2: Consider the point $x_1=0$ on the boundary, with ${\bf x_{10}} \equiv (0,x_2,x_3)$. Applying the toroidal boundary condition $f(0,x_2,x_3) = f(L_1,x_2,x_3)$ to the second derivative in the $x_1$-direction (13) gives

(27)\[\begin{split} \frac{ \partial f^2( {\bf x_{10} } )}{\partial x_1^2} & = & \frac{1 }{h_1^2} \left[ f(h_1,x_2,x_3) + f(-h_1,x_2,x_3) - 2 f(0,x_2,x_3) \right] \\ & = & \frac{1 }{h_1^2} \left[ f(h_1,x_2,x_3) + f(L_1-h_1,x_2,x_3) - 2 f(0,x_2,x_3) \right]\end{split}\]

Bilinear interpolation¶

We will use the following notation in this section.

grid-spacing $h$
$f^n \equiv f({\bf x^{(n)}})$ is the value of the function at the $n\text{th}$ node.
$f^n_x$ ($f^n_y$) is the value of the $x$-derivative ($y$-derivative )at node $n$ (point ${\bf x^{(n)}}$).

$math/img/pixel_grad_interpolation.png$

Figure 1: Node designations for interpolation.

Scalar field interpolation¶

Point $f^0$ lies in the pixel with surrounding nodes: $f^1$, $f^2$, $f^3$, and $f^4$ (Figure 1). We use a 2D Taylor expansions about point $f^0$ to approximate the function at the surrounding nodes

(28)\[\begin{split} f^1 & = & f^0 - \delta_x f^0_x - \delta_y f^0_y + \mathcal{O}(h^2) \\ f^2 & = & f^0 - \delta_x f^0_x + (h - \delta_y) f^0_y + \mathcal{O}(h^2) \\ f^3 & = & f^0 + (h - \delta_x) f^0_x - \delta_y f^0_y + \mathcal{O}(h^2) \\ f^4 & = & f^0 + (h - \delta_x) f^0_x + (h - \delta_y) f^0_y + \mathcal{O}(h^2) .\end{split}\]

Combining the first two and last two equations in (28) gives

(29)\[\begin{split} (h-\delta_y) f^1 + \delta_y f^2 & = & h f^0 - h \delta_x f^0_x + \mathcal{O}(h^3) \\ (h-\delta_y) f^3 + \delta_y f^4 & = & h f^0 + h (h - \delta_x) f^0_x + \mathcal{O}(h^3)\end{split}\]

Combining the equations in (29) gives

(30)\[\begin{split} (h-\delta_x) \left[ (h-\delta_y) f^1 + \delta_y f^2 \right] + \delta_x \left[ (h-\delta_y) f^3 + \delta_y f^4 \right] & = & h^2 f^0 + \mathcal{O}(h^4) .\end{split}\]

Finally,

(31)\[\begin{split} f^0 & \approx & (1-\Delta_x)(1-\Delta_y) f^1 + (1-\Delta_x) \Delta_y f^2 + \Delta_x (1-\Delta_y) f^3 + \Delta_x \Delta_y f^4 + \mathcal{O}(h^2) ,\end{split}\]

where $\Delta \equiv \frac{\delta}{h}$.

Equation (31) can be easily generalized to the 3D case.

Gradient interpolation¶

Applying the bilinear interpolation equation (31) to the $x$-derivative $f_x^0$ at point $x^0$ gives

(32)\[\begin{split} f_x^0 & = & (1-\Delta_x)(1-\Delta_y) f_x^1 + (1-\Delta_x) \Delta_y f_x^2 + \Delta_x (1-\Delta_y) f_x^3 + \Delta_x \Delta_y f_x^4 + \mathcal{O}(h^2) .\end{split}\]

Using the first-order approximations for the derivatives at the node points (5)

(33)\[\begin{split} f^1_x & \approx & \frac{1}{2h} \left( f^{3} - f^{9} \right)\\ & \vdots &\end{split}\]

in (32) gives

(34)\[\begin{split} f_x^0 & = & \frac{1}{2h} \left[ (1-\Delta_x)(1-\Delta_y) \left( f^{3} - f^{9} \right) + (1-\Delta_x) \Delta_y \left( f^{4} - f^{10} \right) \right. \\ & & \; \; \; \left. + \; \Delta_x (1-\Delta_y) \left( f^{11} - f^{1} \right) + \Delta_x \Delta_y \left( f^{12} - f^{2} \right) \right] .\end{split}\]

The derivatives calculated with (34) are continuous across nodes.

Base node $(i,j)$.

(35)\[\begin{split} f_x^0 & = & \frac{1}{2h} \left[ (1-\Delta_x)(1-\Delta_y) \left( f(i+1,j) - f(i-1,j) \right) \right. \\ & & \; \; \; \left. \right] .\end{split}\]

Equation (34) can easily be generalized to derivatives in other dimensions and to the 3D case.

Gradient interpolation at boundaries¶

$math/img/pixel_grad_interpolation_boundary.png$

Figure 2: Node designations for interpolation.

Case 1: Consider the case of gradient interpolation at point $f^0$, which lies in the pixel with surrounding nodes 1,2, 3, an 4, with associated function values $f^1$, $f^2$, $f^3$, and $f^4$ (Figure 2 (a)), where nodes 1 and 2 lie on the boundary.

Applying the bilinear interpolation equation (31) to the $x$-derivative $f_x^0$ at point $x^0$ gives

(36)\[\begin{split} f_x^0 & \approx & (1-\Delta_x)(1-\Delta_y) f_x^1 + (1-\Delta_x) \Delta_y f_x^2 + \Delta_x (1-\Delta_y) f_x^3 + \Delta_x \Delta_y f_x^4 .\end{split}\]

Using the first-order approximations for the derivatives at the node points (5) and (12)

(37)\[\begin{split} f^1_x & \approx & \frac{1}{2h} \left[ - 3 f^{1} + 4 f^{3} - f^{11} \right] \\ f^2_x & \approx & \frac{1}{2h} \left[ - 3 f^{2} + 4 f^{4} - f^{12} \right] \\ f^3_x & \approx & \frac{1}{2h} \left( f^{11} - f^{1} \right) \\ f^4_x & \approx & \frac{1}{2h} \left( f^{12} - f^{2} \right) \\\end{split}\]

in (36) gives

(38)\[\begin{split} f_x^0 & \approx & (1-\Delta_x)(1-\Delta_y) \frac{1}{2h} \left[ - 3 f^{1} + 4 f^{3} - f^{11} \right] + (1-\Delta_x) \Delta_y \frac{1}{2h} \left[ - 3 f^{2} + 4 f^{4} - f^{12} \right] \\ & & + \Delta_x (1-\Delta_y) \frac{1}{2h} \left( f^{11} - f^{1} \right) + \Delta_x \Delta_y \frac{1}{2h} \left( f^{12} - f^{2} \right) .\end{split}\]

Case 2: Consider the case of gradient interpolation at point $f^0$, which lies in the pixel with surrounding nodes 1,2, 3, an 4, with associated function values $f^1$, $f^2$, $f^3$, and $f^4$ (Figure 2 (b)), where nodes 3 and 4 lie on the boundary.

Applying the bilinear interpolation equation (31) to the $x$-derivative $f_x^0$ at point $x^0$ gives

(39)\[\begin{split} f_x^0 & \approx & (1-\Delta_x)(1-\Delta_y) f_x^1 + (1-\Delta_x) \Delta_y f_x^2 + \Delta_x (1-\Delta_y) f_x^3 + \Delta_x \Delta_y f_x^4 .\end{split}\]

Using the first-order approximations for the derivatives at the node points (5) and (12)

(40)\[\begin{split} f^1_x & \approx & \frac{1}{2h} \left( f^{3} - f^{9} \right) \\ f^2_x & \approx & \frac{1}{2h} \left( f^{4} - f^{10} \right) \\ f^3_x & \approx & \frac{1}{2h} \left[ 3 f^{3} - 4 f^{1} + f^{9} \right] \\ f^4_x & \approx & \frac{1}{2h} \left[ 3 f^{4} - 4 f^{2} + f^{10} \right]\end{split}\]

in (39) gives

(41)\[\begin{split} f_x^0 & \approx & (1-\Delta_x)(1-\Delta_y) \frac{1}{2h} \left( f^{3} - f^{9} \right) + (1-\Delta_x) \Delta_y \frac{1}{2h} \left( f^{4} - f^{10} \right) \\ & & + \Delta_x (1-\Delta_y) \frac{1}{2h} \left[ 3 f^{3} - 4 f^{1} + f^{9} \right] + \Delta_x \Delta_y \frac{1}{2h} \left[ 3 f^{4} - 4 f^{2} + f^{10} \right] .\end{split}\]