A kernel-free boundary integral method for reaction-diffusion equations

1.   Introduction

Let $\Omega \subset \mathbb{R} ^{2}$ be a bounded domain. Define the reaction-diffusion equation

subject to both the initial condition

and the Dirichlet boundary condition

where $u({\mathbf p}, t)$ is a function that describes the concentration of substances concerning time $t$ and the two-dimensional space vector ${\mathbf p}$, $D$ is the diffusion coefficient, which is used to express the rate of diffusion of a substance through space, $R(u)$ describes a chemical or biological reaction as a function of concentration $u$, $j({\mathbf p})$ is the initial data in $\Omega$ at time $t = 0$, and $g^{D}({\mathbf p}, t)$ is the boundary data on $\partial \Omega$. The first term on the right-hand side of Eq (1.1) represents the diffusion of a substance through space, and the second term represents the concentration-dependent reaction. The two parts of Eq $(1.1)$ can lead to a singular perturbation problem. One is when $D$ is small and the diffusion process is completed in a very short time, and the other is when the reaction term $R(u)$ is nonlinear, both of which can lead to dramatic oscillations in space. Reaction-diffusion equations for a singular perturbation have been addressed in several fields, including physics^[1,2,3], ecology^[4,5,6], chemistry^[7], and climate science^[8].

For singular perturbation problems, the use of conventional finite-difference methods on uniform grids can lead to significant oscillations near the boundaries. An effective approach is to improve the mesh by enabling automatic adjustments of the mesh resolution based on the properties and requirements of the problem. Gupta et al. ^[9] generated a spatially adaptive mesh using an equidistributed positive monitoring function to automatically detect the position, height, and thickness of the parabolic boundary layer in the solution. Zhang and Liu ^[] conducted convergence analysis on rectangular grids of Bakhvalov type, a special finite element mesh structure. Kaushik et al. ^[11] recursively generated modified hierarchical grids by constructing a mesh distribution function, which served as an upper bound on the optimization error. The Shishkin grid comprises internal and external layers with the former being a denser grid. A transition layer ^[12] exists between the two layers to ensure a smooth transition. Besides, some variations like segmented uniform Shishkin grids ^[13] have been proposed to get more flexibility.

Another effective approach is to improve the numerical methods tailored to the specific problem. The finite difference method (FDM) and the finite element method (FEM) are the most commonly used numerical techniques for solving the singular perturbation problem. In the framework of the FDM, non-isometric finite difference schemes ^[14], hybrid finite difference schemes ^[13], and special FDMs incorporating Euler's method ^[15] have been successively proposed. High-order FEMs ^[11], bilinear FEMs ^[16], FEMs in h-p version ^[17,18], and other novel FEMs have been extensively used. The reaction-diffusion problem can be solved in time and space separately by applying the multiscale methods, where different FDMs^[13,15] are applied to time and space, respectively, based on the problem's characteristics. In addition to this, Galerkin's method ^[19,20,21] as an application of the variational method is often applied to solve singular perturbation problems. When the problem has significant differences in characteristic scales, the numerical asymptotic method ^[22] allows the problem to be divided into different regions, where different numerical methods are used separately.

The high-order compact FDM ^[23,24,25], when combined with multigrid methods ^[26], provides an efficient approach for discretizing equations on both uniform and non-uniform grids. For uniform grids, Zhang ^[27] and Gupta et al. ^[28,29] proposed multigrid methods and high-order compact schemes, which were later extended to sixth-order accuracy by Spotz and Carey ^[30], and Kyei et al. ^[31]. On the other hand, for non-uniform grids, researchers such as Ge ^[32] and Ge et al. ^[33] addressed irregular grid challenges through transformation-free high-order schemes. Similarly, Mayo and Greenbaum ^[34] proposed high-order methods for external integral problems, effectively handling complex grid geometries. Additionally, Saldanha ^[35] and Boisvert ^[36] proposed general-purpose methods that bridge the gap between uniform and non-uniform grids, demonstrating flexibility across various computational scenarios. Hackbusch and Nowak ^[37] proposed an approximate matrix-vector multiplication method for solving linear systems with full matrices in the boundary element method (BEM) and sparse matrices in FEM, which saves a large number of computations and has been subsequently applied. To further improve computational speed and memory efficiency, fast solvers have been proposed and applied to solve the equation on a mesh ^{[38,39,40,41,42]}, which contains a high-order boundary integral equation solver^[43], including the "black box" solver ^[44] and FFT-based solver^{[45,46,47,48,49]}.

As we all know, the traditional boundary integral method (BIM) ^[50] requires the analytical expression of the kernel function for integral evaluation. To circumvent the calculation of complex Green's functions, Ying and Henriquez ^[46] proposed the KFBI method, which transforms elliptic problems into boundary integral equations (BIEs) and approximates the boundary or volume integrals by the FDM. This can be naturally achieved since those integrals are exactly the solutions to some equivalent simple interface problems according to the potential theory^[51]. In this way, analytical expressions of Green's functions are no longer needed to be determined, formulated, or computed in the whole calculation process. This is the true essence of the "kernel-free" method. Besides, the coefficient matrix of the discrete system always remains unchanged for equivalent interface problems, thus the system can be solved with a fast elliptic solver. The BIEs for second-order elliptic problems are the Fredholm integral equations of the second kind and can be solved efficiently by any Krylov subspace iterative method since the condition number of the resulting system is essentially independent of the mesh parameter. These advantages constitute the KFBI method to be an efficient, accurate, and stable solver for elliptic problems in general domains. In this paper, we employ the KFBI method as the spatial discretization method and use the GMRES method ^[45,52] to iteratively solve the corresponding BIEs.

The rest of the paper is organized as follows. Section 2 discretizes the reaction-diffusion equation using the two-step semi-implicit Runge-Kutta method reducing it to the form of the modified Helmholtz equation in each time step. Section 3 introduces the boundary value problem (BVP) and its integral form, reformulating the Dirichlet BVP into the corresponding BIE. Section 4 addresses simple interface problems equivalent to boundary or volume integrals involved in the BIE. Section 5 presents the KFBI method for solving the interface problem. Specific numerical examples are provided in Section 6.

2.   Time integral method

The two-step semi-implicit Runge-Kutta (SIRK) method is selected to discretize the time variable $t$ in the reaction-diffusion equation for the singular perturbation problem. The two-step SIRK method is a numerical integration algorithm that combines the stability of the implicit Runge-Kutta method with the geometric properties of the symplectic method. Consider $\partial u / \partial t = f(u) + g(u)$ where $f(u)$ is a rigid term. To balance the computational effort and stability, both implicit and explicit processing of rigid and non-rigid terms are employed to derive

Here $u^{n} \approx u({\mathbf p}, t^{n})$ represents the approximation at the time $t^{n}$, where $\Delta t = t^{n+1} - t^{n}$ is the time step. The two-stage SIRK method divides time step $\Delta t$ into two stages. The interval is considered $(t^{n}, t^{n}+\frac{1}{2} \Delta t)$ as the first stage, and $(t^{n} + \frac{1}{2} \Delta t, t^{n} + \Delta t)$ for the second stage. $K_{0}$ and $K_{1}$ are intermediate variables for the slope of the first and second stages, respectively. By combining $K_{0}$ and $K_{1}$, a second-order time-integral approximation is constructed. In the SIRK scheme (2.1), the first two steps are semi-implicit, while the third step is fully explicit. This structure allows the method to preserve the symplecticity and compatibility of the system while ensuring numerical stability. Define $v_{0} = u^{n}$, $v_{1} = u^{n} + \frac{1}{2}K_{0}\Delta t \approx u ^{n + \frac{1}{2}}$, and $v_{2} = u^{n} + \frac{1}{2}K_{1}\Delta t \approx u ^{n + \frac {1}{2}}$. Then the simplification of (2.1) reads

Substituting $f(v) = D\Delta v$ and $g(v) = R(v)$ into $(2.2)$ will reduce the first two equations to the form of the modified Helmholtz equation such that

According to the Dirichlet boundary condition (1.3), we have

This gives us two elliptic partial differential boundary value problems:

and

Solutions $v_{1}$ and $v_{2}$ for (2.5) and (2.6) are solved with the KFBI method, and $u^{n+1}$ is updated by substituting into $u^{n+1} = v_{1} + v_{2} - v_{0}$. The specific calculation process for the elliptic boundary value problem with the KFBI method is described in the following sections.

3.   Boundary value problem

Suppose $\Omega \subset \mathbb{R} ^{2}$ is a closed domain and $\Gamma = \partial \Omega$ represents its smooth boundary. We consider the following Dirichlet boundary value problem (BVP):

where $\mathcal{A}$ is a differential operator for the modified Helmholtz equation, $u$ is an unknown function, $\Delta$ denotes the Laplace operator, ${\mathbf p}$ is a spatial variable in $\Omega$, $\kappa$ is a constant, $f$ is a smooth function defined in $\Omega$, and $g^{D}$ represents the Dirichlet boundary data.

Let $\mathcal{B} \subset \mathbb{R}^{2}$ be a large rectangle containing $\Omega$, and let $G$ denote Green's function to the differential operator $\mathcal{A}$ on $\mathcal{B}$. Green's function satisfies the following conditions:

for any ${\mathbf q} \in \mathcal{B}$. Here $\delta ({\mathbf p}-{\mathbf q})$ denotes the Dirac function, and $\mathcal{A}$ is the differential operator defined earlier. According to potential theory^[51], the Dirichlet BVP can be reformulated as:

The solution to the Dirichlet BVP is

where $\varphi({\mathbf q})$ is the density function and ${\mathbf n_{q}}$ is the unit normal vector outward at an interface point ${\mathbf q} \in \Gamma$. Equation $(3.4)$ can be solved either by Krylov subspace methods such as the GMRES method ^[52] or the simple Richardson iteration, which is given by

where ${\mathbf p} \in \Gamma$ is the spatial vector on the boundary, and $\gamma$ is the iteration parameter with values in the range (0, 2). As the iteration converges, $u({\mathbf p})$ can be calculated by formula (3.5).

4.   KFBI method

The KFBI method is usually used to solve problems with smooth irregular surfaces, where the elliptic PDEs are transformed into BIEs using Green's first or second formula. The integrals in these equations are computed by a Cartesian grid-based finite difference method. It does not require the explicit expression of the kernel function, thereby avoiding the direct calculation of complex volume and boundary integrals.

4.1.   Redefinition of integrals

Let $\mathcal{B}\subset \mathbb{R}^{2}$ be a rectangular region in the plane. $\Gamma \subset \mathcal{B}$ is a smooth interface that divides $\mathcal{B}$ into two regions, $\Omega _{i}$ and $\Omega _{e}$, i.e., $\Gamma = \bar{\Omega}_{i}\cap \bar{\Omega }_{e}$ and $\mathcal{B} = \Omega_{i}\cup \Omega_{e}\cup \Gamma$. Let ${\mathbf n}$ denote the outward unit normal vector on $\Gamma$. With Green's function defined in $\mathcal{B}$, the double layer potential is presented by

The volume integral is defined as

By the continuity/discontinuity of the boundary and volume integrals, it is possible to reformulate the integrals as their equivalent simpler interface problems. The double-layer potential $v({\mathbf p}) = \mathcal{M}\varphi({\mathbf p})$ satisfies the following interface problem:

The discontinuity of the double-layer potential is given by

The volume integral $v({\mathbf p}) = \mathcal{G}f({\mathbf p})$ satisfies

Note that the interface conditions above reflect the continuous property of the potential function $v$ and its normal derivative $\partial_{{\mathbf n}} v$.

4.2.   Interface problem

For a general function $u$ defined in $\mathcal{B}$, which is continuous only in $\Omega_{i}$ and $\Omega_{e}$, and discontinuous across the interface $\Gamma$. We define the unknown restriction functions $u_{i} = u|_{{\mathbf p} \in \Omega_i}$ and $u_{e} = u|_{{\mathbf p} \in \Omega_e}$, while $f_{i} = f_{i}({\mathbf p})$ and $f_{e} = f_{e}({\mathbf p})$ are, respectively, the smooth known functions defined in $\Omega_{i}$ and $\Omega_{e}$. We now consider the interface problem

subject to the interface conditions

and the Dirichlet boundary condition

To explain the interface conditions (4.7), let $u^{+}({\mathbf p})$ be the one-sided limit of $u$ at an interface point ${\mathbf p}$ from the subdomain $\Omega_{i}$, and let $u^{-}({\mathbf p})$ be the one-sided limit at ${\mathbf p}$ from $\Omega_{e}$. The difference between the two is then defined as the jump function $\varphi$ on the interface. Similarly, the difference between the limits of their normal partial derivatives is defined as the jump in the partial derivatives $\psi$:

With the notations and properties of the integrals given above, BIE (3.4) can be concisely expressed as

where $b = g^{D}-\left(\mathcal{G}f \right) ^+$ is the right-hand side. The discretization points on the interface $\Gamma$ can be simply selected as quasi-uniformly spaced points through uniformly distributed curve parameters, or intersection points of the interface with Cartesian grid lines if the interface is implicitly defined. Then the unknown densities are discretized on these points, thus the dimension of the discrete system is exactly the number of discrete points on the interface.

Two kinds of iteration methods are feasible to solve the BIE. The Richardson iteration (3.6) mentioned before can be briefly rewritten in the form of

To efficiently solve the discrete system, we consider the generalized minimal residual (GMRES) method ^[45], which is a Krylov subspace method for solving general linear system $Ax = b$. This method iteratively constructs an orthonormal basis of the Krylov subspace $\mathcal{K}_{m} = \mathrm{span}\{r_0, Ar_0, A^{2}r_0, \dots, A^{m-1}r_0\}$, where $r_0 = b - Ax_0$ is the initial residual. At the $k$-th step, the algorithm minimizes the residual norm $\left \| b - Ax_{k} \right \|$ over the subspace $\mathcal{K}_m$. The key steps in the GMRES method are described in Algorithm 1.

In this work, the coefficient matrix $A$ in the linear system is not explicitly given but is defined by the integral operator. To use the GMRES method, BIE (4.10) is regarded as a linear system, where the discrete vector $\varphi$ is unknown to be solved, the discrete boundary integral operator $\left(\mathcal{M}\varphi\right)^+$ is considered as the matrix-vector product $"Ax"$, and $b = g^{D}-\left(\mathcal{G}f \right) ^+$ is the right-hand side computed in advance. This information is substituted into the GMRES procedure. In the iterative process, when a matrix-vector product is encountered, we compute the corresponding integral as its counterpart. The GMRES method is particularly well-suited for large, sparse systems. Specifications for integral computation are presented in the following subsections.

4.3.   Discretization with the compact difference method

Let $\mathcal{B}$ be a square box with sides $(a, b)$ and $(c, d)$ that contains $\Gamma$, $M$ is the number of partitions on each side, and $h = (b-a)/M = (d-c)/M$ is the grid parameter. The grid nodes can then be represented by $(x_{i}, y_{j})$, where $x_{i} = a + ih, y_{j} = c + jh$ with $i, j = 0, 1, 2, \dots, M$. For a general function $v$ defined on $\mathcal{B}$, denote $v_{i, j} \approx v(x_{i}, y_{j})$ as its grid approximation. For a given function $f$, denote $f_{i, j} = f(x_{i}, y _{j})$ as its value on the grid. Based on the traditional five-point finite difference scheme for second-order partial differential approximation^[53], one obtains the following five-point finite difference scheme for the modified Helmholtz equation:

Here $\mathcal{A}_h$ represents the difference operator for the partial differential equations and $\kappa$ is a constant.

In the discretization process, the standard finite difference scheme produces large local truncation errors near the interface. Thus, it is necessary to identify the irregular grid points near the interface before correcting the truncation error. For a given grid point $P$, the four neighboring grid points $P_1, P_2, P_3,$ and $P_4$ along with the centering point $P$ constitute the five-point compact difference stencil, as shown in Figure 1. Grid points whose five-point stencil intersects the interface are classified as irregular points. More concisely, two adjacent grid points located on different sides of the interface are considered to be irregular. For simplicity, we consider only the case where the five-point stencil intersects the interface once.

Suppose $(x_i, y_j)$ is an irregular grid point in the subdomain $\Omega_i$, whose right-hand neighboring grid $(x_{i+1}, y_j)$ lies in $\Omega_e$. The second-order finite difference scheme at $(x_i, y_j)$ is given by

where the superscript "$+$" is used because the centering point $(x_i, y_j)$ is an interior point $\Omega_i$. However, since $(x_{i+1}, y_j)$ lies in $\Omega_e$ as an assumption, this leads to a large local truncation error. Next, we will estimate the leading order term of the local truncation error using a Taylor expansion. Denote $z_{i}$ as the intersection point of the five-point stencil and the interface in the horizontal direction, with $x_{i}\le z_{i} < x_{i+1}$. We perform a Taylor expansion of $v^+(x, y_{j})$ and $v^-(x, y_{j})$ at the intersection point $(z_{i}, y_{j})$, and then

Substituting the Taylor expansion (4.14) into the calculation of the truncation error yields

where $[v] = v^{+}(z_{i}$, $y_{j}) - v^{-}(z_{i}, y_{j})$, $[v_{x}] = v_{x}^{+}(z_{i}$, $y_{j}) - v_{x}^{-}(z_{i}, y_{j})$ and $[v_{xx}] = v_{xx}^{+}(z_{i}, y_{j}) - v_{xx}^{-}(z_{i}, y_{j})$ represent jumps of $v$ and its partial derivatives. If $(x_i, y_j)$ is an irregular grid point in the subdomain $\Omega_e$, then the truncation error is

Thus, we can obtain a first-order approximation $C_{h, x}^{\pm}$ of the truncation error $E_{h, x}^{\pm}$ at the interface

Similarly, the approximation of the truncation error in the vertical direction can be calculated by

where $z_j$ is the intersection point along the $y$-direction, and $y_{j}\le z_{j} < y_{j+1}$. For irregular points, intersections of the five-point stencil and the interface are discussed separately in both horizontal and vertical directions, which are shown in Figure 2. The leading order terms of the truncation error in the $x$- and $y$-direction can be computed according to correction formulas (4.17) and (4.18). Thus, we have the following corrected finite difference scheme.

The modified system has the same coefficient matrix as the standard five-point compact finite difference equations, thus can be solved by the FFT-based fast elliptic solver. Numerical experiments show that correcting the local truncation errors, which are $O(h)$ at irregular grid points, is sufficient to guarantee global second-order accuracy of the discrete solution.

4.4.   Interpolation of the boundary integral values

Let $v({\mathbf p})$ be a smooth function on $\mathcal{B}\setminus \Gamma$, which is discontinuous only across $\Gamma$, and thus we need to recalculate the value of the function on the interface. Here we use polynomial interpolation because polynomial interpolation provides numerical stability. Assume that ${\mathbf q} \in \Gamma$ is a point on the interface, and then the Taylor expansion of $v({\mathbf p})$ around the point ${\mathbf q}$ yields

where $(\zeta, \eta)^{T} = {\mathbf p} - {\mathbf q}$. According to the Taylor expansions (4.20) and (4.21), the solution $v^{+}({\mathbf q})$ is obtained. There are six unknowns in the equation, and for the equation to have a unique solution, it is necessary to find six different points and perform a Taylor expansion at each of them to form a system of linear equations, which is then solved by the GMRES method.

For a given point $Q$ on the interface, the first step is to select the nearest grid point $P_1$ to $Q$. Then the four neighboring grid points $P_2, P_3, P_4,$ and $P_5$ of $P_1$ (above, below, left and right) are chosen as part of the six-point template. Finally, the diagonal grid point of $P_1$ is selected as the sixth point. When $Q$ lies on a grid line, the diagonal grid point to the left of the vector $\overrightarrow{P_1 Q}$, which is closer to $Q$, is chosen as the sixth point. When $Q$ lies exactly on a grid point, the diagonal grid point on the left side of the two outer diagonal grid points of $Q$ is selected as the sixth point, as shown in Figure 3.

For convenience, let

Thus, for each point ${\mathbf p}_{j}$ $(j = 0, 1, 2, 3, 4, 5)$, the Taylor expansion can be rewritten as

Since the points on the interface and the jumps in the partial derivatives of their first- and second-order are known, we can re-express the Taylor expansions for the points ${\mathbf p}_{j}$ outside the interface as

Combining the formulas (4.23) and (4.25) yields a system of linear equations for six points, which is then solved using either the GMRES method or the CG algorithm. The resulting $V^{+}$ is used as the function value at the points on the interface. We note that the jump calculation required in the correction and interpolation module is detailed in ^[42].

5.   Algorithm summary

In the proposed method, the reaction-diffusion equation is discretized by a two-stage SIRK method. At each time step, the resulting spatial equation is always a modified Helmholtz equation and can be solved by the KFBI method. The KFBI method first reformulates the elliptic boundary value problem into BIEs. Integrals involved in the BIEs are regarded as the solutions to some equivalent simple interface problems defined in the regular domain $\mathcal{B}$. Grid solutions of the equivalent interface problems are obtained by a three-step procedure, namely discretization, correction, and a fast solution, while the boundary evaluation of integrals needs a further interpolation process. The specific calculation process is summarized in Algorithm 2.

6.   Numerical example

This section presents numerical examples of the modified Helmholtz equation and the reaction-diffusion equation, both of which are solved using the KFBI method. For convenience, the square domain $\mathcal{B}$ is set to be $[-1.2, 1.2] \times [-1.2, 1.2]$ by default, and interfaces or boundaries are described by parametric curves or quadratic spline curves. Boundary integral equations involved in the algorithm are iteratively solved by the GMRES method and iteration ends when the residuals are less than $10^{-15}$. Tables 1–3 provide numerical results of different reaction-diffusion equations with different diffusion coefficients $D$. In each case, solution error $e_{h}$ is computed in both the maximum norm and the one norm where

Here $N$ is the number of points in the domain $\Omega$. Table 1 shows the results for the curve obtained by four spline interpolations for the scattered points. Tables 2 and 3 show the results based on the known parametric equations.

Tables 4–6 show the results calculated from the modified Helmholtz equation. The first two rows of the table show the infinite norm error $\left \| e_{h} \right \|_{\infty }$ and the one norm error $\left \| e_{h} \right \|_{1 }$, respectively, and the third row shows the CPU times (in seconds) on a MacBook Air laptop computer with an Apple M2 processor. Table 4 computes the error for the four-spline interpolation curve, and Tables 5 and 6 compute the error for the parametric curve. Table 7 shows a numerical example for the Allen-Cahn equation whose exact solution is unknown. We test the numerical accuracy by subtracting the numerical solution on a coarse grid by that on a corresponding refined grid, and computing the discrete absolute difference in both the maximum norm and one norm.

Example 1. This example solves the reaction-diffusion equation $\partial u / \partial t = D \Delta u + g(u)$ on the Dirichlet boundary, with the analytic solution and the corresponding $g(u)$ chosen as

The interface is obtained from 100 scattered points using four spline interpolations and the final numerical results are shown in Table 1. The heat map and contour map of the numerical results are shown in Figure 4, and the heat map of the error is shown in Figure 5.

Example 2. This example solves the reaction-diffusion equation $\partial u / \partial t = D \Delta u + g(u)$ on the Dirichlet boundary, where the solution region is bounded by parametric curves. The parametric equations are

with $c = 0.1$. The analytic solution in this case is

Numerical results are shown in Table 2, the heat map and contour plot of the numerical results are shown in Figure 6, and the heat map of the error is shown in Figure 7.

Example 3. This example solves the reaction-diffusion equation $\partial u / \partial t = D \Delta u + g(u)$ on the Dirichlet boundary, where the solution region is enclosed by parametric curves. The parametric equations are

with $c = 0.2$. The analytic solution is

Numerical results are shown in Table 3, the heat map and contour plot of the numerical results are shown in Figure 8, and the heat map and contour plot of the errors are shown in Figure 9.

Example 4. This example solves the modified Helmholtz equation $\Delta u(x, y) - \kappa u(x, y) = f(x, y)$, $(x, y) \in \Omega$, on the Dirichlet boundary, and the solution region $\Omega$ is obtained by interpolating 100 scattered points using four spline interpolations. The coefficient $\kappa$ is set to be $2.0$. The analytic solution chosen for this example is

Numerical results are shown in Table 4, and heat maps of the numerical results and errors are shown in Figure 10.

Example 5. This example solves the modified Helmholtz equation $\Delta u(x, y) - \kappa u(x, y) = f(x, y)$, $(x, y) \in \Omega$, on the Dirichlet boundary, where the solution region is bounded by parametric curves. The parametric equations are

with $c = 0.3$. The coefficient $\kappa$ is set to be $1.0$. The analytic solution in this case is

Numerical results are shown in Table 5, and the heat maps of the numerical results and errors are shown in Figure 11.

Example 6. This example solves the modified Helmholtz equation $\Delta u(x, y) - \kappa u(x, y) = f(x, y)$, $(x, y) \in \Omega$, on the Dirichlet boundary, where the solution region is bounded by parametric curves, and parametric equations are

with $c = 0.1$. The coefficient $\kappa$ is set to be $1.0$. The analytic solution in this example is the same as in Example 5. Numerical results are shown in Table 6, and the heat maps of the numerical results and errors are shown in Figure 12.

Example 7. This example solves the Allen-Cahn equation $\partial u / \partial t = D \Delta u - f(u)$ with $D = 0.01$ and $f(u) = u - u^3$ being the derivatives of the Ginzburg-Landau double-well potential. Suppose that the solution to the Allen-Cahn equation is subject to the initial condition $u_0(x, y) = 0.05 \sin(x) \sin(y)$, and the homogeneous Dirichlet boundary condition $u|_{\partial \Omega} = 0$. This example has been investigated in ^[54] where the problem domain is set to be a square. Here, we consider a general irregular problem domain. Suppose the domain boundary is given by

with $\alpha = 2.3$ and $\beta = 0.2$. The bounding box for the interface problem is set to be $[0, 2\pi]\times [0, 2\pi]$, consistent with the reference. To verify the numerical accuracy, we take the numerical solution obtained on $1024 \times 1024$ grid mesh with $\Delta t = 5 \times 10^{-3}$ as the exact solution and compare it with the corresponding approximation on a coarse grid and time step. The error in maximum norm and one norm at time $T$ = 1.2 are shown in Table 7, from which it is observed that the scheme achieves second-order accuracy. Figure 13 presents the heat maps of the numerical results and their corresponding errors.

7.   Discussion

This study primarily relies on the KFBI method for solving elliptic partial differential equations on the Dirichlet boundary, which transforms the boundary value problem into an interface problem, avoiding the need for complex integral computations. The entire computational process consists of four steps. The first step is to choose a suitable rectangular region to include the interface and establish a Cartesian grid, the second step is to correct for irregular points at the interface, the third step is to solve the problem, using the FFT to compute values at all grid points, and the fourth step is to interpolate. Polynomial interpolation is used to recalculate the function values at the interface. Finally, the boundary integral equations are solved using the GMRES method. This method is tested to enable numerical solution errors of second-order accuracy. This applies not only to the case where the interface is a parametric curve but also to quadratic spline interpolated curves. On this basis, further consideration is given to solving the reaction-diffusion equation, which is transformed into the form of two modified Helmholtz equations using the two-step semi-implicit Runge-Kutta method, and the final numerical solution error is obtained with the same second-order accuracy.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant No. DMS-12101553), and the Natural Science Foundation of Zhejiang Province (Grant No. LQ22A010017).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.