A first-order computational algorithm for reaction-diffusion type equations via primal-dual hybrid gradient method

To clearly convey our proposed idea, let us first consider the following reaction-diffusion equation as an illustrative example on 2D periodic region $\Omega=\mathbb{T}^{2}$. We assume $\Omega$ is square shaped and denote its side length as $L$.

Here we assume $\lambda$ is a positive constant coefficient; $f:\mathbb{R}\rightarrow\mathbb{R}$ is the nonlinear function depicting the reaction term. Since we assume $\Omega$ to be the periodic region, we use periodic boundary conditions (BC) for equation (10).

Although there are numerous pieces of research on designing numerical schemes for RD equations, to demonstrate how our method works, let us narrow down and focus on the implicit one-step finite difference (FD) scheme. Once we have demonstrated how to implement the method to this implicit scheme, such a method can be easily extended to general numerical schemes.

We discretize the time interval $[0,T]$ into $N_{t}$ equal subintervals with length $h_{t}=T/N_{t}$. Suppose we discretize each side of the region $\Omega$ into $N_{x}$ subintervals with space stepsize $h_{x}=L/N_{x}$, we choose the central difference scheme to discretize the Laplace operator $\Delta$. We denote the discrete Laplace operator with the periodic boundary condition as $\textrm{Lap}_{h_{x}}^{\textrm{P}}$ which is an $N_{x}^{2}\times N_{x}^{2}$ block-circulant matrix possessing the following form

Here $I$ is the $N_{x}$ by $N_{x}$ identity matrix. We denote $U^{k}\in\mathbb{R}^{N_{x}^{2}}$ as the numerical solution of (10) at the $k$th time step. We vectorize the $N_{x}\times N_{x}$ square array along the column to form the 1D vector $U^{k}$. That is, for $l=i\cdot N_{x}+j$ with $1\leq i\leq N_{x}$ and $1\leq j\leq N_{x}$, $U^{k}_{l}$ is the numerical approximation of $u((j-1)h_{x},(i-1)h_{x},kh_{t})$.

Now the implicit one-step FD scheme for (10) is cast as

Here $U_{0}$ denotes the initial condition on mesh grid points. When solving for the numerical solution of (10), one has to sequentially solve a series of nonlinear equations as shown in (11). This is the place in which we should apply the PDHG method. Let us denote the function $F:\mathbb{R}^{N_{x}^{2}}\rightarrow\mathbb{R}^{N_{x}^{2}}$ as

We want to solve $F(U)=0$. As discussed in the introduction, this is equivalent to minimizing $\iota(F(U))$, which can further be cast as the following min-max saddle point problem

where $L$ is defined as $L(U,P)=P^{\top}F(U)$. As a result, solving the equation $F(U)=0$ finally boils down to the min-max problem (13).

Now the PDHG method suggests the following gradient ascent-descent dynamic for solving (13).

Similar to our discussion in the introduction, the third line above involves a nonlinear equation that cannot be directly solved. We thus can replace the term $L(U,\widetilde{P}_{n+1})$ in (16) by the linearization $\widehat{L}(U,P)=L(U_{n},P)+\nabla_{U}L(U_{n},P)(U-U_{n})$, then (16) can be explicitly computed as

Let us denote $U_{*}$ as the solution to $F(U)=0$. It is not hard to tell that $(U_{*},0)$ is a critical point of the functional $L(U,P)$. Furthermore, $(U_{*},0)$ is the equilibrium point of the time-discrete dynamic (14), (15), (17).

To analyze the convergence speed to the equilibrium state $(U_{*},0)$, we first consider the affine case in which $F(U)=AU-b$ with $A$ as an $N_{x}^{2}\times N_{x}^{2}$ symmetric matrix. We have the following result, similar to the analysis carried out in [31].

The proof of the theorem is provided in Appendix A. The explicit form of $\eta_{*}(\kappa)$ and $\gamma_{*}$ are given in remark 2 of Appendix A.

The optimal convergence rate $\gamma_{*}$ will be very close to $1$ if the condition number $\kappa$ is large. Furthermore, one will require $O(\kappa^{2})$ iterations for $U_{n}$ to converge. This can be very expensive when $A$ is a large-scale matrix with a large condition number. For example, we consider the heat equation

with periodic boundary conditions. We apply the one-step implicit scheme to this equation, i.e., we consider solving $(I-\lambda h_{t}\textrm{Lap}_{h_{x}}^{\textrm{P}})U^{n+1}=U^{n}$ at each time step $n$. Thus $A=I-h_{t}\lambda\textrm{Lap}_{h_{x}}^{\textrm{P}}$. One can tell that the eigenvalues of $A$ equal $1+4h_{t}\lambda N^{2}\sin^{2}\left(\frac{\pi k}{N}\right)$ for $1\leq k\leq N.$ When $N$ is even, the condition number $\kappa$ of $A$ equals $1+4\lambda N^{2}h_{t}$. Since we can get rid of the CFL condition by using the implicit scheme, we can pick $h_{t}\gg\frac{1}{N^{2}}$, which leads to $\kappa(A)\gg 1$. The convergence of the primitive PDHG method could be very slow, even for the heat equation.

As discussed in remark 2, $\gamma_{*}$ approaches $0$ if the condition number $\kappa$ drops to $1$. Hence, we need to control the condition number of the matrix $A$ to achieve faster convergence speed. This motivates us to introduce the preconditioning technique to the PDHG method. As suggested in both [25] and [31], we replace the $l^{2}$ norm used in either $\|U-U_{n}\|_{2}$ or $\|P-P_{n}\|_{2}$ by the $G-$norm $\|\cdot\|_{G}$ which is defined as

with $G$ as a symmetric, positive definite matrix. In this work, we mainly focus on substituting the norm $\|P-P_{n}\|_{2}$ with $\|P-P_{n}\|_{G}$. The PDHG method involving $G-$norm in its proximal step is sometimes named $G-$prox PDHG [25]. The dynamic obtained via such $G-$prox PDHG is shown below.

We should pause here to emphasize to the reader that the above three-line dynamic (18), (19), (20) will be the core gadget for our RD equation solver throughout the remaining part of the paper. Before we move on to further details on solving the RD equation (10) via $G-$prox PDHG dynamic, let us provide a little more explanation on how (18) - (20) improve the convergence speed $\gamma_{*}$. Analogous to Theorem 1, we have the following corollary for affine function $F$.

It is now clear that if we can find a matrix $G$ that approximates $AA^{\top}$ well, then $A^{\top}G^{-1}A$ will be reasonably close to the identity matrix $I$. Thus the condition number of $A^{\top}G^{-1}A$ will hopefully remain close to $1$. In such cases, by properly choosing the step size $\tau_{u},\tau_{p}$ such that $\tau_{u}\tau_{p}$ is close to $1$, we can obtain a rather fast convergence rate $\gamma_{*}$.

Up to this stage, although most of our intuition and analysis on $G-$prox PDHG dynamic comes from the case when $F$ is affine, it is natural to extend our treatment to the nonlinear $F(U)$ defined in (12). We may still anticipate the effectiveness of our method since (12) can be recast as

which can be treated as an affine function with a nonlinear perturbation $h_{t}f(U)$ carrying the small $h_{t}$ coefficient. We now discuss several details in applying the $G-$prox PDHG dynamic (18)-(20) to the above $F(U)$. We aim at evolving the following dynamic in order to update the implicit one-step scheme (11).

We summarize our method in the following algorithm. It is not hard to tell that the total complexity of each inner PDHG method is $O(\sharp\{\textrm{PDHG iter}\}\cdot N_{x}^{2}\log N_{x}\})$.

In this section, we mainly focus on the one-step implicit finite difference (FD) scheme to illustrate how we apply the PDHG iterations to update the given FD scheme. But we should emphasize that our method is not restricted to such a scheme. One can extend the PDHG method to various types of numerical schemes by formulating the scheme at a certain time step $k$ as a nonlinear equation $F^{k}(U)=0$, and construct the functional $L^{k}(U,P)=P^{\top}F^{k}(U)$. Then one can apply the dynamic (18) - (20) to update the numerical solution from $U^{k}$ to $U^{k+1}$. In addition, our method is applicable to more general reaction-diffusion equations (systems). Further discussions and details are supplied in section 3.

Theorem 1 suggests that under the linear case, the convergence of PDHG method relies on condition number $\kappa$, which is directly related to $h_{x},h_{t}$ of our discrete scheme. In practice, we fix $h_{x}$ and $\tau_{u},\tau_{p}$ in the algorithm. At every time step $n$ we discover that when time stepsize $h_{t}$ gets larger than a certain threshold value $h_{t}^{*}$ which depends on $h_{x},\tau_{u},\tau_{p}$ and $n$, PDHG method will hardly converge. The method works well when $h_{t}$ is slightly smaller than the threshold. The theoretical study on how $h_{t}^{*}$ guarantees the convergence of our method will be an important future research direction.

Adaptive time stepsize As discussed above, we cannot guarantee the convergence of the PDHG iteration (21) - (23) for any time stepsize $h_{t}$. Since the aforementioned threshold $h_{t}^{*}$ may vary at different time stages, and how $h_{t}^{*}$ varies depends on the nature of the equation as well as the discretization scheme. Given the potential difficulty of choosing the suitable $h_{t}$ that guarantees both the numerical accuracy as well as the convergence of the PDHG iterations, we come up with the strategy of using adaptive stepsize $h_{t}$ throughout the computation. We choose a time stepsize $h_{t}^{0}$ that guarantees the numerical accuracy and will serve as the upper bound of all $h_{t}$ throughout our method, we also set up two integers $N^{*}>N_{*}>0$ as the thresholding integers for enlarging or shrinking the time stepsize $h_{t}$. We also pick a rescaling coefficient $\eta\in(0,1)$. During each time step $n$ of the algorithm, we record the total number of PDHG iterations $M_{\mathrm{PDHG}}$, and reset the stepsize $h_{t}$ for the next time step based on the following rules:

• •

  If $M_{\mathrm{PDHG}}>N^{*}$, we shrink $h_{t}$ by rate $\eta$, $h_{t}=\eta h_{t}$;

• •

  If $N^{*}\geq M_{\mathrm{PDHG}}\geq N_{*}$, we remain $h_{t}$ unchanged;

• •

  If $M_{\mathrm{PDHG}}<N_{*}$, we enlarge $h_{t}$ by rate $\frac{1}{\eta}$, $h_{t}=\frac{h_{t}}{\eta}$, if $\frac{h_{t}}{\eta}\leq h_{t}^{0}$; we remain $h_{t}$ unchanged if $\frac{h_{t}}{\eta}>h_{t}^{0}$.

It is also reasonable to fix $h_{x},h_{t}$ but to only shrink the PDHG stepsizes $\tau_{u},\tau_{p}$ when we encounter difficulties in converging. However, according to our experience, shrinking $\tau_{u},\tau_{p}$ usually requires much more PDHG iterations for convergence, which may cause the algorithm less efficient compared with the aforementioned adaptive $h_{t}$ strategy.

The Allen-Cahn equation [1] is a typical reaction-diffusion equation taking the following form

Here $a,b>0$ are positive coefficients,

is a double-well potential function with its derivative $W^{\prime}(u)=u^{3}-u$. We will always assume periodic boundary conditions in our discussion.

The Allen-Cahn equation can be treated as the $L^{2}-$gradient flow of the following functional $\mathcal{E}(u)$.

We now switch to another well-known reaction-diffusion equation known as the Cahn-Hilliard equation [4], which takes the following form.

Here we assume $a,b>0$, and $W(u)$ defined the same as in the Allen-Cahn equation. In this section, we will restrict our discussion to periodic boundary conditions. Similar to the Allen-Cahn equation, the Cahn-Hilliard equation can be treated as the $H^{-1}-$gradient flow of the functional $\mathcal{E}(u)$ defined in (27). Due to this reason, compared with the Allen-Cahn equation, the Cahn-Hilliard equation involves one extra operator $-\Delta$ on the right-hand side of the equation. This difference leads to several slight modifications to our original algorithm.

The functional $F(U)$ introduced in (12) is now