Latest revision as of 14:33, 4 January 2024

Abstract

This paper introduces a sixth-order Immersed Interface Method (IIM) for addressing 2D Poisson problems characterized by a discontinuous forcing function with straight interfaces. In the presence of this discontinuity, the problem exhibits a non-smooth solution at the interface that divides the domain into two regions. Here, the IIM is employed to compute the solution on a fixed Cartesian grid. This method integrates necessary jump conditions resulting from the interface into the numerical schemes. In order to achieve a sixth-order method, the proposed approach combines implicit finite differences with the IIM. The proposed scheme is efficient because the matrix arising from discretization remains the same as in the smooth problem, and changes are made to the resulting linear system by introducing new terms on the right side. These supplementary terms account for the discontinuities in the solution and its derivatives, with calculations restricted near the interface. The paper demonstrates the accuracy of the proposed method through various numerical examples.

Keywords: Poisson equation, immersed interface method, finite difference, sixth-order of accuracy, implicit finite difference

1. Introduction

Developing advanced algorithms for solving the Poisson equation holds great significance in numerous research domains, including computational fluid dynamics, wave propagation, and theoretical physics [1]. The discontinuous problem emerges in scenarios marked by abrupt changes at interfaces that separate different regions within a given domain [2,3,4]. In this article, we particularly focus on linear interfaces, which are commonly encountered in layered phenomena (for example, as seen in [5,6,7] and the cited references). Additionally, the pursuit of high-order methods for addressing such challenges is highly advantageous since their enhanced precision permits the use of coarser grids, subsequently reducing computational expenses.

This paper presents high-order finite-difference schemes up to sixth-order for the Poisson equation for straight interfaces. The problem is given by [2,3,4]

where ${\textstyle u}$ and ${\textstyle f}$ are the solution and known right-hand side function, respectively. We divide ${\textstyle \Omega }$ in two regions, ${\textstyle \Omega ^{+}}$ and ${\textstyle \Omega ^{-}}$, separated by an immersed interface ${\textstyle \Gamma }$. The computational domain can be one-dimensional (1D) (with several interface points), or a two-dimensional (2D) region with straight interfaces. We use Dirichlet boundary conditions on ${\textstyle \partial \Omega }$. We assume that the solution, the right-hand side function, and their derivatives may have discontinuities at ${\textstyle \Gamma }$. Thus, we require jump conditions as additional inputs. The principal jump conditions ${\textstyle [u]={\mathfrak {p}}}$ and ${\textstyle [u_{n}]={\mathfrak {q}}}$ are known functions and are defined on ${\textstyle \Gamma }$. Here, ${\textstyle u_{n}}$ is the derivative in the normal direction.

Many numerical methods have been proposed to solve accurately the Poisson equation; however most of these methods are limited to smooth solutions. For instance, many developments have been made to get fourth- and sixth-order finite-difference methods for the Poisson equation (1), see for example [8,9,10,11,12,13,14]. On the other hand, to overcome the interface issue, several approaches exist to approximate discontinuous solutions, such as the level set method [15], immersed boundary method [16], ghost fluid method [17], interpolation matched interface method [18], Galerkin finite element method [19], boundary condition capturing [20], and interface neural networks [21]. However, there are not many available high-order discretizations of the Poisson problems with interfaces and many of these algorithms are only second-order accuracy. From these methods, the immersed interface method [4,22,23,24,25] is one popular option to solve Eqs. (1)-(3) accurately by simple modifications of standard finite differences.

There only exist a few high-order immersed interface methods to solve elliptic equations. Fourth-order interface methods for elliptic equations with discontinuous solutions or discontinuous coefficients are investigated in [26,27,28,29,30]. Lately, Feng and Li [31] presented a third-order IIM for elliptic interface problems, but it is limited to straight interfaces lying at grid points. Pan et al. [32] proposed a third-order IIM to solve elliptic problems on irregular domains, and Colnago et al. [33] developed a fourth-order approximation. More recently, Feng et al. [34] presented a sixth-order IIM to solve Poisson interface problem with singular sources based on the undetermined coefficients technique.

This paper presents a new high-order implicit finite-difference immersed interface method (HIFD-IIM) up to sixth-order accuracy to solve the 2D Poisson problem (1)-(3) with straight interfaces. The implicit finite-difference methodology is based on calculating the unknown variable and its corresponding derivatives simultaneously [35,36]. The system is changed by adding new terms on the right side, named jump contributions. These terms include the jumps of the solution and its derivatives, and they are calculated near the interface. In this context, the straight interfaces allows to calculate the jumps contributions directly from principal jump conditions (2) and (3) without any other calculation. The modifications are only performed at grid points where the method's stencil intersects the interface. Moreover, the matrix in the linear system is the same as the smooth problem. This formulation makes the method attractive as it is easy to implement and does not require other convergence restrictions than the ones from standard methods for smooth solutions. To the best authors' knowledge, there is not other methods using the proposed formulation.

The paper is organized as follows. Section 2 introduces the implicit formulation. The main theoretical results are presented in this section. Section 3 deals with the 1D Poisson problem. Section 4 shows how to implement the high order methods for 2D problems with straight interfaces. Section 5 contains several examples to test the algorithm's capacity. Finally, in Section 6, we present the conclusions and future work.

2. Implicit finite-difference formulation

We begin our analysis by considering the 1D finite-difference (FD) scheme for the second derivative of a real-valued function ${\textstyle u}$. In this work, the numerical solution is approximated using a uniform grid. An interval ${\textstyle [{\mathfrak {a}},\,{\mathfrak {b}}]}$ is divided into ${\textstyle N}$ sub-intervals, using ${\textstyle x_{i}={\mathfrak {a}}+(i-1)h}$, ${\textstyle i=1,2,\dots ,N,N+1}$, where the grid size is given by ${\textstyle h=({\mathfrak {b}}-{\mathfrak {a}})/N}$. For simplicity, we set ${\textstyle x=x_{\alpha }}$ between two consecutive grid points ${\textstyle x_{I}\leq x_{\alpha }{<}x_{I+1}}$. We called ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ irregular points; meanwhile the reminder points are referred as regular. Besides ${\textstyle h}$, the discretization needs the definitions of ${\textstyle h_{L}=x_{I}-x_{\alpha }}$, and ${\textstyle h_{R}=x_{I+1}-x_{\alpha }}$. Note that ${\textstyle h_{R}}$ is a positive value, meanwhile ${\textstyle h_{L}}$ is negative (Figure 1).

Figure 1. Stencil of the immersed interface method at ${\displaystyle x_{I}}$ and ${\displaystyle x_{I+1}}$

Finally, we denote as ${\textstyle \varphi _{i}=\left.\varphi \right|_{i}=\varphi (x_{i})}$ the evaluation of function ${\textstyle \varphi }$ at the ${\textstyle i}$th-point of the grid.

2.1 Approximation on regular points

It is well-known that the central finite difference gives us a second-order approximation, and can be written as

where

for a given small value ${\textstyle h>0}$.

One way to increase the order of accuracy in Eq. (4) is using a FD operator with more grid points. However, there is another way to get high-order approximations without changing the length of the FD operator (5). It is by considering a high-order implicit finite-difference (HIFD) formulation, as presented in Zapata and Balam [13]. Thus, using Taylor series expansions, it directly follows that the sixth-order formula is given by

where the FD operator ${\textstyle \delta ^{2}}$ is given by Eq. (5) and the partial operator is defined as

${\textstyle b=1/12}$ and ${\textstyle e=1/360}$. The above formulation (6) turns in a fourth-order method if we choose ${\textstyle b=1/12}$ and ${\textstyle e=0}$. Moreover, if we select ${\textstyle b=0}$ and ${\textstyle e=0}$ the standard second-order approximation is recovered.

We remark that Eq. (6) can only be applied in cases where the problem's solution has enough regularity [37]. To overcome this issue, we combine the HIFD with the IIM to solve problems with discontinuous solutions, as described in the next section.

2.2 Approximation on irregular points

This paper proposes a new formulation named HIFD-IIM that is based on the combination of HIFD and IIM. To derive high-order schemes, the IIM requires additional conditions at ${\textstyle i=I,I+1}$. These are known as jump conditions at ${\textstyle x_{\alpha }}$. Thus, for a function ${\textstyle \varphi }$, they are given by

The IIM contribution requires to include high-order jump derivatives. However, we can reduce the number of jumps by applying a less accurate scheme at the irregular points. As other IIMs proposed by different authors [22,24,38,39,40], the global order is ${\textstyle O(h^{2})}$, even if the local truncation error at ${\textstyle i=I}$ and ${\textstyle i=I+1}$ is one order lower, i.e., ${\textstyle O(h)}$ [41]. Recently, Pan et al. [32] proposed a global third-order IIM using ${\textstyle O(h^{3})}$ and ${\textstyle O(h^{2})}$ local truncation errors at regular and irregular grid points, respectively. Following similar ideas, the main result of this paper is presented in Theorem 1.

Theorem 1: HIFD-IIM. Let us consider the known jump conditions

on ${\textstyle x_{\alpha }}$ such that ${\textstyle x_{I}\leq x_{\alpha }{<}x_{I+1}}$. Then ${\textstyle u_{xx}}$ can be approximated at ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ by the finite-difference scheme

where

We obtain a fifth-order scheme for at ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ following similar ideas as the generalized Taylor expansion proposed by Xu and Wang [38] and the IIM for elliptic interface problems with straight interfaces proposed by Feng and Li [39].

We initially consider extended solutions of ${\textstyle u}$, named ${\textstyle u_{\ell }}$ and ${\textstyle u_{r}}$, as shown in Figure 2. The idea is to have smooth functions such that we can apply the standard central scheme to ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$. These functions are defined as

where ${\textstyle u^{-}}$ and ${\textstyle u^{+}}$ are defined in Eq. (8).

Figure 2. Extended solutions ${\displaystyle u_{\ell }}$ and ${\displaystyle u_{r}}$

Taylor series expansions of ${\textstyle u_{I+1}}$ around ${\textstyle x_{\alpha }}$ yields

Using the definition of jumps (9), it follows

Next, we use ${\textstyle u_{\ell }}$ definition in (13) to obtain

Thus,

where ${\textstyle C_{I}}$ is defined as Eq. (11). Substituting (15) into a standard central scheme for the second-order derivative, it yields

On the other hand, using Taylor series of ${\textstyle (u_{\ell })_{xx}}$, we have

Thus,

Finally, we get ${\textstyle \left.\partial ^{4}u_{xx}\right|_{I}=\left.\delta _{x}^{2}u\right|_{I}+C_{I}+O(h^{5})}$. This completes the proof. The same procedure can be applied for the proof at ${\textstyle x_{I+1}}$ using definition (14).

It is important to remark that ${\textstyle C_{I}}$ and ${\textstyle C_{I+1}}$ are constants computed from the jump derivatives of ${\textstyle u}$ and we assume that those values are known. To emphasize that the contribution includes all jump derivatives up to sixth-order we write ${\textstyle C_{I}^{6}}$ instead ${\textstyle C_{I}}$.

Remark 1: We can rewrite the contributions (11) and (12) as follows

where ${\textstyle b}$ and ${\textstyle e}$ are defined as in Eq. (7).

Remark 2: If the solution is smooth, then all contributions ${\textstyle C_{I}^{6}}$ and ${\textstyle C_{I+1}^{6}}$ are equal to zero and the standard sixth-order method [37] is recovered in Eq. (6).

Corollary 1: If we consider ${\textstyle b=1/12}$ and ${\textstyle e=0}$ in Eqs. (1) and (17), then we obtain a third-order scheme for ${\textstyle i=I,I+1}$ (global fourth-order method) as follows

where

Corollary 2: If ${\textstyle b=e=0}$ in Eqs. (1) and (17), the method represents an explicit finite-difference scheme of first-order of accuracy for ${\textstyle i=I,I+1}$ (global second-order method) as follows

where

In this case, we only require to explicitly know jump conditions ${\textstyle [u]}$, ${\textstyle \left[u_{x}\right]}$, and ${\textstyle \left[u_{xx}\right]}$.