A new sixth-order immersed interface method for solving Poisson equations with straight interfaces

<!-- metadata commented in wiki content
==A new sixth-order immersed interface method for solving Poisson equations with straight interfaces==

'''R. Itza Balam<sup>a,b</sup>, M. Uh Zapata<sup>a,b</sup>, J. Montalvo-Urquizo<sup>c</sup>'''

{|
|-
|<sup>a</sup> Centro de Investigación en Matemáticas A. C., CIMAT Unidad Mérida, México
|}

{|
|-
|<sup>b</sup> Consejo Nacional de Humanidades, Ciencias y Tecnologías, CONAHCYT, México
|}

{|
|-
|<sup>c</sup> School of Science and Engineering, Tecnológico de Monterrey
|}
-->
==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 <span id='citeF-1'></span>[[#cite-1|[1]]]. The discontinuous problem emerges in scenarios marked by abrupt changes at interfaces that separate different regions within a given domain <span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-2|[2,3,4]]]. In this article, we particularly focus on linear interfaces, which are commonly encountered in layered phenomena (for example, as seen in <span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span>[[#cite-5|[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 <span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-2|[2,3,4]]]

<span id="eq-1"></span>
<span id="eq-2"></span>
<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\Delta u({x}) = f({x}), \quad {x}\in \Omega \setminus \Gamma , \quad \Omega = \Omega ^{+} \cup \Omega ^{-} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \left[u\right]_\Gamma   = \mathfrak{p}({x}),  \quad {x}\in \Gamma , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math> \left[u_{{n}} \right]_\Gamma = \mathfrak{q}({x}) , \quad {x}\in \Gamma{.}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math display="inline">u</math> and <math display="inline">f</math> are the solution and known right-hand side function, respectively. We divide <math display="inline">\Omega </math> in two regions, <math display="inline">\Omega ^{+}</math> and <math display="inline">\Omega^{-}</math>, separated by an immersed interface <math display="inline">\Gamma </math>. 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 <math display="inline">\partial \Omega </math>. We assume that the solution, the right-hand side function, and their derivatives may have discontinuities at <math display="inline">\Gamma </math>. Thus, we require jump conditions as additional inputs. The principal jump conditions <math display="inline">[u]=\mathfrak{p}</math> and <math display="inline">[u_{{n}}]=\mathfrak{q}</math> are known functions and are defined on <math display="inline">\Gamma </math>. Here, <math display="inline">u_{{n}}</math> 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 [[#eq-1|(1)]], see for example <span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span>[[#cite-8|[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 <span id='citeF-15'></span>[[#cite-15|[15]]], immersed boundary method <span id='citeF-16'></span>[[#cite-16|[16]]], ghost fluid method <span id='citeF-17'></span>[[#cite-17|[17]]], interpolation matched interface method <span id='citeF-18'></span>[[#cite-18|[18]]], Galerkin finite element method <span id='citeF-19'></span>[[#cite-19|[19]]], boundary condition capturing <span id='citeF-20'></span>[[#cite-20|[20]]], and interface neural networks <span id='citeF-21'></span>[[#cite-21|[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 <span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-4'></span><span id='citeF-25'></span>[[#cite-22|[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 <span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-26|[26,27,28,29,30]]]. Lately, Feng and Li <span id='citeF-31'></span>[[#cite-31|[31]]] presented a third-order IIM for elliptic interface problems, but it is limited to straight interfaces lying at grid points. Pan et al. <span id='citeF-32'></span>[[#cite-32|[32]]] proposed a third-order IIM to solve elliptic problems on irregular domains, and Colnago et al. <span id='citeF-33'></span>[[#cite-33|[33]]] developed a  fourth-order approximation. More recently, Feng et al. <span id='citeF-34'></span>[[#cite-34|[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 [[#eq-1|(1)]]-[[#eq-3|(3)]] with straight interfaces. The implicit finite-difference methodology is based on calculating the unknown variable and its corresponding derivatives simultaneously  <span id='citeF-35'></span><span id='citeF-36'></span>[[#cite-35|[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 [[#eq-2|(2)]] and [[#eq-3|(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 <math display="inline">u</math>. In this work, the numerical solution is approximated using a uniform grid.  An interval <math display="inline">[\mathfrak{a},\,\mathfrak{b}]</math> is divided into <math display="inline">N</math> sub-intervals, using <math display="inline">x_i = \mathfrak{a} + (i-1)h</math>, <math display="inline">i=1, 2, \dots , N,N+1</math>, where the grid size is given by <math display="inline">h=(\mathfrak{b}-\mathfrak{a})/N</math>. For simplicity, we set <math display="inline">x=x_\alpha </math> between two consecutive grid points <math display="inline">x_I\leq x_\alpha{<}x_{I+1}</math>. We called <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> irregular points; meanwhile the reminder points are referred as regular. Besides <math display="inline">h</math>, the discretization needs the definitions of <math display="inline">h_{L} = x_{I}-x_\alpha </math>, and <math display="inline">h_{R} = x_{I+1}-x_\alpha </math>. Note that <math display="inline">h_{R}</math> is a positive value, meanwhile <math display="inline">h_{L}</math> is negative (Figure [[#img-1|1]]).

<div id='img-1'></div>
<div id='img-1'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Balam_338770048-Stencil.png|350px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 1'''. Stencil of the immersed interface method at <math>x_{I}</math> and <math>x_{I+1}</math>
|}


Finally, we denote as <math display="inline">\varphi _i = \left.\varphi \right|_i = \varphi (x_i)</math> the evaluation of function <math display="inline">\varphi </math> at the <math display="inline">i</math>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

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\left.u_{xx}\right|_i = \left.\delta ^2u\right|_i + O(h^2),   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where
<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\left.\delta ^2u\right|_i := \frac{u_{i-1}-2u_{i}+u_{i+1}}{h^2},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

for a given small value <math display="inline">h > 0</math>.

One way to increase the order of accuracy in Eq. [[#eq-4|(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 [[#eq-5|(5)]]. It is by considering a high-order implicit finite-difference (HIFD) formulation, as presented in Zapata and Balam <span id='citeF-13'></span>[[#cite-13|[13]]]. Thus, using Taylor series expansions, it directly follows that the sixth-order formula is given by

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\left.\partial ^4u_{xx}\right|_i = \left.\delta ^2 u\right|_i + O(h^6), \qquad i\ne I, I+1.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where the FD operator <math display="inline">\delta ^2</math> is given by Eq. [[#eq-5|(5)]] and the partial operator is defined as

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\partial ^4(\cdot )  :=  (\cdot )  + bh^2(\cdot )_{xx} + eh^4 (\cdot )_{xxxx},   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

<math display="inline">b=1/12</math> and <math display="inline">e = 1/360</math>. The above formulation [[#eq-6|(6)]] turns in a fourth-order method if we choose <math display="inline">b=1/12</math> and <math display="inline">e=0</math>. Moreover, if we select <math display="inline">b=0</math> and <math display="inline">e=0</math> the standard second-order approximation is recovered.

We remark that Eq. [[#eq-6|(6)]] can only be applied in cases where the problem's solution has enough regularity <span id='citeF-37'></span>[[#cite-37|[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 <math display="inline">i=I, I+1</math>. These are known as jump conditions at <math display="inline">x_\alpha </math>. Thus, for a function <math display="inline">\varphi </math>, they are given by

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\left[\varphi \right]= \varphi ^{+} - \varphi ^{-},  \qquad  \varphi ^{+} = \varphi (x_\alpha ^{+} )= \lim _{x\to x_\alpha ^+}\varphi (x), \qquad  \varphi ^{-} = \varphi (x_\alpha ^{-}) = \lim _{x\to x_\alpha ^-}\varphi (x).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

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 <span id='citeF-22'></span><span id='citeF-24'></span><span id='citeF-38'></span><span id='citeF-39'></span><span id='citeF-40'></span>[[#cite-22|[22,24,38,39,40]]], the global order is <math display="inline">O(h^2)</math>, even if the local truncation error at <math display="inline">i=I</math> and <math display="inline">i=I+1</math> is one order lower, i.e., <math display="inline">O(h)</math> <span id='citeF-41'></span>[[#cite-41|[41]]]. Recently, Pan et al. <span id='citeF-32'></span>[[#cite-32|[32]]] proposed a global third-order IIM using <math display="inline">O(h^{3})</math> and <math display="inline">O(h^{2})</math> local truncation errors at regular and irregular grid points, respectively. Following similar ideas, the main result of this paper is presented in Theorem [[#theorem-thm:xj_right_left|1]].

<span id='theorem-thm:xj_right_left'></span>'''Theorem 1:  HIFD-IIM'''. Let us consider the known jump conditions

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\left[u \right],\quad  \left[u_{x} \right],\quad  \left[u_{xx} \right],\quad  \left[u_{xxx} \right],\quad  \left[u_{xxxx} \right],\quad  \left[u_{xxxxx} \right],\quad  \left[u_{xxxxxx} \right],\quad   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

on <math display="inline">x_\alpha </math> such that <math display="inline">x_I\leq x_\alpha{<}x_{I+1}</math>. Then <math display="inline">u_{xx}</math> can be approximated at <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> by the finite-difference scheme

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>\left.\partial ^4 u_{xx}\right|_i = \left.\delta ^2 u\right|_i + C_i + O(h^5), \quad i= I, I+1,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

where
<span id="eq-11"></span>
<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math> C_I     :=  - \dfrac{1}{h^2}\left\{[u] + h_R\left[u_{x}\right]+ \dfrac{h_R^2}{2}\left[u_{xx} \right]+ \dfrac{h_R^3}{3!}\left[u_{xxx} \right]+\dfrac{h_R^4}{4!}\left[u_{xxxx} \right]+ \dfrac{h_R^5}{5!}\left[u_{xxxxx} \right]+ \dfrac{h_R^6}{6!}\left[u_{xxxxxx} \right]\right\}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math> C_{I+1} :=    \dfrac{1}{h^2}\left\{[u]+h_L\left[u_{x} \right]+\dfrac{h_L^2}{2}\left[u_{xx} \right]+ \dfrac{h_L^2}{3!}\left[u_{xxx} \right]+\dfrac{h_L^4}{4!}\left[u_{xxxx} \right]+ \dfrac{h_L^6}{5!}\left[u_{xxxxx} \right]+\dfrac{h_L^6}{6!}\left[u_{xxxxxx} \right]\right\}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}


We obtain a fifth-order scheme for at <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> following similar ideas as the generalized Taylor expansion proposed by Xu and Wang <span id='citeF-38'></span>[[#cite-38|[38]]] and the IIM for elliptic interface problems with straight interfaces proposed by Feng and Li <span id='citeF-39'></span>[[#cite-39|[39]]].

We initially consider extended solutions of <math display="inline">u</math>, named <math display="inline">u_\ell </math> and <math display="inline">u_r</math>, as shown in Figure [[#img-2|2]]. The idea is to have smooth functions such that we can apply the standard central scheme to <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math>. These functions are defined as

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_\ell (x) = \begin{cases}u(x), & x\leq x_\alpha ,\\ \displaystyle u^-+h_R u_{x}^{-}+\frac{h_R^2}{2}u_{xx}^{-} +\dfrac{h_R^3}{3!}u_{xxx}^{-}+\dfrac{h_R^4}{4!}u_{xxxx}^{-}+\dfrac{h_R^5}{5!}u_{xxxxx}^{-}+\dfrac{h_R^6}{6!}u_{xxxxxx}^{-} , & x_\alpha{<}x,\\ \end{cases}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_r(x) = \begin{cases}\displaystyle u^++h_Lu_{x}^{+} +\frac{h_L^2}{2}u_{xx}^{+} +\dfrac{h_L^3}{3!}u_{xxx} ^{+}+\dfrac{h_L^4}{4!}u_{xxxx}^{+}+\dfrac{h_L^5}{5!}u_{xxxxx}^{+}+\dfrac{h_L^6}{6!}u_{xxxxxx}^{+}, & x<x_\alpha ,\\ u(x), & x_\alpha \leq x, \end{cases}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">u^{-}</math> and <math display="inline">u^{+}</math> are defined as in Eq.  [[#eq-8|(8)]].

<div id='img-2'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|align="center" | 
{|style="margin: 0em auto 0.1em auto;width:auto;" 
|+
|-
|style="text-align: center;padding:10px;"| [[Image:Draft_Balam_338770048-Function_ul.png|344px|]]
|style="text-align: center;padding:10px;"| [[Image:Draft_Balam_338770048-Function_ur.png|344px|]]
|}
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 2'''. Extended solutions <math>u_\ell </math> and <math>u_r</math>
|}


Taylor series expansions of <math display="inline">u_{I+1}</math> around <math display="inline">x_\alpha </math> yields

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_{I+1} = u^{+}+h_Ru_x^{+}+\dfrac{h_R^2}{2}u_{xx}^{+}+\dfrac{h_R^3}{3!}u_{xxx} ^{+}+\dfrac{h_R^4}{4!}u_{xxxx}^{+} +\dfrac{h_R^5}{5!}u_{xxxxx}^{+} +\dfrac{h_R^6}{6!}u_{xxxxxx}^{+}+O(h^7). </math>
|}
|}

Using the definition of jumps [[#eq-9|(9)]], it follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_{I+1} = \left(u^-+h_Ru_x^-+\dfrac{h_R^2}{2}u_{xx} ^-+\dfrac{h_R^3}{3!}u_{xxx} ^-+\dfrac{h_R^4}{4!}u_{xxxx}^-+\dfrac{h_R^5}{5!}u_{xxxxx}^-+\dfrac{h_R^6}{6!}u_{xxxxxx}^-\right)</math>
|-
| style="text-align: center;" | <math> +\;[u]+h_R\left[u_{x} \right]+\dfrac{h_R^2}{2}\left[u_{xx} \right]+\dfrac{h_R^3}{3!}\left[u_{xxx} \right]+\dfrac{h_R^4}{4!}\left[u_{xxxx} \right]+\dfrac{h_R^5}{5!}\left[u_{xxxxx} \right]+\dfrac{h_R^6}{6!}\left[u_{xxxxxx} \right]+O(h^7). </math>
|}
|}

Next, we use <math display="inline">u_\ell </math> defintion in [[#eq-13|(13)]] to obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_{I+1} = {u_\ell }_{I+1}+[u]+h_R\left[u_{x} \right]+\dfrac{h_R^2}{2}\left[u_{xx} \right]+\dfrac{h_R^3}{3!}\left[u_{xxx} \right]+\dfrac{h_R^4}{4!}\left[u_{xxxx} \right]+\dfrac{h_R^5}{5!}\left[u_{xxxxx} \right]+\dfrac{h_R^6}{6!}\left[u_{xxxxxx} \right]+O(h^7). </math>
|}
|}

Thus,

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_{I+1} = {u_\ell }_{I+1} - h^2C_I +O(h^7),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">C_I</math> is defined as equation [[#eq-11|(11)]]. Substituting [[#eq-15|(15)]] into a standard central scheme for the second-order derivative, it yields

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\delta ^2 u\right|_I   =\frac{1}{h^2}u_{I+1}-\dfrac{2}{h^2}u_{I}+\frac{1}{h^2}u_{I-1} = \frac{1}{h^2} u_{\ell I+1}-\frac{2}{h^2}{u_\ell }_I +\frac{1}{h^2}{u_\ell }_{I-1} - C_{I} + O(h^5). </math>
|}
|}

On the other hand, using Taylor series of <math display="inline">(u_\ell )_{xx}</math>, we have

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.(u_\ell )_{xx}\right|_I  =\frac{1}{h^2} u_{\ell I+1}-\frac{2}{h^2}{u_\ell }_I +\frac{1}{h^2}{u_\ell }_{I-1} - bh^2(u_\ell )_{xxxx}- eh^4(u_\ell )_{xxxxxx} +O(h^6).  </math>
|}
|}

Thus,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\delta ^2 u\right|_I   =  \left.(u_\ell )_{xx} \right|_I + \left.bh^2 (u_\ell )_{xxxx}\right|_I + \left.eh^4 (u_\ell )_{xxxxxx}\right|_I - C_I+O(h^5). </math>
|}
|}


Finally, we get <math display="inline">\left.\partial ^4 u_{xx}\right|_I = \left.\delta ^2_x u\right|_I + C_I + O(h^5)</math>. This completes the proof.  The same procedure can be applied for the proof at <math display="inline">x_{I+1}</math> using definition [[#eq-14|(14)]].

It is important to remark that <math display="inline">C_I</math> and <math display="inline">C_{I+1}</math> are constants computed from the jump derivatives of <math display="inline">u</math> and we assume that those values are known. To emphasize that the contribution includes all jump derivatives up to sixth-order we write <math display="inline">C_I^6</math> instead <math display="inline">C_I</math>.

''Remark 1'': We can rewrite the contributions [[#eq-11|(11)]] and [[#eq-12|(12)]] as follows

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>C^6_I     :=- \dfrac{1}{h^2}\left\{[u] + h_R\left[u_{x}\right]+ \dfrac{h_R^2}{2}\left[u_{xx} \right] + b\left(2h_R^3\left[u_{xxx} \right]+\dfrac{h_R^4}{2}\left[u_{xxxx} \right]\right)\right.</math>
|-
| style="text-align: center;" | <math>  \qquad  \quad \left.+\,e\left(3h_R^5\left[u_{xxxxx} \right]+ \dfrac{h_R^6}{2}\left[u_{xxxxxx} \right]\right)\right\}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math> C^6_{I+1} :=    \dfrac{1}{h^2}\left\{[u] + h_L\left[u_{x}\right]+ \dfrac{h_L^2}{2}\left[u_{xx} \right]+ b\left(2h_L^3\left[u_{xxx} \right]+\dfrac{h_L^4}{2}\left[u_{xxxx} \right]\right)\right.</math>
|-
| style="text-align: center;" | <math>   \qquad \quad \left.+\,e\left(3h_L^5\left[u_{xxxxx} \right]+ \dfrac{h_L^6}{2}\left[u_{xxxxxx} \right]\right)\right\},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where <math display="inline">b</math> and <math display="inline">e</math> are defined as in Eq. [[#eq-7|(7)]].

''Remark 2'': If the solution is smooth, then all contributions <math display="inline">C^6_I</math> and <math display="inline">C^6_{I+1}</math> are equal to zero and the standard sixth-order method <span id='citeF-37'></span>[[#cite-37|[37]]] is recovered in Eq. [[#eq-6|(6)]].

''Corollary 1'': If we consider <math display="inline">b=1/12</math> and <math display="inline">e=0</math> in Eqs. [[#theorem-def:HIFD-IIM_CI_be|(1)]] and [[#eq-17|(17)]], then we obtain a third-order scheme for <math display="inline"> i= I, I+1</math> (global fourth-order method) as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.u_{xx}\right|_i  + \left.bh^2 (u_{xx})_{xx}\right|_i = \left.\delta ^2 u\right|_i + C^4_i + O(h^3), \quad i=I,I+1,    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

where
<span id="eq-19"></span>
<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>C^4_I     = - \dfrac{1}{h^2}\left\{[u] + h_R\left[u_{x} \right]+ \dfrac{h_R^2}{2}\left[u_{xx} \right]+ \dfrac{h_R^3}{3!}\left[u_{xxx} \right]+ \dfrac{h_R^4}{4!}\left[u_{xxxx} \right]\right\}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math> C^4_{I+1} =    \dfrac{1}{h^2}\left\{[u]+h_L\left[u_{x} \right]+ \dfrac{h_L^2}{2}\left[u_{xx} \right]+ \dfrac{h_L^3}{3!}\left[u_{xxx} \right]+ \dfrac{h_L^4}{4!}\left[u_{xxxx} \right]\right\}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

''Corollary 2'': If <math display="inline">b = e = 0</math> in Eqs. [[#theorem-def:HIFD-IIM_CI_be|(1)]] and [[#eq-17|(17)]], the method represents an explicit finite-difference scheme of first-order of accuracy for <math display="inline"> i= I, I+1</math> (global second-order method) as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.u_{xx}\right|_i = \left.\delta ^2 u\right|_i + C^2_i+ O(h), \quad i=I,I+1,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

where
<span id="eq-22"></span>
<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math>C^2_I    = - \dfrac{1}{h^2}\left\{[u]+h_R\left[u_{x} \right]+\dfrac{h_R^2}{2}\left[u_{xx} \right]\right\}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;" 
|-
| style="text-align: center;" |  <math> C^2_{I+1} =    \dfrac{1}{h^2}\left\{[u]+h_L\left[u_{x} \right]+\dfrac{h_L^2}{2}\left[u_{xx} \right]\right\}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

In this case, we only require to explicitly know jump conditions <math display="inline">[u]</math>, <math display="inline">\left[u_{x} \right]</math>, and <math display="inline">\left[u_{xx} \right]</math>.

Observe that previous corollaries include the superscript <math display="inline">2</math> and <math display="inline">4</math> to emphasize contributions upto second- and fourth-order derivatives, respectively.

==3. One-dimensional problem==

In this section, we use <math display="inline">u_{xx}</math> approximation from Theorem [[#theorem-thm:xj_right_left|1]] to study the 1D Poisson problem given by

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_{xx}(x) = f(x),\quad x\in \Omega = (\mathfrak{a},x_\alpha ) \cup (x_\alpha , \mathfrak{b}),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

where <math display="inline">u</math> and <math display="inline">f</math> can be discontinuous functions at a given point <math display="inline">x=x_\alpha </math>, and the principal jump conditions <math display="inline">\left[u\right] = \mathfrak{p}</math>, <math display="inline">\left[u_{x}\right] = \mathfrak{q}</math> are known values at <math display="inline">x_\alpha </math>. For the boundary conditions, we impose the Dirichlet type. Although this technique can be applied to several interface points, we only focus on one point <math display="inline">x_\alpha \in [x_{I},x_{I+1})</math> to simplify our exposition.

If we apply the partial operator [[#eq-7|(7)]] at both sides of [[#eq-24|(24)]], then we get

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\partial ^4 u_{xx}\right|_i= \partial ^4 f\left.\right|_i, \qquad i=2, \dots ,I,I+1,\dots , N.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

Substituting Eqs. [[#eq-6|(6)]] and [[#eq-10|(10)]] into the left hand-side of Eq. [[#eq-25|(25)]], we get

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\partial ^4 u_{xx}\right|_i = \left.\delta ^2 u\right|_i+ C^6_i\{ u\} , \qquad i=2, \dots ,I,I+1,\dots , N,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

where <math display="inline">C^6_i</math> corresponds to the contribution term of <math display="inline">u</math> at <math display="inline">x_i</math> given by

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>C^6_i\{ u\} =  \begin{cases}- \dfrac{1}{h^2}\left\{[u] + h_R\left[u_{x}\right]+ \dfrac{h_R^2}{2}\left[u_{xx} \right]+ \dfrac{h_R^3}{3!}\left[u_{xxx} \right]+\dfrac{h_R^4}{4!}\left[u_{xxxx} \right]\right.\\ \left.\quad \quad \quad + \dfrac{h_R^5}{5!}\left[u_{xxxxx} \right]+ \dfrac{h_R^6}{6!}\left[u_{xxxxxx} \right]\right\},& i=I,\\[4mm]   \dfrac{1}{h^2}\left\{[u]+h_L\left[u_{x} \right]+\dfrac{h_L^2}{2}\left[u_{xx} \right]+ \dfrac{h_L^2}{3!}\left[u_{xxx} \right]+\dfrac{h_L^4}{4!}\left[u_{xxxx} \right]\right.\\ \left.\quad \quad \quad + \dfrac{h_L^6}{5!}\left[u_{xxxxx} \right]+\dfrac{h_L^6}{6!}\left[u_{xxxxxx} \right]\right\},& i=I+1,\\    0, & \mathrm{otherwise.}  \end{cases}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}


Here, we introduce the notation <math display="inline">C^6_i\{ u\} </math> to emphasize that the contribution depends on the high-order jump derivatives and it is computed from the function <math display="inline">u</math>. It is necessary because in next section we will require to obtain contributions from different functions.

If we explicitly know the function <math display="inline">f</math> and its derivatives, then the right-hand side of Eq. [[#eq-25|(25)]] can be calculated as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\partial ^4 f\right|_i=f_i + bh^2 f_{xx}|_i+eh^4 f_{xxxx}|_i. </math>
|}
|}

On the other hand, if we only know values of <math display="inline">f</math> on the grid, we have to approximate the second- and fourth-order derivative of <math display="inline">f</math> at <math display="inline">x_i</math>.  As with other implicit schemes <span id='citeF-13'></span><span id='citeF-36'></span>[[#cite-13|[13,36]]], the right-hand side derivatives are calculated using a central finite-difference method. The discretization of the right-hand side for <math display="inline">i\neq I, I+1</math> is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\partial ^4 f\right|_i=f_i +\left.bh^2 f_{xx}\right|_i+\left.eh^4 f_{xxxx}\right|_i</math>
|-
| style="text-align: center;" | <math> =\,df_{i-2}+ (b-4d)f_{i-1}+\,(1-2b+6d)f_i+\,(b-4d)f_{i+1}+\,df_{i+2}+O(h^5), </math>
|}
|}

where <math display="inline">d =: e-b^2 = -1/240</math>. For <math display="inline">i = I, I+1</math>, we need to compute the derivatives using the IIM technique which is described as follows. Notice that the second derivative of <math display="inline">f</math> in <math display="inline">bh^2 f_{xx}</math> requires a discretization of third-order accuracy to obtain a local error of <math display="inline">O(h^5)</math> because it is already multiplied by a factor <math display="inline">h^2</math>. Thus, using Eq. [[#theorem-eqn:HIFDIIM_4th|(1)]], we obtain

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.f_{xx}\right|_i = \left.\delta ^2f\right|_i- \left.bh^2f_{xxxx}\right|_i + C^4_i\{ f\} +O(h^3).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

On the other hand, the fourth-order derivative of <math display="inline">f</math> in <math display="inline">eh^4 f_{xxxx}</math> only requires an approximation of the first-order accuracy because its coefficient includes the term <math display="inline">h^4</math>; thus, we still have a local <math display="inline">O(h^5)</math> to keep a global sixth-order accurate method. Now applying Eqs. [[#theorem-eqn:HIFDIIM_4th|(1)]] and [[#theorem-eqn:HIFDIIM_2nd|(2)]] we get

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.f_{xxxx}\right|_i =\delta ^2\left.f_{xx}\right|_i + C^2_i\{ f_{xx}\} + O(h)</math>
|-
| style="text-align: center;" | <math> =\delta ^2\left(\left[\delta ^2f -bh^2\partial _{xx}f_{xx}\right|_i + C^4_i\{ f\} + O(h^3)\right)+ C^2_i\{ f_{xx}\} + O(h), </math>
|}
|}

then,

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.f_{xxxx}\right|_i =\left.\delta ^4f\right|_i + \delta ^2C^4_i\{ f\} + C^2_i\{ f_{xx}\} + O(h).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta ^4f = \dfrac{1}{h^4}\left(f_{i+2}-4f_{i+1}+6f_i-4f_{i-1}+f_{i-2}\right). </math>
|}
|}

Note that now the contribution notation includes from which function <math display="inline">f</math> or <math display="inline">f_{xx}</math> it comes. We also remark that term <math display="inline">bh^2\delta ^2f_{xxxx}</math> will be small if we guarantee <math display="inline">\delta ^2f_{xxxx}\approx f_{xxxxxx}</math> is bounded, which is always true in our analysis. Now, using identities [[#eq-28|(28)]] and [[#eq-29|(29)]], we obtain the discretization of the right-hand side as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_i + bh^2 f_{xx}|_i+eh^4 f_{xxxx}|_i=\,df_{i-2}+ (b-4d)f_{i-1}+\,(1-2b+6d)f_i+\,(b-4d)f_{i+1}+\,df_{i+2} </math>
|-
| style="text-align: center;" | <math> +\; bh^2C^4_i\{ f\} + dh^4\left(C^2_i\{ f_{xx}\} + \delta ^2C^4_i\{ f\} \right)+ O(h^5). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Finally, the HIFD-IIM for the 1D Poisson equation [[#eq-24|(24)]] at <math display="inline">x_i</math> is given by

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left.\delta ^2 u\right|_i = \,df_{i-2}+ (b-4d)f_{i-1}+\,(1-2b+6d)f_i+\,(b-4d)f_{i+1}+\,ef_{i+2} + \;\mathcal{C} +  O(h^5),   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

where <math display="inline">\mathcal{C}=(C_f)_i-C_i^6\{ u\} </math> and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(C_f)_i:=  bh^2C^4_i\{ f\} + dh^4C^2_i\{ f_{xx}\} + dh^2\left(C^4_{i+1}\{ f\} -2C^4_i\{ f\} +C^4_{i-1}\{ f\} \right), </math>
|}
|}

for <math display="inline">i=2, \dots , N</math>. Note that Dirichlet boundary conditions are directly applied at <math display="inline">i=1</math> and <math display="inline">i=N+1</math>.

Using definitions [[#eq-19|(19)]], [[#eq-20|(20)]], [[#eq-22|(22)]], and [[#eq-23|(23)]], we can simplify contribution <math display="inline">(C_f)_i</math> as

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(C_f)_i  =  \begin{cases}-d\left([f]+h_R\left[f_x\right]+\dfrac{h_R^2}{2}\left[f_{xx}\right]+\dfrac{h_R^3}{3!}\left[f_{xxx}\right]+\dfrac{h_R^2}{4!}\left[f_{xxxx}\right]\right)&i=I-1\\[2mm] -\left(a_0[f]+r_1\left[f_x\right]+r_2\left[f_{xx}\right]+r_3\left[f_{xxx}\right]+r_4\left[f_{xxxx} \right]\right)& i=I,  \\[2mm]    a_0[f]+l_1\left[f_x\right]+l_2\left[f_{xx}\right]+l_3\left[f_{xxx}\right]+l_4\left[f_{xxxx} \right],& i=I+1,\\[2mm]    d\left([f]+h_L\left[f_x\right]+\dfrac{h_L^2}{2}\left[f_{xx}\right]+\dfrac{h_L^3}{3!}\left[f_{xxx}\right]+\dfrac{h_L^4}{4!}\left[f_{xxxx} \right]\right)& i=I+2,\\[2mm]    0 & \mathrm{otherwise}, \end{cases}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

where <math display="inline">a_0 = b-d</math>, and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{ll}r_1 = bh_R + d(h_L-2h_R),  &  l_1 = bh_L + d(h_R-2h_L), \\[2mm] r_2 = \dfrac{1}{2}\left(bh^2_R + d(h_L^2-2h_R^2+2h^2)\right),  & l_2 = \dfrac{1}{2}\left(bh^2_L + d(h_R^2-2h_L^2+2h^2)\right),\\[2mm] r_3 = \dfrac{1}{3!}\left(bh^3_R + d(h_L^3-2h_R^3+6h^2h_R)\right),  & l_3 = \dfrac{1}{3!}\left(bh^3_L + d(h_R^3-2h_L^3+6h^2h_L)\right),\\[2mm] r_4 = \dfrac{1}{4!}\left(bh^4_R + d(h_L^4-2h_R^4+24h^2h_R^2)\right),  & l_4 = \dfrac{1}{4!}\left(bh^4_L + d(h_R^4-2h_L^4+24h^2h_L^2)\right). \end{array} </math>
|}
|}

===3.1 1D HIFD-IIM and the principal jump conditions===

Note that <math display="inline">\left[u\right]</math>, <math display="inline">\left[u_{x}\right]</math>, <math display="inline">\left[u_{xx}\right]</math>, <math display="inline">\left[u_{xxx}\right]</math>, <math display="inline">\left[u_{xxxx}\right]</math>, <math display="inline">\left[u_{xxxxx}\right]</math> and <math display="inline">\left[u_{xxxxxx}\right]</math> must be known to apply the proposed sixth-order HIFD-IIM. Thus, it seems that more jump conditions of <math display="inline">u</math> rather than the principal jump conditions are required to have a sixth-order accurate method. However, we can use the Poisson equation [[#eq-24|(24)]] to obtain relations between <math display="inline">u</math>, <math display="inline">f</math> and their derivatives as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[u_{xx}\right]=[f], \quad  \left[u_{xxx}\right]=\left[f_{x}\right], \quad  \left[u_{xxxx}\right]=\left[f_{xx}\right], \quad  \left[u_{xxxxx}\right]=\left[f_{xxx}\right], \quad  \left[u_{xxxxxx}\right]=\left[f_{xxxx}\right].  </math>
|}
|}

Thus, the total jump contribution for the 1D problem, <math display="inline">\mathcal{C}</math>, is given by 

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathcal{C}=  \begin{cases}  -d\left([f]+h_R\left[f_x\right]+\dfrac{h_R^2}{2}\left[f_{xx}\right]+\dfrac{h_R^3}{3!}\left[f_{xxx}\right]+\dfrac{h_R^2}{4!}\left[f_{xxxx}\right]\right), &i=I-1,\\[3mm]   \dfrac{1}{h^2}\left\{[u] + h_R\left[u_{x}\right]\right\}-\left(R_0[f]+R_1\left[f_x\right]+R_2\left[f_{xx}\right]+R_3\left[f_{xxx}\right]+R_4\left[f_{xxxx} \right]\right), & i=I,  \\[3mm] -\dfrac{1}{h^2}\left\{[u] + h_L\left[u_{x}\right]\right\}+L_0[f]+L_1\left[f_x\right]+L_2\left[f_{xx}\right]+L_3\left[f_{xxx}\right]+L_4\left[f_{xxxx} \right],& i=I+1,\\[3mm]    d\left([f]+h_L\left[f_x\right]+\dfrac{h_L^2}{2}\left[f_{xx}\right]+\dfrac{h_L^3}{3!}\left[f_{xxx}\right]+\dfrac{h_L^4}{4!}\left[f_{xxxx} \right]\right),& i=I+2,\\[3mm]    0, & \mathrm{otherwise}, \end{cases} </math>
|}
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>R_0 = a_0-\dfrac{h_R^2}{2h^2},  \quad R_k =  r_k-\dfrac{h_R^{k+2}}{(k+2)!h^2}, \quad L_0 = a_0-\dfrac{h_L^2}{2h^2}, \quad L_k = l_k-\dfrac{h_L^{k+2}}{(k+2)!h^2}. </math>
|}
|}

for <math display="inline">k=1,2,3,4</math>. Thus, the contribution <math display="inline">\mathcal{C}=(C_f)_i-C^6_i\{ u\} </math> depends only on the principal jump conditions and right-hand side jumps <math display="inline">[f]</math>, <math display="inline">\left[f_{x}\right]</math>, <math display="inline">\left[f_{xx}\right]</math>, <math display="inline">\left[f_{xxx}\right]</math>, and <math display="inline">\left[f_{xxxx} \right]</math>.

''Remark 3'': For the 1D Poisson problem, a second-order IIM (<math display="inline">b=0</math>, <math display="inline">d=0</math>) only requires knowing the principal jump conditions <math display="inline">[u]</math>, <math display="inline">\left[u_{x}\right]</math>, and <math display="inline">\left[f\right]</math>. On the other hand, a fourth-order IFD-IIM (<math display="inline">d=0</math>) is obtained knowing the principal jump conditions <math display="inline">[u]</math>, <math display="inline">\left[u_x\right]</math>, <math display="inline">\left[f \right]</math>, and additionally <math display="inline">\left[f_x \right]</math> and <math display="inline">\left[f_{xx} \right]</math>. However, the extra jumps are from the right-hand side, which is already known analytically or can be approximated using the current values of <math display="inline">f</math>. In this paper, we will assume that we know them.

''Remark 4'': We can achieve sixth or fourth-order approximation in some particular grids even if we do not know high-order derivative jumps. For example, for the sixth-order HIFD-IIM, if <math display="inline">h_R/h=1</math>, then <math display="inline">h_L=0</math>, and both weight terms next to the fourth-order derivative jump of <math display="inline">f</math> are equal to zero. Thus, we do not require to know the jump condition <math display="inline">\left[f_{xxxx}\right]</math> to obtain a sixth-order method when the interface is located at a grid point. Similarly, we can get a fourth-order scheme even if we do not know jump condition <math display="inline">\left[f_{xx}\right]</math> when the interface is located at a grid point.

==4. Two-dimensional problem==

We now apply the methodology to 2D problems and straight interfaces. We study the 2D Poisson problem given by

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>u_{xx}(x,y) +u_{yy}(x,y) = f(x,y),\quad (x,y)\in \Omega = (\mathfrak{a},x_\alpha ) \cup (x_\alpha , \mathfrak{b})\times (\mathfrak{c},\mathfrak{d}),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

where <math display="inline">u</math> and <math display="inline">f</math> can be discontinuous functions at <math display="inline">x=x_\alpha </math>, and the principal jump conditions are known functions in the variable <math display="inline">y</math> and specified as follows

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[u\right]= \mathfrak{p}(x_\alpha ,y),  \quad \mathrm{and} \quad  \left[u_{{n}} \right]= \mathfrak{q}(x_\alpha ,y) , \qquad \mathfrak{c} \leq y \leq \mathfrak{d}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

The numerical domain is discretized using an uniform mesh, <math display="inline">h = (\mathfrak{b}-\mathfrak{a})/N</math> (<math display="inline">N</math> sub-intervals in the <math display="inline">x</math>-direction) and assuming that <math display="inline">\mathfrak{d}-\mathfrak{c} = Mh</math> (<math display="inline">M</math> sub-intervals in the <math display="inline">y</math>-direction). We denote as <math display="inline">u_{ij}=u(x_i,y_j)</math> and <math display="inline">f_{ij}=f(x_i,y_j)</math> where the <math display="inline">ij</math>th point is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(x_i,y_j)=(\mathfrak{a}+(i-1)h,\mathfrak{c}+(j-1)h), \qquad i=1,\dots , N, N+1, j=1, \dots , M, M+1. </math>
|}
|}

The grid points are also classified into two types: regular and irregular. If the straight interface, <math display="inline">x = x_\alpha </math>, intersects the FD stencil surrounding the <math display="inline">ij</math>th point, the center point is called irregular; otherwise, the grid point is regular (Figure [[#img-3|3]]).

<div id='img-3'></div>
{| class="wikitable" style="margin: 0em auto 0.1em auto;border-collapse: collapse;width:60%;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Balam_338770048-Fig_IRPoints.png|420px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figure 3'''. Two-dimensional computational domain with a uniform mesh showing regular and irregular grid points. On the right-hand side, we list different types of stencils used to discretize the 2D Poisson problem
|}


Finally, let us define the discrete operators <math display="inline">\delta ^2_x</math> and <math display="inline">\delta ^4_x</math> at the <math display="inline">ij</math>th point as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta ^2_x u(x_i,y_j) = \dfrac{u_{i+1j}-2u_{ij}+u_{i-1j}}{h^2},</math>
|-
| style="text-align: center;" | <math> \delta ^4_x u(x_i,y_j) = \delta ^2_x\delta ^2_x u(x_i,y_j) = \dfrac{u_{i+2j}-4u_{i+1j}+6u_{ij}-4u_{i-1j}+u_{i-2j}}{h^4}. </math>
|}
|}

Formulas <math display="inline">\delta ^2_y</math> and <math display="inline">\delta ^4_y</math> are defined similarly.

We develop the discretization around the <math display="inline">ij</math>th point. However, to simplify our exposition we drop the evaluation at <math display="inline">ij</math>th point. We start the discretization by applying the HIFD operator [[#eq-7|(7)]] in <math display="inline">x</math>- and <math display="inline">y</math>-direction to the Poisson equation [[#eq-33|(33)]], as

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _y^4 u_{xx} + \partial _y^4 u_{yy} = \partial _y^4 f, \quad \mathrm{and} \quad  \partial _x^4 u_{xx} + \partial _x^4 u_{yy} = \partial _x^4 f, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

where

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _x^4(\cdot ):= (\cdot ) + bh^2(\cdot )_{xx} + eh^4(\cdot )_{xxxx}, \quad \textrm{ and } \quad \partial _y^4(\cdot ):= (\cdot ) + bh^2(\cdot )_{yy} + eh^4(\cdot )_{yyyy}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Adding both equations in Eq. [[#eq-35|(35)]],

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial ^4_x u_{xx} + \partial ^4_y u_{yy} + \partial ^4_y u_{xx} + \partial ^4_x u_{yy} =  \partial _x^4 f + \partial _y^4 f.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

In the following sections, we will describe the discretization of Eq. [[#eq-37|(37)]] for regular and irregular grid points.

===4.1 The HIFD-IIM at regular points===

Using one-dimensional formula [[#eq-6|(6)]], Eq. [[#eq-36|(36)]] and some algebraic simplifications, the sixth-order implicit method applied on the regular grid points is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right)+ (e-2b^2)h^4\left(\delta _x^4\delta _y^2u+\delta _y^4\delta _x^2u\right)</math>
|-
| style="text-align: center;" | <math>  =  f + b h^2 \left(f_{xx} + f_{yy}\right)+ eh^4\left(f_{xxxx} + f_{yyyy}\right)+ O(h^6). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

The implicit methods have not only high accuracy, but they are also more efficient in terms of the number of iterations required to solve the linear system of the discretization <span id='citeF-13'></span>[[#cite-13|[13]]]. In the next section we work with the irregular grind points.

===4.2 The HIFD-IIM at irregular points===

As the interface is a vertical line, the left-hand side terms without cross derivatives of [[#eq-37|(37)]] can be approximated using the HIFD-IIM formulation as follows

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial ^4_x u_{xx} = \delta _x^2u + C^6\{ u\} +O(h^5),\quad \textrm{ and } \quad  \partial ^4_y u_{yy} = \delta _y^2u +O(h^6). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

The left-hand side terms with cross derivatives of [[#eq-37|(37)]] require more attention. Using the Poisson equation and Theorem [[#theorem-thm:xj_right_left|1]], we can show that the following equation holds

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial ^4_y u_{xx}  + \partial ^4_x u_{yy}  =  f + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right)+ (e-2b^2)h^4\left(\delta _x^4\delta _y^2u+\delta _y^4\delta _x^2u\right)</math>
|-
| style="text-align: center;" | <math> +\;bh^2\delta _y^2C^6\{ u\} +bh^2\left(C^4\{ u_{yy}\} + bh^2C^2\{ u_{yyyy}\} \right)</math>
|-
| style="text-align: center;" | <math> +\;(e-2b^2)h^4\left(\delta _x^2C^4\{ u_{yy}\} + bh^2\delta _x^2C^2\{ u_{yyyy}\} + C^2\{ u_{yyxx}\}  + \delta _y^4C^6\{ u\} \right)+ O(h^5).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

On the other hand, using the definition of operator <math display="inline">\partial ^4</math> given by [[#eq-36|(36)]], the right-hand side can be written as follows

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _x^4 f + \partial _y^4 f = 2f + b h^2 \left(f_{xx} + f_{yy}\right)+ eh^4\left(f_{xxxx} + f_{yyyy}\right).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

Substituting equations [[#eq-39|(39)]]-[[#eq-41|(41)]] into [[#eq-37|(37)]] we obtain

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right)+ (e-2b^2)h^4\left(\delta _x^4\delta _y^2u+\delta _y^4\delta _x^2u\right)</math>
|-
| style="text-align: center;" | <math> =  f + b h^2 \left(f_{xx} + f_{yy}\right)+ eh^4\left(f_{xxxx} + f_{yyyy}\right) +\;\mathcal{C}_u + O(h^5),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{C}_u  = -C^6\{ u\} -bh^2\delta _y^2C^6\{ u\} -bh^2\left(C^4\{ u_{yy}\} + bh^2C^2\{ u_{yyyy}\} \right)</math>
|-
| style="text-align: center;" | <math> -(e-2b^2)h^4\left(\delta _x^2C^4\{ u_{yy}\} + bh^2\delta _x^2C^2\{ u_{yyyy}\} + C^2\{ u_{yyxx}\}  + \delta _y^4C^6\{ u\} \right) + O(h^5). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

All above terms are evaluated at the <math display="inline">ij</math>th point (omitted to avoid misunderstandings). If we known explicitly <math display="inline">f</math>, then the discretization is complete. However, there are problems which the Poisson implementation only allows to known the right-hand side at the discretization points. In the following section, we focus on this case.

===4.3 Right-hand side approximation===

It is possible that <math display="inline">f</math> is only know at the grid points. If this is the case, we only require to compute terms <math display="inline">f_{xx}</math> and <math display="inline">f_{xxxx}</math> using HIFD-IIM because the interface is a vertical line, see Fig.&nbsp;[[#img-3|3]]. Moreover, formulas [[#eq-28|(28)]] and [[#eq-29|(29)]] are valid in this case. Then,

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f + b h^2 \left(f_{xx} + f_{yy}\right)+ eh^4\left(f_{xxxx} + f_{yyyy}\right) =f + bh^2\left(\delta _x^2f + \delta _y^2f\right)+ dh^4\delta _x^4f + \; eh^4\delta _y^4f</math>
|-
| style="text-align: center;" | <math> +\;bh^2C^4\{ f\} +dh^4\left(\delta _x^2C^4\{ f\} + C^2\{ f_{xx}\} \right)+O(h^5),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where <math display="inline">d = e-b^2</math>. Finally, putting together equations [[#eq-42|(42)]] and [[#eq-44|(44)]], we obtain the ''fully discretized 2D problem'' as follows

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right)+ (e-2b^2)h^4\left(\delta _x^4\delta _y^2u+\delta _y^4\delta _x^2u\right)</math>
|-
| style="text-align: center;" | <math> = f + bh^2\left(\delta _x^2f + \delta _y^2f\right)+ dh^4\delta _x^4f+ eh^4\delta _y^4f +\;\mathcal{C}_f +\mathcal{C}_u  +O(h^5),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

where <math display="inline">\mathcal{C}_f = bh^2C^4\{ f\} + dh^4\left(\delta _x^2C^4\{ f\} + C^2\{ f_{xx}\} \right)</math>.

If the 2D Poisson problem admits a smooth solution, then all contribution terms in [[#eq-45|(45)]] are equal to zero, and a sixth-order method is recovered. In the presence of an interface, the contributions <math display="inline">C^2</math>, <math display="inline">C^4</math>, and <math display="inline">C^6</math> on regular grid points are always zero. This scheme requires a stencil type C, see Fig. [[#img-3|3]]

<span id='theorem-eqn:2dPoissonfullDscrt4th'></span>''Remark 5'': Taking <math display="inline">e=0</math> into equations [[#eq-36|(36)]] and only considering approximations up to four-order accurate, we obtain the IFD-IIM

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right) =  f + bh^2\left(\delta _x^2f + \delta _y^2f\right)+ \widehat{\mathcal{C}_f} + \widehat{\mathcal{C}_u}+O(h^3),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

where <math display="inline">\widehat{\mathcal{C}_f}=C^4\{ f\} </math> and <math display="inline">\widehat{\mathcal{C}_u} = -C^4\{ u\} -bh^2\delta _y^2C^4\{ u\} -bh^2C^2\{ u_{yy}\} </math>. In addition, if the 2D Poisson problem admits a smooth solution, all contribution terms in [[#theorem-eqn:2dPoissonfullDscrt4th|(5)]] vanish, and we get a fourth-order implcit scheme. In fact, the contributions <math display="inline">C^k</math> are always zero on regular grid points. This scheme requires a stencil type A, see Fig.&nbsp;[[#img-3|3]]. Furthermore, if we consider that <math display="inline">b=0</math> in [[#theorem-eqn:2dPoissonfullDscrt4th|(5)]] then we obtain the standard explicit immersed interface method, called IIM, which is second order accurate.

''Remark 6'': It is possible to develop a discrete formulation for the 2D Poisson problem when the straight interface is horizontal. The resulting scheme is

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right)+ (e-2b^2)h^4\left(\delta _x^4\delta _y^2u+\delta _y^4\delta _x^2u\right)</math>
|-
| style="text-align: center;" | <math> =f + bh^2\left(\delta _x^2f + \delta _y^2f\right)+ dh^4\delta _x^4f+ eh^4\delta _y^4f+\;\mathcal{C} +O(h^5), </math>
|}
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{C} = bh^2C^4\{ f\} + dh^4\left(\delta _y^2C^4\{ f\} + C^2\{ f_{yy}\} \right)</math>
|-
| style="text-align: center;" | <math> -\;C^6\{ u\} -bh^2\delta _x^2C^6\{ u\} -bh^2\left(C^4\{ u_{xx}\} + bh^2C^2\{ u_{xxxx}\} \right)</math>
|-
| style="text-align: center;" | <math> -\;(e-2b^2)h^4\left(\delta _y^2C^4\{ u_{xx}\} + bh^2\delta _x^2C^2\{ u_{xxxx}\} + C^2\{ u_{xxyy}\}  + \delta _x^4C^6\{ u\} \right). </math>
|}
|}

The constants <math display="inline">C^2</math>, <math display="inline">C^4</math>, and <math display="inline">C^6</math> are computed using jump derivatives in the <math display="inline">y</math>-direction.

===4.4 Boundary treatment===

The grid points close to the interface require special treatment because the HIFD-IIM formulation of the 21-point stencil C cannot be applied there. For regular grid points, we use the following fourth-order method

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u + b h^2 \left(\delta _y^2\delta _x^2u+\delta _x^2\delta _y^2u\right) =  f + bh^2\left(\delta _x^2f + \delta _y^2f\right)+ O(h^4). </math>
|}
|}

The deduction of this scheme, stencil type A, can be found in Zapata and Balam <span id='citeF-13'></span>[[#cite-13|[13]]].

It is necessary to develop a new fourth-order method for the grid points near both boundary and interface. Following the same ideas to develop the sixth-order method, we obtain a new scheme with stencil of type B

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _x^2u+ \delta _y^2u  + b h^2\left(\delta _x^2\delta _y^2u + \delta _y^2\delta _x^2u\right)+ (e-2b^2)h^4\delta _x^4\delta _y^2u  =  f + bh^2\left(\delta _x^2f + \delta _y^2f\right)+ h^4\delta _x^4f </math>
|-
| style="text-align: center;" | <math> -\;b^2h^4\delta _y^4f+\mathcal{C}_{boundary} + O(h^4), </math>
|}
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{C}_{boundary} =  bh^2C^4\{ f\} + dh^4\left(\delta _x^2C^4\{ f\} + C^2\{ f_{xx}\} \right) -C^6\{ u\} -bh^2\delta _y^2C^6\{ u\} </math>
|-
| style="text-align: center;" | <math> -\;bh^2\left(C^4\{ u_{yy}\} + bh^2C^2\{ u_{yyyy}\} \right)</math>
|-
| style="text-align: center;" | <math> -\;(e-2b^2)h^4\left(\delta _x^2C^4\{ u_{yy}\} + bh^2\delta _x^2C^2\{ u_{yyyy}\} + C^2\{ u_{yyxx}\} \right)+ O(h^4). </math>
|}
|}

''Remark 7'': We emphasized that the proposed methods are high-order accuracy regardless of the position of the interface concerning the grid. Moreover, the scheme does not assume restrictions in the jumps, such as the natural jump conditions (<math display="inline">[u_x]=0</math>, <math display="inline">[u_y]=0</math>). This characteristic is a significant advantage of the proposed HIFD-IIM, besides the higher-order, compared with the fourth-order simplified immersed interface method developed by Feng et al.<span id='citeF-39'></span>[[#cite-39|[39]]].

==5. Numerical results==

This section tests the HIFD-IIM for different one- and two-dimensional examples with straight interfaces. In the following simulations, we numerically solve the Poisson equation for a given right-hand side function and compare it with its analytic solution. We present different examples to test the HIFD-IIM capabilities. First, we investigate the method's accuracy for one-dimensional problems. Next, we validate the HIFD-IIM for two-dimensional solutions with straight interfaces.

The errors are reported utilizing the <math display="inline">L_\infty </math>-norm, as <math display="inline">\left\|e \right\|_\infty = \max _{i,j} \left|u_{ij}-U_{ij} \right|</math>, where <math display="inline">u_{ij}</math> and <math display="inline">U_{ij}</math> corresponds to the exact and numerical  solution at <math display="inline">(x_i,y_j)</math>, respectively. The estimated order of accuracy is computed as

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\textrm{Order} := \log \left(\left\|e \right\|_{N_1}/ \left\|e \right\|_{N_2}\right)/\log \left(N_2/N_1 \right), </math>
|}
|}

where <math display="inline">N_1</math> and <math display="inline">N_2</math> indicates the different number of sub-intervals. In all tables, the last row shows the numerical order calculated by the regression-line slope based on a Least Squares Method (LSM).

===5.1 One-dimensional examples===

We initially consider the 1D problem [[#eq-24|(24)]] to analyze the order of the proposed implicit methods. The numerical method is tested using three different tests. Example 1.A was designed to verify the high-order implicit method for smooth solutions. Example 1.B studies a Poisson equation with a discontinuous solution in a single interface point. Example 1.C presents a discontinuous problem with multiple interface points. Thus, the proposed solution is taken from the following list of functions:

<span id="eq-47"></span>
<span id="eq-48"></span>
<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\textrm{Example 1.A:}  \quad u(x) = e^{-4\pi (x-1/4)^2}.  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|-
| style="text-align: center;" | <math> \textrm{Example 1.B:}  \quad u(x) = \begin{cases}\sin (\pi x)& x\leq x_\alpha ,\\ \cos (\pi x)& x   > x_\alpha{.} \end{cases}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|-
| style="text-align: center;" | <math> \textrm{Example 1.C:}  \quad  u(x) = \begin{cases}\sin (2\pi x)& x\leq x_{\alpha _1},\\ \cos (2\pi x)& x_{\alpha _1}<x \leq x_{\alpha _2}, \\ \sin (2\pi x)& x_{\alpha _2} < x. \end{cases}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}
|}

where <math display="inline">x_\alpha </math>, <math display="inline">x_{\alpha _1}</math> and <math display="inline">x_{\alpha _2}</math> are known values corresponding to the interface location. The right-hand side function, <math display="inline">f</math>, is obtained directly from [[#eq-33|(33)]] and <math display="inline">u</math>. For all cases, we impose the Dirichlet boundary conditions according to <math display="inline">u</math>. The computational domain is the interval <math display="inline">[0,1]</math>, and the grid spacing is <math display="inline">h = 1 /N</math> for different <math display="inline">N</math> numbers.

For Example 1.A, due to the solution's regularity, the jump contributions are equal to zero. Table [[#table-1|1]] presents the convergence analysis of Example 1.A for different grid resolutions. Note that the IIM, IFD-IIM, and HIFD-IIM achieve their corresponding order of accuracy. These results match with the ones obtained using an implicit methodology as presented in <span id='citeF-13'></span>[[#cite-13|[13]]].


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>'''Table. 1''' Convergence analysis of Example 1.A for IFD-IIM and HIFD-IIM.
|-  
| style="text-align: right;" |      
|       
|    
|-
| <math display="inline">N</math>  
| <math display="inline">L_\infty </math>-norm 
|  Order 
| <math display="inline">L_\infty </math>-norm 
|  Order
| <math display="inline">L_\infty </math>-norm 
|  Order
|-
| 10 
|  7.52e-02 
|  &#8211;-   
|  9.30e-03 
|  &#8211;-  
|  1.39e-04 
|  &#8211;-   
|-
| 20 
|  1.70e-02 
|  2.15  
|  5.05e-04 
|  4.20 
|  1.52e-06 
|  6.51  
|-
| 40 
|  4.15e-03 
|  2.03  
|  3.06e-05 
|  4.04 
|  1.77e-08 
|  6.43  
|-
| 80 
|  1.03e-03 
|  2.01  
|  1.90e-06 
|  4.01 
|  2.34e-10 
|  6.24  
|-
| 160 
|  2.57e-04 
|  2.00  
|  1.18e-07 
|  4.00 
|  3.24e-12 
|  6.18  
|-
| LSM 
|    
| '''2.04''' 
| 
|   '''4.06''' 
| 
|   '''6.34''' 
|-
|
|}
Example 1.B shows the capacity of the proposed method to solve a single interface problem located at <math display="inline">x = x_\alpha </math>. We test two different interface points: <math display="inline">x_\alpha=0.40</math> and <math display="inline">x_\alpha=0.63</math>. We initially select the mesh grid given by <math display="inline">N=10\times{2}^{n}</math>, thus the first interface is always located on a grid point (<math display="inline">h_R/h=1</math>). For the second case, we have different <math display="inline">h_R/h</math> values for the same <math display="inline">N</math> numbers. Fig.&nbsp;[[#img-4|4]] shows the numerical and exact solution  using <math display="inline">N = 40</math>. As expected, the exact solution is accurately recovered for both cases.

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa1B_1D_sol.png|480px|Numerical and exact solution of Example 1.B using N = 40 using (a) x_α= 0.4, and (b) x_α= 0.63.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Numerical and exact solution of Example 1.B using <math>N = 40</math> using (a) <math>x_\alpha = 0.4</math>, and (b) <math>x_\alpha = 0.63</math>.
|}
Table [[#table-2|2]] shows the convergence analysis for Example 1.B. As expected, the desired order of accuracy are obtained for the two <math display="inline">x_\alpha </math> values. Observe that high-order methods do not depend on the location of the interface. However, their error magnitude presents minor variations due to the interface position. Errors with mesh size close to <math display="inline">N = 160</math> have a random behavior due the effect of arithmetic operations close to the machine precision. Fig.&nbsp;[[#img-5|5]] shows the error analysis corresponding to interface locations <math display="inline">x_\alpha = 0.40</math> and <math display="inline">x_\alpha = 0.63</math> for <math display="inline">N = 10,11,12, \dots ,100</math>. Note that errors of IFD-IIM give a good behavior even if the <math display="inline">h_R/h</math> varies and they are close to fourth order. As expected, errors of HIFD-IIM are also near sixth order.


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>'''Table. 2''' Convergence analysis of Example 1.B using the IFD-IIM and HIFD-IIM.
|-     
|       
|    
|-
| <math display="inline">N</math>  
| <math display="inline">L_\infty </math>-norm 
|  Order 
| <math display="inline">L_\infty </math>-norm 
|  Order
| <math display="inline">L_\infty </math>-norm 
|  Order
|-
| 10 
|  1.69e-02 
|  &#8211;-  
|  5.75e-05 
|  &#8211;-  
|  6.42e-07 
|  &#8211;-  
|-
| 20 
|  4.21e-03 
|  2.00 
|  3.58e-06 
|  4.00 
|  1.01e-08 
|  6.42 
|-
| 40 
|  1.05e-03 
|  2.00 
|  2.24e-07 
|  4.00 
|  1.59e-10 
|  6.21 
|-
| 80 
|  2.63e-04 
|  2.00 
|  1.40e-08 
|  4.00 
|  2.48e-12 
|  6.11 
|-
| 160 
|  6.57e-05 
|  2.00 
|  8.74e-10 
|  4.00 
|  4.55e-14 
|  5.82 
|-
| LSM 
| 
|  '''2.00''' 
| 
|  '''4.00''' 
| 
|  '''5.95'''
|-
| 
|      
|       
|    
|-
| <math display="inline">N</math>  
| <math display="inline">L_\infty </math>-norm 
|  Order 
| <math display="inline">L_\infty </math>-norm 
|  Order
| <math display="inline">L_\infty </math>-norm 
|  Order
|-
| 10 
|  6.63e-03 
|  &#8211;-  
|  3.42e-05 
|  &#8211;-  
|  6.38e-07 
|  &#8211;-   
|-
| 20 
|  1.48e-03 
|  2.16 
|  1.97e-06 
|  4.12 
|  9.89e-09 
|  6.44  
|-
| 40 
|  4.46e-04 
|  1.73 
|  1.39e-07 
|  3.82 
|  1.58e-10 
|  6.18  
|-
| 80 
|  9.37e-05 
|  2.25 
|  7.75e-09 
|  4.16 
|  2.40e-12 
|  6.15  
|-
| 160 
|  2.75e-05 
|  1.77 
|  5.37e-10 
|  3.85 
|  7.31e-14 
|  5.08  
|-
| LSM 
| 
|  '''1.98''' 
| 
|  '''3.99''' 
| 
|  '''5.81'''
|-
|
|}
<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa1B_1D_CloudOrder4.png|480px|]]
|[[Image:Draft_Balam_338770048-Fig_Exa1B_1D_CloudOrder6.png|480px|Convergence analysis of Example 1.B for N = 10,\dots ,100 using (a),(b) IFD-IIM  and (c),(d) HIFD-IIM for x_α= 0.40 and x_α= 0.63.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Convergence analysis of Example 1.B for <math>N = 10,\dots ,100</math> using (a),(b) IFD-IIM  and (c),(d) HIFD-IIM for <math>x_\alpha = 0.40</math> and <math>x_\alpha = 0.63</math>.
|}
We remark that, the contribution formula includes jumps <math display="inline">[u]</math>, <math display="inline">[u_x]</math>, <math display="inline">[u_{xx}]</math>, <math display="inline">[u_{xxx}]</math>, and <math display="inline">[u_{xxxx}]</math> to obtain a fourth-order accurate method. Fig.&nbsp;[[#img-6|6]] shows that if we add additional jumps of high-order derivatives to <math display="inline">\mathcal{C}</math>, such as <math display="inline">[u_{xxxxx}]</math>, we observe that the error oscillation decreases compared to Fig.&nbsp;[[#img-5|5]] results. It is expected because now the method is <math display="inline">O(h^4)</math> for the whole computational domain, including the irregular points. Thus, we can mitigate error oscillations due to interface position by adding high-order jumps. A similar behavior is observed for the sixth-order HIFD-IIM if we include the seventh derivative jump <math display="inline">[u_{xxxxxxx}]</math>.

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa1B_1D_OrderCloud5th.png|480px|Convergence analysis of Example 1.B for N = 10,\dots ,100 using (a) x_α= 0.40, and (b) x_α= 0.63. The contribution term includes jumps up to fifth-order ([uₓₓₓₓₓ]=[fₓₓₓ]).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Convergence analysis of Example 1.B for <math>N = 10,\dots ,100</math> using (a) <math>x_\alpha = 0.40</math>, and (b) <math>x_\alpha = 0.63</math>. The contribution term includes jumps up to fifth-order (<math>[u_{xxxxx}]=[f_{xxx}]</math>).
|}
Finally, Example 1.C investigates the method's capacity to solve a multiple interface problem. We only focus on two interface points located at <math display="inline">x_{\alpha _1}= 0.30</math> and <math display="inline">x_{\alpha _2} = 0.70</math>. However, the methodology could be applied for several interfaces by doing minor modifications in the implementation. Fig.&nbsp;[[#img-7|7]] presents the analytical and numerical solution using <math display="inline">N=40</math>. This figure also shows the corresponding absolute error. Fig.&nbsp;[[#img-8|8]] shows the error analysis for high-order IIM using different grid resolutions <math display="inline">N = 10, 20, 80 ,160</math>. As expected, the IFD-IIM and HIFD-IIM are fourth- and sixth-order accurate methods, respectively.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa1C_1D_SolError6th.png|480px|(a) Numerical and exact solution of Example 1.C with multiple interfaces using N = 40; (b) Absolute error of the numerical solution using the HIFD-IIM.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' (a) Numerical and exact solution of Example 1.C with multiple interfaces using <math>N = 40</math>; (b) Absolute error of the numerical solution using the HIFD-IIM.
|}
<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa1C_1D_Cloud4th6th.png|480px|Convergence error analysis of (c) the IFD-IIM, and (d) HIFDM-IIM, using different grid resolutions N = 10, 20, 40, 80 ,160.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Convergence error analysis of (c) the IFD-IIM, and (d) HIFDM-IIM, using different grid resolutions <math>N = 10, 20, 40, 80 ,160</math>.
|}

===5.2 Two-dimensional examples===

In this section we study the method's capacities to solve two-dimensional problems with different jump-contribution characteristics. Here we use straight interfaces located at <math display="inline">x = x_\alpha </math>. First, we consider smooth problem to verify the accurate implementation of the implicit method. Next, we analyze several discontinuous problems and finally, we include a more complex test where the jump derivatives increase rapidly.

For all 2D examples, the right-hand side function and jump conditions are computed from the corresponding exact solution, the computational domain is <math display="inline">\Omega = [-1,1]\times [-1,1]</math> and the grid size in both directions is the same using <math display="inline">N=10</math>, <math display="inline">20</math>, <math display="inline">40</math>, <math display="inline">80</math>, <math display="inline">160</math>. We impose Dirichlet boundary conditions. This paper uses the Successive-Over Relaxation (SOR) method with <math display="inline">\epsilon=10^{-14}</math> and <math display="inline">\omega=1.9</math> as tolerance and relaxation parameters, respectively.

====5.2.1 2D Poisson equation with smooth solution====

In the first 2D example, we solve the Poisson equation [[#eq-33|(33)]] for a smooth solution to show the correct implementation of the high-order implicit methods. In this case, the exact solution is given by

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\textrm{Example 2.A:} \qquad u(x,y) = \sin (3x)\sin (3y).  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

Table [[#table-3|3]] shows the convergence analysis of Example 2.A using different grid resolutions for the fourth- and sixth-order implicit formulation. As expected for a smooth solution, the implicit formulation improves the precision of the standard second-order numerical solution.


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>'''Table. 3''' Convergence analysis of Example 2.A.
|- 
| style="text-align: right;" |         
|        
|  
|-
| <math display="inline">N</math> 
| <math display="inline">L_\infty </math>-norm   
|  Order       
| <math display="inline">L_\infty </math>-norm   
|  Order 
| <math display="inline">L_\infty </math>-norm   
|  Order
|-
| 10 
|  2.65e-02 
|  &#8211;-  
|  3.07e-04 
|  &#8211;-  
|  4.64e-05 
|  &#8211;-  
|-
| 20 
|  7.06e-03 
|  1.91 
|  2.10e-05 
|  3.87 
|  7.54e-07 
|  5.94 
|-
| 40 
|  1.76e-03 
|  2.00 
|  1.32e-06 
|  3.99 
|  1.17e-08 
|  6.01 
|-
| 80 
|  4.40e-04 
|  2.00 
|  8.24e-08 
|  4.00 
|  1.82e-10 
|  6.00
|-
| 160 
|  1.10e-04 
|  2.00 
|  5.16e-09 
|  4.00 
|  2.95e-12 
|  5.95
|-
| LSM 
|  
| '''1.98''' 
|  
| '''3.97''' 
|  
| '''5.98''' 
|-
|
|}

====5.2.2 2D Poisson equation with straight interfaces====

In this section, we solve the 2D Poisson problem using the following set of functions:

<span id="eq-51"></span>
<span id="eq-51"></span>
<span id="eq-51"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\textrm{Example 2.B:}  \quad u(x,y) = \begin{cases}\sin (x)e^y & x\leq x_\alpha ,\\ -x^2e^y + 2& x_\alpha{<}x,\\ \end{cases}    \qquad \quad \quad \quad \textrm{where}    x_\alpha = 0.  </math>
|-
| style="text-align: center;" | <math> \textrm{Example 2.C:}  \quad u(x,y) = \begin{cases}\sin (x)\sin (y) & x\leq x_\alpha ,\\ e^x\sin (3y) + 1 & x_\alpha{<}x,\\ \end{cases}   \quad \quad \quad  \textrm{where}   x_\alpha = 0.4.  </math>
|-
| style="text-align: center;" | <math> \textrm{Example 2.D:}  \quad  u(x,y) = \begin{cases}\cos (3x)\cos (3y) & x\leq x_\alpha ,\\ \sin (3y)\sin (3y) & x_\alpha{<}x,\\ \end{cases}   \qquad \quad \textrm{where}   x_\alpha = 0.5.  </math>
|}
|}

Examples 2.B, 2.C and 2.D investigate the influence on the absolute error and the accuracy over several assumptions of the jump derivatives. Example 2.B analyzes the solution where jump derivatives in the <math display="inline">y</math>-direction vanish at <math display="inline">x_\alpha </math>. In Example 2.C, the jump derivative in the <math display="inline">y</math>-direction changes slower for <math display="inline">x<x_\alpha </math> than the one for <math display="inline">x>x_\alpha </math>. Example 2.D studies the problem without any assumption about the jump derivatives. Note that the interface location <math display="inline">x_\alpha </math> is also different for each example.

Fig.&nbsp;[[#img-9|9]] shows the numerical solution of these examples using a grid resolution of <math display="inline">N = 80</math>. As expected, the HIFD-IIM solves the problem accurately for each case.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa2_2D_ABC_solB.png|594px|Numerical solution of 2D Poisson equation with a straight interface corresponding to Examples 2.B, 2.C, and 2.D using N=80.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Numerical solution of 2D Poisson equation with a straight interface corresponding to Examples 2.B, 2.C, and 2.D using <math>N=80</math>.
|}
Table [[#table-4|4]] shows the convergence analysis for Example 3 using the implicit methods. As expected, the IFD-IIM and  HIFD-IIM schemes are close to fourth- and sixth-order, respectively. This table also confirms the standard IIM is second-order accurate. Observe that the order corresponding to the HIFD-IIM with <math display="inline">N=160</math> in Example 2.B is reduced due to the arithmetic operations close to the machine precision. Fig.&nbsp;[[#img-10|10]] shows more details of the convergence analysis using a cloud of points from <math display="inline">N = 10</math> to <math display="inline">N=160</math>.


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-4'></span>'''Table. 4''' Convergence analysis of Examples 2.B, 2.C, and 2.D using a straight interface.
|-
| 
| style="text-align: right;" |         
|        
|  
|-
| <math display="inline">N</math> 
| <math display="inline">L_\infty </math>-norm   
|  Order       
| <math display="inline">L_\infty </math>-norm   
|  Order 
| <math display="inline">L_\infty </math>-norm   
|  Order
|-
| 10 
|  5.76e-04 
|  &#8211;-  
|  5.06e-06 
|  &#8211;-  
|  2.10e-06 
|  &#8211;-   
|-
| 20 
|  1.50e-04 
|  1.94 
|  3.19e-07 
|  3.99 
|  3.78e-08 
|  5.79  
|-
| 40 
|  3.77e-05 
|  1.99 
|  2.00e-08 
|  3.99 
|  6.34e-10 
|  5.90  
|-
| 80 
|  9.45e-06 
|  2.00 
|  1.25e-09 
|  4.00 
|  1.03e-11 
|  5.95  
|-
| 160 
|  2.36e-06 
|  2.00 
|  7.84e-11 
|  4.00 
|  1.17e-12 
|  3.13  
|-
| LSM 
|  
| '''1.99''' 
|  
| '''4.00''' 
|  
| '''5.34''' 
|-
| 
|         
|        
|  
|-
| <math display="inline">N</math> 
| <math display="inline">L_\infty </math>-norm   
|  Order       
| <math display="inline">L_\infty </math>-norm   
|  Order 
| <math display="inline">L_\infty </math>-norm   
|  Order
|-
| 10 
|  2.65e-02 
|  &#8211;-  
|  8.40e-04 
|  &#8211;-  
|  3.01e-04 
|  &#8211;-  
|-
| 20 
|  6.60e-03 
|  2.00 
|  6.58e-05 
|  3.67 
|  6.74e-06 
|  5.48  
|-
| 40 
|  1.61e-03 
|  2.04 
|  4.52e-06 
|  3.86 
|  1.25e-07 
|  5.75  
|-
| 80 
|  3.98e-04 
|  2.01 
|  2.96e-07 
|  3.93 
|  2.14e-09 
|  5.87  
|-
| 160 
|  9.90e-05 
|  2.01 
|  1.90e-08 
|  3.97 
|  3.50e-11 
|  5.93  
|-
| LSM 
|  
| '''2.02''' 
|  
| '''3.87''' 
|  
| '''5.77''' 
|-
| 
|         
|        
|  
|-
| <math display="inline">N</math> 
| <math display="inline">L_\infty </math>-norm   
|  Order       
| <math display="inline">L_\infty </math>-norm   
|  Order 
| <math display="inline">L_\infty </math>-norm   
|  Order
|-
| 10 
|  4.62e-02 
|  &#8211;-  
|  5.46e-04 
|  &#8211;-  
|  2.62e-04 
|  &#8211;-   
|-
| 20 
|  1.13e-02 
|  2.03 
|  3.94e-05 
|  3.79 
|  2.91e-06 
|  6.49  
|-
| 40 
|  2.81e-03 
|  2.01 
|  2.63e-06 
|  3.91 
|  5.40e-08 
|  5.75  
|-
| 80 
|  7.04e-04 
|  2.00 
|  1.69e-07 
|  3.96 
|  9.12e-10 
|  5.89  
|-
| 160 
|  1.76e-04 
|  2.00 
|  1.07e-08 
|  3.98 
|  1.48e-11 
|  5.95  
|-
| LSM 
|  
| '''2.01''' 
|  
| '''3.91''' 
|  
| '''5.98''' 
|-
|
|}
<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa2_2D_ABC_cloud4thB.png|588px|]]
|[[Image:Draft_Balam_338770048-Fig_Exa2_2D_ABC_cloud6thB.png|588px|Convergence analysis of Examples 2.B, 2.C, and 2.D using the (a)-(c) IFD-IIM and (d)-(f) HIFD-IIM from N=10 to 100.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 10:''' Convergence analysis of Examples 2.B, 2.C, and 2.D using the (a)-(c) IFD-IIM and (d)-(f) HIFD-IIM from <math>N=10</math> to <math>100</math>.
|}

====5.2.3 2D Poisson solution with large jump derivatives ====

Finally, we construct the following example to analyze the variations in the errors due to jump derivative magnitudes. Thus, we have

<span id="eq-51"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\textrm{Example 2.E:}  u(x,y) = \begin{cases}0 & x\leq x_\alpha ,\\ 2+\ln (x(1+y^2)) & x_\alpha{<}x,\\ \end{cases} \quad  \textrm{where}   x_\alpha = 0.5.   </math>
|}
|}

Here, the jump derivatives for all points at the interface are <math display="inline">[u_x] = 2</math>, <math display="inline">[u_{xx}] = -4</math>, <math display="inline">[u_{xxx}] =16</math>, <math display="inline">[u_{xxxx}] = -96</math>, <math display="inline">[u_{xxxxx}] =  768</math>, and <math display="inline">[u_{xxxxxx}] = -7680</math>. Is important to remark that opposite to the previous examples, the jump derivatives increase rapidly. This behavior makes the problem challenging to solve.

Table [[#table-5|5]] and Fig.[[#img-11|11]] show the convergence analysis for Example 2.E. Numerical results show that Example 2.E has more variability in error than previous examples. However, the order of each technique is close to the proposed one. These findings also confirm that local truncation error depends not only on <math display="inline">h_R</math> and <math display="inline">h_L</math> but also on the jump magnitudes.


{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-5'></span>'''Table. 5''' Convergence analysis of Example 2.E.
|-      
| style="text-align: right;" |         
|        
|  
|-
| <math display="inline">N</math> 
| <math display="inline">L_\infty </math>-norm   
|  Order       
| <math display="inline">L_\infty </math>-norm   
|  Order 
| <math display="inline">L_\infty </math>-norm   
|  Order
|-
| 10 
|  1.11e-02 
|  &#8211;-  
|  2.72e-04 
|  &#8211;-  
|  1.38e-04 
|  &#8211;-   
|-
| 20 
|  1.80e-03 
|  2.62 
|  5.61e-05 
|  2.28 
|  2.48e-06 
|  5.80  
|-
| 40 
|  5.64e-04 
|  1.67 
|  3.62e-06 
|  3.96 
|  9.43e-08 
|  4.72  
|-
| 80 
|  1.57e-04 
|  1.85 
|  2.29e-07 
|  3.98 
|  2.61e-09 
|  5.17  
|-
| 160 
|  4.11e-05 
|  1.93 
|  1.43e-08 
|  3.99 
|  6.10e-11 
|  5.42  
|-
| LSM 
|  
| '''1.97''' 
|  
| '''3.64''' 
|  
| '''5.21''' 
|-
|
|}
<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Balam_338770048-Fig_Exa4_2D_ABC_solcloud4th6thB.png|588px|Numerical and exact solution of Example 2.E using N = 80; (b),(c) convergence error analysis of IFD-IIM and HIFD-IIM respectively, using different grid resolutions from N=10 to 100.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Numerical and exact solution of Example 2.E using <math>N = 80</math>; (b),(c) convergence error analysis of IFD-IIM and HIFD-IIM respectively, using different grid resolutions from <math>N=10</math> to <math>100</math>.
|}

==6 Conclusions==

The present paper introduces a new sixth-order immersed interface combined with an implicit finite difference to solve 2D Poisson problems with straight interfaces. The resulting numerical method is <math display="inline">O(h^6)</math> at regular points, and <math display="inline">O(h^5)</math> at irregular points. Furthermore, a fourth-order immersed interface method is obtained as a particular case of the proposed scheme. This paper also presents a numerical technique to handle the boundaries in the Poisson problem. The global accuracy of the sixth-order was demonstrated using several numerical examples. As expected, this approach does not depend on the interface position. For future work, the proposed approximation will be used to solve more general elliptic equations and interface shapes, and time-dependent problems in higher dimensions.

==References==

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''    H. J. Diersch, Fletcher, C.A.J. ''Computational Techniques for Fluid Dynamics. Vol. I: Fundamental and General Techniques. Vol. II: Specific Techniques for Different Flow Categories''. Springer-Verlag, 1998. <div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''    E. Javierre, C. Vuik, F.J. Vermolen, S. Van der Zwaag, A comparison of numerical models for one-dimensional Stefan problems, J. Comput. Appl. Math. 192 (2) (2006) 445&#8211;459. <div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''   Y.E. Shi, R.K. Ray, K.D. Nguyen, A projection method-based model with the exact C-property for shallow-water ﬂows over dry and irregular bottom using unstructured ﬁnite-volume technique, Comput. Fluids 76 (2013) 178–195. <div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''   Li, Z. & Ito, K. ''The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains'', SIAM: Frontiers in Applied Mathematics, 2006. <div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Z. Li, D. McTigue, and J. Heine, A numerical method for diffusive transport with moving boundaries and discontinuous material properties, Int J Numer Anal Meth Geomech 21 (1997), 653&#8211;662. <div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  Z. Li and J. Zou, Theoretical and numerical analysis on a thermo-elastic system with discontinuities, J Comput Appl Math 91 (1998), 1&#8211;22. <div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  S. Tsynkov, G. Baruch, G. Fibich, & E. Turkel, Fourth order schemes for time-harmonic wave equations with discontinuous coefﬁcients, Commun Comput Phys 5 (2009), 442&#8211;455. <div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''    Li M., Fornberg B. & Tang T. A compact fourth order finite difference scheme for the steady incompressible Navier-Stokes equations, ''Int. J. Numer. Methods Fluids'' '''20''' (1995), 1137&#8211;1151. <div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''      Zhang, J. Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization. ''J. Comput. Phys.'' '''179(1)''' (2002), 170&#8211;179. <div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''    Nabavi, M., Siddiqui, M.H.K., & Dargahi J. A new 9-point sixth-order accurate compact finite difference method for the Helmholtz equation, ''J. Sound Vib.'' '''307''' (2007), 972&#8211;982. <div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''    Wang Y., &  Zhang J. Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D poisson equation, ''J. Comput. Phys.'' '''228''' (2009), 137&#8211;146. <div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''      Zhai, S., Feng, X., & He, Y. A new method to deduce high-order compact difference schemes for two-dimensional Poisson equation. ''Applied Mathematics and Computation'' '''230''' (2014), 9&#8211;26. <div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  Uh Zapata, M., & Itzá Balam, R. High-order implicit finite difference schemes for the two-dimensional Poisson equation. ''Applied Mathematics and Computation'' '''309''' (2017), 222&#8211;244. <div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  Itzá Balam, R., & Uh Zapata, M. A new eighth-order implicit finite difference method to solve the three-dimensional Helmholtz equation. ''Computers & Mathematics with Applications'' '''80(5)''' (2020), 1176&#8211;1200.  <div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''   Sethian, J. A. ''Level Set Methods and Fast Marching Methods: Evolving Interfaces in Geometry, Fluid Mechanics'', Computer Vision, and Materials Sciences, Cambridge University Press, 1999. <div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''   Peskin, C. S. The immersed boundary method, ''Acta Numer.'' '''11''' (2002), 479&#8211;517. <div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  Liu, X. D., & Soderis, T. C. Convergence of the ghost fluid method for elliptic equations with interfaces, ''J. Math. Comp.'' '''72''' (2003), 1731&#8211;1746. <div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  Hu, H., Pan, K., & Tan, Y. An interpolation matched interface and boundary method for elliptic interface problems, ''J. Comput. Appl. Math.'' '''234''' (2010), 73&#8211;94. <div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''   Mu, L., Wang, J., Ye, X., & Zhao, S. A new weak Galerkin finite element method for elliptic interface problems, ''J. Comput. Phys.'' '''325''' (2016), 157&#8211;173. <div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  Cho, H., Han, H., Lee, B., Ha, Y.,  & Kang, M. A second-order boundary condition capturing method for solving the elliptic interface problems on irregular domains, ''Journal of Scientific Computing'' '''81(3)''' (2019), 217&#8211;251. <div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  Itzá Balam, R., Hernandez-Lopez, F., Trejo-Sánchez, J., & Uh Zapata, M. An immersed boundary neural network for solving elliptic equations with singular forces on arbitrary domains. ''Mathematical Biosciences and Engineering: MBE'' '''18(1)''' (2020), 22&#8211;56. <div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''   Leveque R.J., & Li Z. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, ''SIAM J. Numer. Anal.'' '''31(4)''' (1994), 1019-1044. <div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''   Wiegmann A., &  Bube K.P. The explicit-jump immersed interface method: finite difference methods for PDEs with piecewise smooth solutions, ''SIAM J. Numer. Anal.'' '''37(3)''' (2000), 827&#8211;862. <div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''   Berthelsen, P. A. (2004). A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions. J. Comput. Phys. 197(1), 364&#8211;386. <div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''   Seo, J. H. & Mittal, R. A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries, ''J. Comput. Phys.'' '''230''' (2011), 1000&#8211;1019. <div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''    Ito, K., Li, Z., & Kyei, Y. Higher-order, Cartesian grid based finite difference schemes for elliptic equations on irregular domains. ''SIAM Journal on Scientific Computing'' '''27(1)''' (2005), 346-367. <div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''   Gibou, F., & Fedkiw, R. A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem. ''J. Comput. Phys.'' '''202(2)''' (2005), 577-601. <div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''   Linnick, M. N., & Fasel, H. F. A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains. ''J. Comput. Phys.'' '''204(1)''' (2005), 157-192. <div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''    Zhou, Y. C., Zhao, S., Feig, M., & Wei, G. W. High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources. ''J. Comput. Phys.'' '''213(1)''' (2006), 1&#8211;30.  <div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''   Zhong, X. A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity. ''J. Comput. Phys.'' '''225(1)''' (2007), 1066-1099. <div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''   Feng, X., Li, Z., & Qiao, Z. High order compact finite difference schemes for the Helmholtz equation with discontinuous coefficients. ''J. of Comput. Math.'' (2011) 324&#8211;340. <div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''    Pan, K., He, D., & Li, Z. A high order compact FD framework for elliptic BVPs involving singular sources, interfaces, and irregular domains. ''J. Sci. Comput.'' '''88(3)''' (2021), 1&#8211;25. <div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''    Colnago, M., Casaca, W., & de Souza, L. F. A high-order immersed interface method free of derivative jump conditions for Poisson equations on irregular domains. ''J. Comput. Phys.'' '''423''' (2020), 109791. <div id="cite-34"></div>
'''[[#citeF-34|[34]]]'''    Feng, Q., Han, B., & Minev, P. Sixth order compact finite difference schemes for Poisson interface problems with singular sources. ''Computers & Mathematics with Applications'' '''99''' (2021), 2-25. <div id="cite-35"></div>
'''[[#citeF-35|[35]]]'''           Claerbout, J.F. The craft of wave-field extrapolation, in Imaging the Earth's Interior, ''Blackwell Scientific Publications, Oxford'' (1985), 260-265. <div id="cite-36"></div>
'''[[#citeF-36|[36]]]'''    Liu, Y., & Sen, M. K. A practical implicit finite-difference method: examples from seismic modeling. ''Journal of Geophysics and Engineering'' '''6(3)''' (2009),  231. <div id="cite-37"></div>
'''[[#citeF-37|[37]]]'''  W.F. Spotz, High-order compact finite difference schemes for computational mechanics Doctoral dissertation, The University of Texas at Austin, 1995. <div id="cite-38"></div>
'''[[#citeF-38|[38]]]'''   Xu S., & Wang Z.J. Systematic Derivation of Jump Conditions for the Immersed Interface Method in Three-Dimensional Flow Simulation, ''J. Sci. Comput.'' '''27(6)''' (2006), 1948-1980. <div id="cite-39"></div>
'''[[#citeF-39|[39]]]'''  Feng, X., & Li, Z. Simplified immersed interface methods for elliptic interface problems with straight interfaces. ''Num. Meth. for Par. Diff. Eqs.'' '''28(1)''' (2012), 188&#8211;203.  <div id="cite-40"></div>
'''[[#citeF-40|[40]]]'''  Li, Z. Immersed interface methods for moving interface problems. ''Numerical algorithms'' '''14(4)''' (1997), 269-293. <div id="cite-41"></div>
'''[[#citeF-41|[41]]]'''  Huang, H., & Li, Z. (1999). Convergence analysis of the immersed interface method. IMA Journal of Numerical Analysis, 19(4), 583-608.