# Fast stable finite difference schemes for nonlinear cross-diffusion

The dynamics of cross-diffusion models leads to a high computational complexity for implicit difference schemes, turning them unsuitable for tasks that require results in real-time. We propose the use of two operator splitting schemes for nonlinear cross-diffusion processes in order to lower the computational load, and establish their stability properties using discrete L^2 energy methods. Furthermore, by attaining a stable factorization of the system matrix as a forward-backward pass, corresponding to the Thomas algorithm for self-diffusion processes, we show that the use of implicit cross-diffusion can be competitive in terms of execution time, widening the range of viable cross-diffusion coefficients for on-the-fly applications.

## Authors

• 2 publications
• ### Stable Backward Diffusion Models that Minimise Convex Energies

Backward diffusion processes appear naturally in image enhancement and d...
03/08/2019 ∙ by Leif Bergerhoff, et al. ∙ 2

• ### An energy stable and positivity-preserving scheme for the Maxwell-Stefan diffusion system

We develop a new finite difference scheme for the Maxwell-Stefan diffusi...
05/16/2020 ∙ by Xiaokai Huo, et al. ∙ 0

• ### Fundamental Schemes for Efficient Unconditionally Stable Implicit Finite-Difference Time-Domain Methods

This paper presents the generalized formulations of fundamental schemes ...
11/28/2020 ∙ by Eng Leong Tan, et al. ∙ 0

• ### Theoretical analysis and simulation methods for Hawkes processes and their diffusion approximation

Oscillatory systems of interacting Hawkes processes with Erlang memory k...
03/24/2020 ∙ by Julien Chevallier, et al. ∙ 0

• ### Computation of the self-diffusion coefficient with low-rank tensor methods: application to the simulation of a cross-diffusion system

Cross-diffusion systems arise as hydrodynamic limits of lattice multi-sp...
11/22/2021 ∙ by Jad Dabaghi, et al. ∙ 0

• ### Optimal and Low-Memory Near-Optimal Preconditioning of Fully Implicit Runge-Kutta Schemes for Parabolic PDEs

Runge-Kutta (RK) schemes, especially Gauss-Legendre and some other fully...
12/23/2020 ∙ by Xiangmin Jiao, et al. ∙ 0

• ### Extending nonstandard finite difference schemes rules to systems of nonlinear ODEs with constant coefficients

In this paper, we present a reformulation of Mickens' rules for nonstand...
07/12/2021 ∙ by Marc Songolo, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Cross-diffusion models consist of evolutionary systems of diffusion type for at least two real-valued functions, where the evolution of each function is not independent of the others. Their use is widespread in areas like population dynamics (see [1] and references therein), but has recently attracted some interest in the image processing community [2, 3, 5, 7] as a natural extension to the complex diffusion proposed by Gilboa et al. [15], where the image is represented by a complex function and the filtering process is governed by a nonlinear PDE of diffusion type with a complex-valued diffusion coefficient. This equation can be written as a cross-diffusion system for the real and imaginary parts of the image, and enhanced imaging possibilities emerge via the extra nonlinear diffusion coefficients [5, 7].

It is known that implicit schemes have better stability properties when compared to their explicit counterparts [9], allowing larger time steps and, in the particular case of cross-diffusion, a wider range of diffusion coefficients [5, 7]. However, the computational complexity of traditional implicit finite difference schemes turns their use impracticable for tasks that aim for on-the-fly results, such as image processing in medical contexts. In this field, normally used strategies to speed-up computations include operator splitting techniques, multigrid schemes or, more recently, fast explicit diffusion [23, 24]. Several splitting techniques have been proposed since the middle of the previous century. First introduced by Douglas, Peaceman and Rachford [13, 12], they were further explored by soviet mathematicians such as Marchuk or Yanenko [25, 20]. Overall, their aim is the reduction of a multidimensional problem to a sequence of one-dimensional problems, in order to obtain simplified linear systems that can be efficiently processed by computers. In the case of the heat equation, the factorization of a point stencil spatial operator leads to a series of tridiagonal systems that can be processed by Thomas algorithm [21, 17].

In this report we propose the adjustment of two splitting techniques commonly used in the imaging community, usually referred to as Additive Operator Splitting (AOS) and Additive Multiplicative Operator Splitting (AMOS), to the nonlinear cross-diffusion case. AOS schemes are very popular since their introduction to the PDE-related imaging field by Weickert [23], but their source dates back a few years to the works of Tai and Neitaanmaki [22]. These authors noted that the classical splitting-up methods can not be used for parallel processors as the computation of the current fractional step always required the knowledge of the previous fractional step. They proposed a new splitting where the computation of the fractional steps are independent of each other, and therefore parallelizable. Further work, now also with Lu [19]

, explored convergence estimates for steady and nonsteady states nonlinear problems and for linear and quasilinear evolution equations. AMOS schemes belong to the family of multiplicative locally one-dimensional schemes, the ones where each fractional step computation requires the value of the previous fractional step. They rely on the same factorization of AOS, but each operator is now applied on top of each other (being called multiplicative). As, in general, the split operators do not commute, the final result depends on the order of the one-dimensional operators. This is a major drawback in some applications, particularly in image processing, as one usually aims for rotation invariant filters. In order to overcome this problem Barash and Kimmel

[6] added a symmetric setting by applying a multiplicative scheme in all possible orders and averaging the results.

The immediate application of these techniques to our problem still has some drawbacks, as the numerical system to solve remains dense due to the interdependence of the evolution functions inherent to the cross-diffusion model. To overcome this issue we reorder the unknowns in order to obtain a tridiagonal by blocks system matrix, allowing for a block LU factorization solver to act as a forward-backward pass, similar to the process given by Thomas algorithm [17].

The text is organized in the following way. In Section 2 we describe the main model and the standard finite difference -method for two-dimensional problems, and we derive its stability properties. In Section 3 we introduce the AOS and AMOS implementations in cross-diffusion, and we extend the analogous stability results and discuss their computational implementations. Still in this section, we describe how to design the system matrix as tridiagonal by blocks, and discuss the stability conditions for such factorization. In Section 4 we extend the results to the three-dimensional case. Finally, in Section 5 we present the theoretical speedups and some numerical experiments.

## 2 Model and standard scheme

A nonlinear cross-diffusion with semi-linear reaction process can be represented by a two-component vector field,

, satisfying the system

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ut=∇⋅(d1(u,v,t)∇u+d2(u,v,t)∇v)−λ1(u,v,t)(u−u0) in Ω×R+,vt=∇⋅(d3(u,v,t)∇u+d4(u,v,t)∇v)−λ2(u,v,t)(v−v0) in Ω×R+,u(x,0)=u0(x), v(x,0)=v0(x) in Ω,uη=0, vη=0 on Γ×R+, (1)

where is the domain of interest, and are the given initial conditions for and and denotes the outward normal vector to the boundary . The cross-diffusion matrix of the model is given by

 D(u,v,t)=[d1(u,v,t)d2(u,v,t)d3(u,v,t)d4(u,v,t)], (2)

with , the so called influence functions, and non-negative real valued bounded functions. Before moving on to the theoretical results for finite difference implementations of (1), we start by grounding ourselves with some definitions. More precisely, we state our understanding of positive definite matrices:

###### Definition 1.

Let be a square matrix of size . We say that is positive (semi-)definite if

 x⊤Ax>0(≥0),∀x∈Cn.

We will make an intensive use of a similar class of matrices, where the implicit necessity for hermiticity in positive definiteness is no longer required:

###### Definition 2.

Let be a real valued square matrix of size . We say that is positive (semi-)definite, not necessarily symmetric, if

 x⊤Ax>0(≥0),∀x∈Rn.

We also unambiguously define the notion of stability finite difference schemes following the line of Jovanović and Suli [9].

###### Definition 3.

Let be the -step of a finite difference scheme with space-mesh parameter and time step applied to problem (1). Let be a mesh-dependent norm involving mesh points of . We say that the scheme is unconditionally stable if there is independent of the mesh size such that

 ∥(¯Um,¯Vm)∥ϕh≤C∥(U0,V0)∥ϕh

holds for any choice of parameters and , where are the initializations and taken at the mesh points. If it holds only for , where is a function of , we say that the scheme is conditionally stable.

#### Explicit and Semi-implicit implementations and stability results

Let the domain be discretized by the points , where

 xj1=a1+h1j1,xj2=a2+h2j2,jk=0,1,…,Nk,
 hk=bk−akNk,k=1,2,

for two given integers , and . This spatial mesh on is denoted by and . Points halfway between two adjacent grid points are denoted by , , where is the canonical basis, that is, is the standard basis unit vector in the th direction.

For the discretization in time, we consider a mesh with time step , , where .

We denote by the value of a mesh function at the point . For the formulation of the finite difference approximations, we use the centered finite difference quotients in the th spatial direction, for ,

 δkZj=Zj+(1/2)ek−Zj−(1/2)ekhk,
 δkZj+(1/2)ek=Zj+ek−Zjhk.

An initial distribution is required, given by two real-valued functions . Let , such that . Given the initial solution , the numerical solution of (1) at the time is obtained considering the following finite difference scheme

The numerical solution of (1) at the time can be obtained considering the following finite difference scheme:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Um+1j−UmjΔt=2∑k=1δk(d1(Wm)m+θjδkUm+θj+d2(Wm)m+θjδkVm+θj)−λ1(Wm)m+θj(Um+θj−U0j),Vm+1j−VmjΔt=2∑k=1δk(d3(Wm)m+θjδkUm+θj+d4(Wm)m+θjδkVm+θj)−λ2(Wm)m+θj(Vm+θj−V0j), (3)

where

 dℓ(Wm)m+θj±(1/2)ek=dℓ(Wmj,tm+θ)+dℓ(Wmj±ek,tm+θ)2,ℓ=1,…,4, (4) λi(Wm)m+θj=λi(Wmj,tm+θ),i=1,2, Zm+θj=θZm+1j+(1−θ)Zmj,Z=U,V,

and corresponding to, respectively at the extremes, explicit and semi-implicit implementations, and to a semi-Crank-Nicholson implementation when .

For each , we define the rectangle and denote by the measure of . We consider the discrete inner products

 (U,V)h =∑□j⊂Ω|□j|4(Uj1,j2Vj1,j2+Uj1+1,j2Vj1+1,j2 (5) +Uj1,j2+1Vj1,j2+1+Uj1+1,j2+1Vj1+1,j2+1),
 (U,V)h∗1=∑□j⊂Ω|□j|2(Uj1+1/2,j2Vj1+1/2,j2+Uj1+1/2,j2+1Vj1+1/2,j2+1),

and

 (U,V)h∗2=∑□j⊂Ω|□j|2(Uj1,j2+1/2Vj1,j2+1/2+Uj1+1,j2+1/2Vj1+1,j2+1/2).

Their correspondent norms are denoted by , and , respectively. For we define .

To simplify the notation and when it is clear from the context, we will write instead of or , for .

Theorem 1 states the stability conditions for this finite difference scheme, which is an extension of a previous result by the authors presented in [7]. Define as the supreme of both and .

###### Theorem 1.

If the cross-diffusion matrix (2) is such that, for any ,

 2∑k=1(d1δkU,δkU)h∗k+(d2δkV,δkU)h∗k+(d3δkU,δkV)h∗k+(d4δkV,δkV)h∗k≥0 (6)

then the scheme (3) is unconditionally stable for .

If the cross-diffusion matrix (2) is such that, for any ,

 2∑k=1 (d1δkU,δkU)h∗k+(d2δkV,δkU)h∗k+(d3δkU,δkV)h∗k+(d4δkV,δkV)h∗k (7) −(4Δth2k((1+η1k)(∥d1δkU∥2h∗k+∥d2δkV∥2h∗k+2(d1δkU,d2δkV)h∗k)

for some , , then the scheme (3) is unconditionally stable for the case .

###### Proof.

Let be any iteration of the finite difference scheme (3). Multiply both members by and , respectively, according to the discrete inner product , use summation by parts, the non-negativity of the reaction multiplier, and the discrete Duhamel principle [11] to obtain, for the case ,

 ∥Wm+1∥2h≤e2(θ2+(1−θ)2)~ϵ−1ϵtm+1(1+tm+1λmaxϵ−1~ϵ−1)∥W0∥2h,

with , for some . For the case , take the same steps and square (3) to obtain

 ∥Wm+1∥2h≤eaϵtm+1(1+tm+1bϵ)∥W0∥2h, (8)

with and , where . ∎

We supplement this result with a series of corollaries related to common instances of cross-diffusion matrices [2, 3, 4, 8, 15].

###### Corollary 1.

If the cross-diffusion matrix (2) is semi-positive definite, not necessarily symmetric, then the method is unconditionally stable for any .

###### Proof.

Notice that we can write (6) as

 2∑k=1∑□j⊂Ω∑i=0,1|□j|2[δkUj+ek/2+ielδkVj+ek/2+iel]⊤[d1(W)j+ek/2+ield2(W)j+ek/2+ield3(W)j+ek/2+ield4(W)j+ek/2+iel][δkUj+ek/2+ielδkVj+ek/2+iel],

where is the direction orthogonal to . If the cross-diffusion matrix is semi-positive definite, taking (4) into account gives that the summand is non negative. ∎

###### Corollary 2.

If the cross-diffusion matrix (2) is of the form

 [d1(W)d2(W)d3(W)d4(W)]=g(W)M,

with a non-negative function and a semi-positive definite matrix, not necessarily symmetric, then the method is unconditionally stable for any .

###### Proof.

Follows immediately from the previous Corollary. ∎

###### Corollary 3.

If the functions and satisfy

 d1(x,y,t)≥12|d2(x,y,t)+d3(x,y,t)|,∀(x,y,t)∈R2×R+, (9) d4(x,y,t)≥12|d2(x,y,t)+d3(x,y,t)|,∀(x,y,t)∈R2×R+,

then the scheme (3) is unconditionally stable for any .

###### Proof.

Follows from the first corollary by realizing that a matrix is positive semi-definite in if and only if its symmetric part is positive semi-definite, and that conditions (9) force the symmetric part to be diagonally dominant by rows, a sufficient condition for positive semi-definiteness [21]. ∎

## 3 Operator splittings

The stability conditions suggest that there is a class of functions for which explicit implementations of cross-diffusion processes require unpractical small time steps. However, implicit methods are numerically expensive to solve due to the dimension of the spatial variable. As mentioned before, we will focus on Additive Operator Schemes (AOS) and Additive Multiplicative Operator Schemes(AMOS).

### 3.1 Splitting models and L2 stability

Applying the AOS technique [23] to our finite difference scheme (3) translates into splitting the spatial operators in each direction, calculate each directional time step and then average the resulting fractional steps. We shall henceforth refer to it as the AOS-CD scheme:

 (10a) for k=1,2, and Um+1j=122∑k=1Um+1,kj,Vm+1j=122∑k=1Vm+1,kj. (10b)

Before moving on to the numerical considerations for this new AOS-CD scheme, we introduce the other splitting technique that we will use to fasten the implicit part of (3). The motivation is the fact that despite their efficiency and stability, AOS schemes have limited accuracy. An application of AMOS technique [6] to our problem leads to the following AMOS-CD schemes

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Um+1,∗j−UmjΔt=δ1(d1(Wm)m+θjδ1Um+θj+d2(Wm)m+θjδ1Vm+θj)−12λ1(Wm)m+θj(Um+θj−U0j)Vm+1,∗j−VmjΔt=δ1(d3(Wm)m+θjδ1Um+θj+d4(Wm)m+θjδ1Vm+θj)−12λ2(Wm)m+θj(Vm+θj−V0j)Um+1,⋆j−UmjΔt=δ2(d1(Wm)m+θjδ2Um+θj+d2(Wm)m+θjδ2Vm+θj)−12λ1(Wm)m+θj(Um+θj−U0j)Vm+1,⋆j−VmjΔt=δ2(d3(Wm)m+θjδ2Um+θj+d4(Wm)m+θjδ2Vm+θj)−12λ2(Wm)m+θj(Vm+θj−V0j) (11a) for the first step, ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Um+1,∗∗j−Um+1,∗jΔt=δ2(d1(Wm)m+θjδ2Um+θj+d2(Wm)m+θjδ2Vm+θj)−12λ1(Wm)m+θj(Um+θ,∗j−U0j)Vm+1,∗∗j−Vm+1,∗jΔt=δ2(d3(Wm)m+θjδ2Um+θj+d4(Wm)m+θjδ2Vm+θj)−12λ2(Wm)m+θj(Vm+θ,∗j−V0j)Um+1,⋆⋆j−Um+1,⋆jΔt=δ1(d1(Wm)m+θjδ1Um+θj+d2(Wm)m+θjδ1Vm+θj)−12λ1(Wm)m+θj(Um+θ,⋆j−U0j)Vm+1,⋆⋆j−Vm+1,⋆jΔt=δ1(d3(Wm)m+θjδ1Um+θj+d4(Wm)m+θjδ1Vm+θj)−12λ2(Wm)m+θj(Vm+θ,⋆j−V0j) (11b) for the second step, and finally Um+1j=12(Um+1,∗∗j+Um+1,⋆⋆j),Vm+1j=12(Vm+1,∗∗j+Vm+1,⋆⋆j). (11c)

Just a quick remark regarding the factorization steps (11a) and (11b). In each bracket there are two independent systems of difference equations (in fact there are several of them, as we will discuss in the next section). We chose to display them in this way in order to clarify that systems (11b) require the solutions of systems (11a).

In the next result we investigate the numerical stability of both schemes:

###### Theorem 2.

If the cross-diffusion matrix (2) is such that, for any ,

 (d1δ1U,δ1U)h∗1+(d2δ1V,δ1U)h∗1+(d3δ1U,δ1V)h∗1+(d4δ1V,δ1V)h∗1≥0, (12)

then the schemes (10) and (11) are unconditionally stable for .

If the cross-diffusion matrix (2) is such that, for any ,

 (d1δ1U,δ1U)h∗1+(d2δ1V,δ1U)h∗1+(d3δ1U,δ1V)h∗1+(d4δ1V,δ1V)h∗1 (13) −(4Δth21((1+η1)(∥d1δ1U∥2h∗1+∥d2δ1V∥2h∗1+2(d1δ1U,d2δ1V)h∗1)

for some , and for , then the schemes (10) and (11) are unconditionally stable for the case

###### Proof.

We start with the AOS-CD scheme, and proceed as in Theorem 1 for each direction, taking into account that the action in direction on and is the same as the action in direction on and and the triangle inequality to obtain, for ,

 ∥Wm+1∥2h≤e4(θ2+(1−θ)2)~ϵ−1ϵtm+1(1+2¯ϵ−1tm+1λmax)∥W0∥2h,

with for some , while for the case the same steps as in Theorem 1 yield

 ∥Wm+1∥2h≤eaϵtm+1(1+tm+1bϵ)∥W0∥2h, (14)

with and , where and . For the AMOS-CD scheme, follow the same steps to get similar inequalities for each factorization step (11a) and (11b) and finish with the triangle inequality in (11c). For the case we obtain

 ∥Wm+1∥2h≤e8(θ2+(1−θ)2)~ϵ−1ϵtm+1(1+tm+1λmaxϵ−1~ϵ−1)∥W0∥2h,

with for some , while for the case we obtain

 ∥Wm+1∥2h≤eaϵtm+1(1+tm+1bϵ)∥W0∥2h, (15)

with and

### 3.2 Computational considerations

We now discuss the computational implementations for each scheme. The matrix formulation for the -method (3) can be obtained by arranging the elements of in a vector defining

 w0=[U1,1U2,1U3,1⋯UN1,N2V1,1V2,1V3,1⋯VN1,N2]⊤,

in what we call an ordering of pixels prioritizing the first direction (order first in direction and second in direction ). As such, setting

 As11=2∑ϑk=1ϑlkD1k(wm)sϑrk,As12=2∑ϑk=1ϑlkD2k(wm)sϑrk, (16) As22=2∑ϑk=1ϑlkD4k(wm)sϑrk,As21=2∑ϑk=1ϑlkD3k(wm)sϑrk,

for , where denote the backward difference operators with respect to direction , denote the forward difference operators with respect to direction , and is a diagonal matrix with , for , following the backward average (4) in the corresponding direction , and

 As:=[As11As12As21As22], (17)

and

 Λi(wm)s=diag{λi(wm)s1,1,λi(wm)s2,1,…,λi(wm)sN1,N2,N3},i=1,2,s=m,m+1, Λs=diag{Λ1(vm)s,λ2(vm)s},s=m,m+1,

the -method iteration in (3) reads

 (I+(1−θ)ΔtAm+(1−θ)θΔtΛm)wm (18) +Δt(θΛm+1+(1−θ)Λm)w0,

for . On the other hand, for each of the factorization steps (10a), (11a) and (11b), the iteration reads as

 (I−θrAm+1k+θrΛm+1k)wm+1k=(I+(1−θ)rAmk−(1−θ)rΛmk)wmk+r(θΛm+1k+(1−θ)Λmk)w0k, (19)

where is vector with an orientation prioritizing direction , follows the same orientation, for the AOS-CD splitting scheme and for the AMOS-CD scheme, and matrix is defined as