==General foundations of high-order immersed interface methods   to solve interface problems==

​

'''H. Escamilla Puc<sup>a</sup>, R. Itza Balam<sup>b,c</sup>, M. Uh Zapata<sup>b,c</sup>'''

​

{|

|-

|<sup>a</sup> Facultad de Matemáticas; Anillo Periférico Norte, Tablaje Cat. 13615, Colonia Chuburná Hidalgo Inn, Mérida, México

|}

​

{|

|-

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

|}

​

{|

|-

|<sup>c</sup> Consejo Nacional de Humanidades, Ciencias y Tecnologías, CONAHCYT, México

|}

-->

​

==Abstract==

​

This work forms the foundation for addressing high-order immersed interface methods to solve interface problems and enables us to conduct in-depth examination of this theory. Here, we focus on the introduction a fourth-order finite-difference formulation to approximate the second-order derivative of discontinuous functions. The approach is based on the combination of a high-order implicit formulation and the immersed interface method. The idea is to modify the standard schemes by introducing additional contribution terms based on jump conditions. These contributions are calculated only at grid points where the stencil intersects with the interface. Here, we discuss the issues of implementing the one-dimensional Poisson equation and the heat conduction equation with discontinuous solutions as a three-point stencil for each grid point on the computational domain. In both cases, the resulting discretization approach yields a tridiagonal linear system with matrix coefficients identical to those employed for smooth solutions. We present several numerical experiments to verify the feasibility and accuracy of the method. Thus, this high-order method provides an attractive numerical framework that can efficiently lead to the solution to more complex problems.

​

'''keywords'''

​

Immersed interface method, implicit finite difference, fourth-order accuracy, Poisson equation, heat conduction equation

​

==1 Introduction==

​

High-order numerical solutions to differential equations arising from discontinuous solutions find extensive utility across various research domains <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-1|[1,2,3,4,5,6]]]. In the case of smooth solutions, the standard central finite-difference method requires a significant number of grid points to achieve a high level of accuracy in its numerical results. As a result, over the past few decades, several schemes have been developed to obtain fourth- and sixth-order finite-difference methods <span id='citeF-7'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-7|[7,8,9,10,11]]], including those one based on the implicit finite-difference (IFD) formulation <span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span>[[#cite-12|[12,13,14,15]]].

​

On the other hand, although several methods have been proposed to address discontinuous problems <span id='citeF-2'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-2|[2,16,17,18,19,20,21]]], the Immersed Interface Method (IIM) <span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-25'></span>[[#cite-22|[22,23,24,25]]] stands out as a highly accurate option that requires minimal adjustments to the standard finite-difference formulation. However, these methods typically achieve second-order accuracy. For instance, there are limited implementations of a few third-, fourth- and sixth-order IIMs available for solving Poisson equations with discontinuous solutions <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><span id='citeF-31'></span><span id='citeF-32'></span><span id='citeF-33'></span>[[#cite-26|[26,27,28,29,30,31,32,33]]].

​

This paper focuses on the basic ideas of combining the implicit finite-difference and immersed interface method (IFD-IIM) to achieve high-order approximations for second-order derivatives of both continuous and discontinuous real-valued functions. The IFD scheme offers a highly accurate numerical method <span id='citeF-34'></span><span id='citeF-35'></span>[[#cite-34|[34,35]]], while the IIM handles discontinuities through minimal adjustments made exclusively at grid points where the stencil intersects the interface <span id='citeF-36'></span><span id='citeF-37'></span>[[#cite-36|[36,37]]], yielding additional terms known as jump contributions.

​

We illustrate the implementation of the IFD-IIM approach with two examples: the one-dimensional Poisson equation for static cases and the heat conduction equation with a fixed interface for time-evolving scenarios. Our proposed method offers several advantages. Notably, the resulting tridiagonal matrix coefficient of the finite-difference scheme remains the same as those for smooth solutions, with the additional terms arising from the jumps located in the right-hand side vector. Consequently, our algorithm is straightforward to implement, employing the efficient Thomas' algorithm.

​

We have organized our study as follows. In Section 2, we introduce a fourth-order implicit finite-difference method capable of handling second-order derivatives, both in smooth and discontinuous scenarios. Section 3 demonstrates the application of this implicit scheme in approximating solutions to the one-dimensional Poisson equation. Section 4 shows the combination of the IFD-IIM with the Crank-Nicolson method to solve the heat conduction equation. Sections 5 and 6 provides a series of numerical examples to illustrate the algorithm's accuracy for both equations. Lastly, Section 7 offers our conclusions and outlines directions for future research.

​

==2 Implicit finite-difference formulation with discontinuities==

​

In this section, we outline the key attributes of the implicit finite difference formulation, demonstrating how the scheme can be adapted for addressing discontinuous problems through the utilization of the immersed interface method.

​

We approximate the numerical solution on the domain <math display="inline">[\mathfrak{a},\,\mathfrak{b}]</math> that is divided into <math display="inline">N</math> sub-intervals, as follows

​

<span id="eq-1"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>x_i = \mathfrak{a} + (i-1)h, \qquad i=1, 2, \dots , N,N+1,   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1)

|}

​

where the grid size is given by <math display="inline">h=(\mathfrak{b}-\mathfrak{a})/N</math>. We employ <math display="inline">U_i</math> and <math display="inline">u_i = u(x_i)</math> to denote the approximate and exact solution at the <math display="inline">i</math>th-point of the grid, respectively. Here, the interface is at <math display="inline">x=x_\alpha </math> located between grid points <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> as <math display="inline">x_I\leq x_\alpha{<}x_{I+1}</math>, see Fig.&nbsp;[[#img-1|1]]. The distances of the closest grid points to the interface are defined as <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> and <math display="inline">h_{L}</math> are positive and negative values, respectively.

​

<div id='img-1'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:Draft_Uh Zapata_643844164-Fig_Scheme.png|594px|Example of a discontinuous function u with an interface located between the points x<sub>I</sub> and x<sub>I+1</sub>.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 1:''' Example of a discontinuous function <math>u</math> with an interface located between the points <math>x_{I}</math> and <math>x_{I+1}</math>.

|}

​

On the other hand, the jump for <math display="inline">u</math> at <math display="inline">x_\alpha </math> is defined as

​

<span id="eq-2"></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]= u^{+} - u^{-}, \quad \mathrm{where} \quad u^{+} = u(x_\alpha ^{+} )= \lim _{x\to x_\alpha ^+}u(x),\quad \mathrm{and} \quad  u^{-} = u(x_\alpha ^{-}) = \lim _{x\to x_\alpha ^-}u(x).  </math>

|}

|}

​

We employ a similar definition for the jumps such as the ones of the right-hand side and the derivatives of <math display="inline">u</math>.

​

In this paper, we designate <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> as ''irregular points'', while considering the rest as ''regular points''. This classification holds significant importance as distinct schemes are applied to each category, as presented in the following theorems.

​

===2.1 Second-order derivative approximation for regular and irregular grid points===

​

The following two Theorems state the main results to approximate the second-order derivative using high-order schemes for regular and irregular points.

​

<span id='theorem-teo:regular'></span>Theorem 1:  '''Regular points''' <span id='citeF-34'></span><span id='citeF-35'></span>[[#cite-34|[34,35]]]. Let us consider a real-valued function <math display="inline">u</math> with an interface <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 by the implicit finite-difference (IFD) scheme

​

<span id="eq-2"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\mathfrak{D}^2 u_{xx}\right|_i = \left.\delta ^2 u\right|_i + O(h^4), \quad i \neq I, I+1,    </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2)

|}

​

where

​

<span id="eq-3"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\mathfrak{D}^2(\cdot )  :=  (\cdot )  + bh^2 (\cdot )_{xx},   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (3)

|}

​

<math display="inline">b=1/12</math>, and central finite-difference formula <math display="inline">\delta ^2</math> is given by

​

<span id="eq-4"></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 := \frac{1}{h^2} \left\{u_{i-1}-2u_i+u_{i+1}\right\}.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|}

​

From Taylor series expansions and under some simplifications, the second-order derivative at any regular point <math display="inline">x_i</math> can be written in terms of the centered finite-difference operator, as follows

​

<span id="eq-5"></span>

{| 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 = \frac{1}{h^2} \left\{u_{i-1}-2u_i+u_{i+1}\right\}- \frac{1}{12}h^2\left.u_{xxxx}\right|_i +O(h^4),   </math>

|}

|}

​

Previous equation holds due to the solution is smooth on a neighborhood around <math display="inline">x_{i}</math>. Thus, we obtain a fourth-order IFD using a stencil of three nodes by moving the fourth-order term to the left-hand side. We get

​

{| 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 +  \frac{1}{12}h^2 \left(\left.u_{xx} \right)_{xx}\right|_i = \frac{1}{h^2} \left\{u_{i-1}-2u_i+u_{i+1}\right\}+ O(h^4).  </math>

|}

|}

​

Finally, the proof is completed by using the definition of operator <math display="inline">\mathfrak{D}^2</math>, <math display="inline">\delta ^2</math>, and <math display="inline">b</math>.

​

It is important to remark that formula [[#eq-2|(2)]] is applicable exclusively for regular points. In order to address this limitation for the two irregular points near the interface, we introduce a modified implicit finite-difference scheme using the IIM, specifically tailored to handle discontinuous solutions. Furthermore, instead of having a fourth-order local truncation error for the irregular points, we proceed as other IIMs <span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-29'></span><span id='citeF-33'></span><span id='citeF-32'></span>[[#cite-22|[22,23,24,29,33,32]]] by taking one order lower at these points. We will numerically show that the global order of convergence can be still <math display="inline">O(h^4)</math> even if the local truncation error at <math display="inline">i=I</math> and <math display="inline">i=I+1</math> is <math display="inline">O(h^{3})</math>.

​

<span id='theorem-teo:irregular'></span>Theorem 2:  '''Irregular points''' <span id='citeF-14'></span>[[#cite-14|[14]]]. Let us consider the known jump conditions

​

<span id="eq-5"></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],\quad  \left[u_{x} \right],\quad  \left[u_{xx} \right],\quad  \left[u_{xxx} \right],\quad  \left[u_{xxxx} \right],\quad   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (5)

|}

​

at <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 implicit finite-difference immersed interface method (IFD-IIM) given by

​

<span id="eq-6"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\mathfrak{D}^2 u_{xx}\left.\right|_i = \left.\delta ^2 u\right|_i + C_i + O(h^3), \quad i= I, I+1,    </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|}

​

where <math display="inline">\mathfrak{D}^2</math> and <math display="inline">\delta ^2</math> are defined in [[#eq-3|(3)]] and [[#eq-4|(4)]], respectively, and

​

<span id="eq-7"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>C_i =  \begin{cases}- \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\}, & i=I, \\[3mm]    \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\}, & i=I+1, \end{cases}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|}

​

and where <math display="inline">b=1/12</math>.

​

We obtain a third-order scheme for at <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> following similar ideas as the ones developed for the generalized Taylor expansion proposed by Xu & Wang <span id='citeF-36'></span>[[#cite-36|[36]]] and the IIM for elliptic interface problems with straight interfaces proposed by Feng & Li <span id='citeF-37'></span>[[#cite-37|[37]]]. The idea is to consider extended smooth solutions of <math display="inline">u</math> such that we can apply the standard central scheme to <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math>. For instance, a function based on the original left solution is defined as

​

<span id="eq-8"></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}{6}u_{xxx}^{-}+\dfrac{h_R^4}{24}u_{xxxx}^{-} , & x_\alpha{<}x,\\ \end{cases}  </math>

|}

|}

​

Using Taylor series expansions of <math display="inline">u_{I+1}</math> around <math display="inline">x_\alpha </math>, the definition of jumps [[#eq-5|(5)]], and some simplification yield

​

{| 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.\delta ^2 u_{\ell }\right|_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}{6}\left[u_{xxx} \right]+\dfrac{h_R^4}{24}\left[u_{xxxx} \right]\right)+O(h^3). </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 + bh^2 \left.(u_\ell )_{xxxx} \right|_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}{6}\left[u_{xxx} \right]+\dfrac{h_R^4}{24}\left[u_{xxxx} \right]\right)+O(h^3). </math>

|}

|}

​

Finally, we get [[#eq-6|(6)]] which complete the proof. The same procedure can be applied for the proof at <math display="inline">x_{I+1}</math>. We refer to the reader to the work of Itza Balam and Uh Zapata <span id='citeF-14'></span>[[#cite-14|[14]]] for more details.

​

<span id='theorem-EFD'></span>Remark 1: If <math display="inline">b=0</math>, then Theorem [[#theorem-teo:regular|1]] yields the standard second-order finite-difference method for regular points given by

​

{| 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 + O(h^2), \quad i \neq I, I+1,  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (8)

|}

​

and Theorem [[#theorem-teo:regular|1]] results in an IIM of first-order for irregular points <math display="inline"> i= I, I+1</math> as follows

​

<span id="eq-9"></span>

{| 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_i+ O(h), \quad i=I,I+1,   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (9)

|}

​

where

​

<span id="eq-10"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>c_i =  \begin{cases}- \dfrac{1}{h^2}\left\{[u]+h_R\left[u_{x} \right]+\dfrac{h_R^2}{2}\left[u_{xx} \right]\right\}, & i=I, \\[3mm]    \dfrac{1}{h^2}\left\{[u]+h_L\left[u_{x} \right]+\dfrac{h_L^2}{2}\left[u_{xx} \right]\right\}, & i=I+1. \end{cases}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (10)

|}

​

Note that, 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>.

​

===2.2 Implicit approximation for a real-valued function===

​

Before to apply the previous results for approximations to differential equations, it would be useful to express with finite differences the operator <math display="inline">\left.\mathfrak{D}^2u\right|_i=u_i+\left.bh^2u_{xx}\right|_i</math> for a real-valued function <math display="inline">u</math> and not its second-order derivative <math display="inline">u_{xx}</math>. The finite-difference formula is obtained by approximating <math display="inline">\left.u_{xx}\right|_i</math> with the central finite differences, as presented in [[#theorem-EFD|(1)]]-[[#eq-10|(10)]] for regular and irregular points, respectively.

​

For regular points (<math display="inline">i\neq I,I+1</math>), it follows from equation [[#theorem-EFD|(1)]]:

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\mathfrak{D}^2u\right|_i    =   u_i+\left.bh^2 u_{xx}\right|_i    =   u_i + bh^2\left\{\left.\delta ^2 u \right|_i  + O(h^2)\right\}   =   b u_{i-1}+\,(1-2b) u_i+\,b u_{i+1}  + O(h^4). </math>

|}

|}

​

Thus, if we define as <math display="inline">\mathfrak{d}^2</math> the resulting finite-difference

​

<span id="eq-11"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\mathfrak{d}^2 u \right|_i  : =  b u_{i-1}+\,(1-2b) u_i+\,b u_{i+1},  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (11)

|}

​

then, we have that

​

<span id="eq-12"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\mathfrak{D}^2 u \right|_i = \left.\mathfrak{d}^2\phi \right|_i + O(h^4), \quad i\neq I,I+1.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (12)

|}

​

However, we still have the same issue to overcome for a discontinuous function <math display="inline">u</math>. We remark that this second-order derivative of <math display="inline">u</math> is already multiplied for <math display="inline">h^2</math>. Consequently, we can use the IIM technique where the contribution term is only first-order accurate to <math display="inline">u</math>; thus, we still have a local <math display="inline">O(h^3)</math> to keep a global fourth-order accurate method. Then for irregular points, the IIM applied for this term follows

​

<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>\left.\mathfrak{D}^2 u \right|_i = \left.\mathfrak{d}^2 u \right|_i + (c_u)_i + O(h^3), \quad i=I,I+1,  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (13)

|}

​

where <math display="inline">\mathfrak{d}^2</math> is given by finite difference [[#eq-11|(11)]] and

​

<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>(c_u)_i  = \begin{cases}- b\left([u]+h_R\left[u_x\right]+\dfrac{h_R^2}{2}\left[u_{xx} \right]\right), & i=I, \\[3mm]    b\left([u]+h_L\left[u_x \right]+\dfrac{h_L^2}{2}\left[u_{xx} \right]\right), & i=I+1. \end{cases}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (14)

|}

​

Finally, for regular and irregular points, we have high-order finite-difference approximations of the implicit operator <math display="inline">\mathfrak{D}^2</math> applied to a real-valued function <math display="inline">u</math> as in [[#eq-12|(12)]] and [[#eq-13|(13)]], and to its second-order derivative <math display="inline">u_{xx}</math> as in [[#eq-2|(2)]] and [[#eq-6|(6)]]. Now, we proceed to implement them to the solution of differential equations, as presented in the following two sections.

​

==3 Poisson equation==

​

In this section, we developed a fourth-order finite difference scheme for the Poisson equation. Let us consider <math display="inline">u</math> and <math display="inline">f</math> as the solution of the problem and known right-hand side function, respectively. Thus the interface problem is given by

​

<span id="eq-15"></span>

<span id="eq-15"></span>

<span id="eq-15"></span>

<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>\Delta u(\boldsymbol{x}) = f(\boldsymbol{x}), \quad \boldsymbol{x}\in \Omega \setminus \Gamma , \quad \Omega = \Omega ^{+} \cup \Omega ^{-} </math>

|-

| style="text-align: center;" | <math> u(\boldsymbol{x}) = g(\boldsymbol{x}),  \quad \boldsymbol{x}\in \partial \Omega , </math>

|-

| style="text-align: center;" | <math> \left[u\right]_\Gamma   = \mathfrak{p}(\boldsymbol{x}),  \quad \boldsymbol{x}\in \Gamma , </math>

|-

| style="text-align: center;" | <math> \left[u_{\boldsymbol{n}} \right]_\Gamma = \mathfrak{q}(\boldsymbol{x}) , \quad \boldsymbol{x}\in \Gamma{.}  </math>

|}

|}

​

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>. Dirichlet boundary conditions are defined on <math display="inline">\partial \Omega </math>. We assume that <math display="inline">u</math> and <math display="inline">f</math> may have discontinuities at the interface <math display="inline">\Gamma </math>. Thus, we require additional conditions known as jumps. Note that the principal jump conditions, <math display="inline">\left[u\right]_\Gamma </math> and <math display="inline">\left[u_{\boldsymbol{n}} \right]_\Gamma </math>, are known functions defined on <math display="inline">\Gamma </math>. Here, <math display="inline">u_{\boldsymbol{n}}</math> is the derivative in the normal direction.

​

In the context of the general problem, the computational domain can be considered into multiple dimensions. Nevertheless, since the primary objective of this paper is to illustrate the fundamental attributes of the proposed implicit high-order method, we concentrate on investigating the one-dimensional (1D) Poisson problem as defined by

​

<span id="eq-15"></span>

<span id="eq-16"></span>

<span id="eq-17"></span>

<span id="eq-18"></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 (\mathfrak{a},x_\alpha ) \cup (x_\alpha , \mathfrak{b}), </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (15)

|-

| style="text-align: center;" | <math> u(x) = g(x),  \quad \boldsymbol{x}\in \{ \mathfrak{a},\mathfrak{b}\} , </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (16)

|-

| style="text-align: center;" | <math> \left[u\right] = \mathfrak{p}, </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (17)

|-

| style="text-align: center;" | <math> \left[u_{x}\right] = \mathfrak{q},  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (18)

|}

|}

​

Here, <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 principal jump conditions <math display="inline">\left[u\right]=\left[u\right]_{x_\alpha }</math> and <math display="inline">\left[u_{x}\right]=\left[u_{x}\right]_{x_\alpha }</math> are known values at <math display="inline">x_\alpha </math>.

​

===3.1 IFD-IIM for the 1D Poisson problem===

​

Let us consider that <math display="inline">x_\alpha </math> is located between the adjacent grid points <math display="inline">x_I</math> and <math display="inline">x_{I+1}</math> as <math display="inline">x_I \leq x_\alpha < x_{I+1}</math>. If we apply operator <math display="inline">\mathfrak{D}^2(\cdot ) = (\cdot ) +bh^2 (\cdot )_{xx}</math> at both sides of [[#eq-15|(15)]], then we get

​

<span id="eq-19"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\mathfrak{D}^2 u_{xx}\right|_{i} = \left.\mathfrak{D}^2 f\right|_{i}, \qquad i=2,\dots , N.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (19)

|}

​

For <math display="inline">i=1</math> and <math display="inline">i=N+1</math>, the grid points are at the boundary, thus the Dirichlet boundary condition can be directly applied. Thus, using the IFD scheme [[#eq-2|(2)]] and approximation [[#eq-12|(12)]] in [[#eq-19|(19)]] at <math display="inline">x_i</math>, we have that the finite-difference scheme for regular points is given by

​

<span id="eq-20"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\delta ^2u\right|_{i} = \left.\mathfrak{d}^2f\right|_{i} + O(h^4), \qquad i\neq I,I+1.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (20)

|}

​

Similarly, using formulas [[#eq-6|(6)]] and [[#eq-13|(13)]], the scheme for irregular points is given by

​

<span id="eq-21"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left.\delta ^2u\right|_{i} + C_i = \left.\mathfrak{d}^2f\right|_{i} + (c_f)_i + O(h^3), \qquad i=I,I+1,  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (21)

|}

​

where <math display="inline">C_i</math> and <math display="inline">(c_f)_i</math> are given by [[#eq-7|(7)]] and [[#eq-14|(14)]], respectively. Thus, combining the results for all grid points, the IFD-IIM for the 1D Poisson equation [[#eq-15|(15)]] at <math display="inline">x_i</math> is given by

​

<span id="eq-22"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\frac{1}{h^2}U_{i-1}\,\,-\frac{2}{h^2}U_{i}\,+\,\frac{1}{h^2}U_{i+1}  =  bf_{i-1}+\,(1-2b)f_i+\,bf_{i+1} + \mathfrak{C}_i, \quad i=2, \dots , N,  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (22)

|}

​

where total contribution is <math display="inline">\mathfrak{C}_i = (c_f)_i - C_i</math> with

​

<span id="eq-23"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>C_i =  \begin{cases}- \dfrac{1}{h^2}\left([u]+h_R\left[u_{x}\right]+\dfrac{h_R^2}{2} + b \left(2 h_R^3\left[u_{xxx}\right]+\dfrac{h_R^4}{2}\left[u_{xxxx}\right]\right)\right), & i=I,\\[4mm]   \dfrac{1}{h^2} \left([u]+h_L\left[u_{x}\right]+\dfrac{h_L^2}{2} + b \left(2h_L^3\left[u_{xxx}\right]+\dfrac{h_L^4}{2}\left[u_{xxxx}\right]\right)\right), & i=I+1,\\[2mm]    0 & \textrm{else where,}  \end{cases}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23)

|}

​

and

​

<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>(c_f)_i  = \begin{cases}- b\left([f]+h_R\left[f_x\right]+\dfrac{h_R^2}{2}\left[f_{xx} \right]\right), & i=I, \\[3mm]    b\left([f]+h_L\left[f_x \right]+\dfrac{h_L^2}{2}\left[f_{xx} \right]\right), & i=I+1,\\[2mm]     0 & \textrm{else where.} \end{cases}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (24)

|}

​

We remark that scheme [[#eq-22|(22)]]-[[#eq-24|(24)]] results in an approximation with local truncation error of fourth- and third-order for regular and irregular grid points, respectively. Thus, a global <math display="inline">O(h^4)</math> method is expected.

​

===3.2 IFD-IIM and 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>, and <math display="inline">\left[u_{xxxx}\right]</math> must be known to apply the proposal fourth-order IFD-IIM. Thus, it seems that more jump conditions of <math display="inline">u</math> rather than the principal jump conditions [[#eq-17|(17)]] and [[#eq-18|(18)]] are required to have a fourth-order accurate method. However, we can use the Poisson equation [[#eq-15|(15)]] to obtain them 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], \qquad  \left[u_{xxx}\right]=\left[f_{x}\right], \qquad  \left[u_{xxxx}\right]=\left[f_{xx}\right].  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (25)

|}

​

Thus, the total jump contribution for the one-dimensional Poisson problem is given by

​

<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>\mathfrak{C}_i =  \begin{cases}  \dfrac{1}{h^2}\left([u]+h_R\left[u_{x}\right]+\dfrac{h_R^2}{2} \right)- b \left\{[f] + h_R\left(- \dfrac{2h_R^2}{h^2}+1\right)\left[f_x\right]+\dfrac{h_R^2}{2} \left(-\dfrac{h_R^2}{h^2}+1\right)\left[f_{xx}\right]\right\}, & i=I, \\[2mm] - \dfrac{1}{h^2} \left([u]+h_L\left[u_{x}\right]+\dfrac{h_L^2}{2} \right)+ b \left\{[f] + h_L\left(- \dfrac{2h_L^2}{h^2} +1\right)\left[f_x\right]+ \dfrac{h_L^2}{2} \left(-\dfrac{h_L^2}{h^2}+1\right)\left[f_{xx}\right]\right\}, & i=I+1,\\[2mm] 0, & \textrm{else where}. \end{cases}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (26)

|}

​

Thus contribution <math display="inline">\mathfrak{C}_i</math> depends only on the principal jump conditions, <math display="inline">[u]</math> and <math display="inline">[u_x]</math>, and right-hand side jumps: <math display="inline">[f]</math>, <math display="inline">\left[f_{x}\right]</math>, and <math display="inline">\left[f_{xx} \right]</math>. The additional jumps derivatives from the right-hand side can be approximated using the current values of <math display="inline">f</math>. In this paper, we will assume that we know them.

​

Remark 2: For the 1D Poisson problem, a global second-order IIM (<math display="inline">b=0</math>) only requires to know 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>.

​

Remark 3: If <math display="inline">h_R/h=1</math>, then <math display="inline">h_L=0</math> and both weight terms next to second-order derivative jump of <math display="inline">f</math> are equal to zero in [[#eq-26|(26)]]. Thus, we do not require to know jump condition <math display="inline">\left[f_{xx}\right]</math> to obtain a fourth-order method when the interface is located at a grid point.

​

==4 Heat equation with a fixed interface==

​

For the second differential equation, we consider the heat conduction equation with initial, boundary, and principal jump conditions, as follows

​

<span id="eq-27"></span>

<span id="eq-27"></span>

<span id="eq-27"></span>

<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>u_t(\boldsymbol{x},t) = \Delta u(\boldsymbol{x},t)- f(\boldsymbol{x},t), \quad \boldsymbol{x}\in \Omega \setminus \Gamma \times [0,T], \quad \Omega = \Omega ^{+} \cup \Omega ^{-} </math>

|-

| style="text-align: center;" | <math> u(\boldsymbol{x},0) = u_0(\boldsymbol{x}),  \quad \boldsymbol{x} \in \Omega ,</math>

|-

| style="text-align: center;" | <math> u(\boldsymbol{x},t) = g(\boldsymbol{x},t),  \quad (\boldsymbol{x},t) \in \partial \Omega{.} </math>

|-

| style="text-align: center;" | <math> \left[u\right]_\Gamma   = \mathfrak{p}(\boldsymbol{x},t),  \quad \boldsymbol{x}\in \Gamma \times [0,T], </math>

|-

| style="text-align: center;" | <math> \left[u_{\boldsymbol{n}} \right]_\Gamma = \mathfrak{q}(\boldsymbol{x},t) , \quad \boldsymbol{x}\in \Gamma \times [0,T].  </math>

|}

|}

​

Here, the source <math display="inline">f</math> and initial value <math display="inline">u_0</math> may be discontinuous or singular across a fixed interface <math display="inline">\Gamma </math>. The interface <math display="inline">\Gamma </math> is immersed in the domain <math display="inline">\Omega </math> and divides into two parts, <math display="inline">\Omega ^{+}</math> and <math display="inline">\Omega ^{-}</math>. As the Poisson equation, this paper only focuses on the one-dimensional problem given by

​

<span id="eq-27"></span>

<span id="eq-28"></span>

<span id="eq-29"></span>

<span id="eq-30"></span>

<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>u_t = u_{xx}-f, \quad (x,t) \in (\mathfrak{a},x_\alpha ) \cup (x_\alpha , \mathfrak{b}) \times (0,T],  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (27)

|-

| style="text-align: center;" | <math> u(x,0) = u_0(x), \quad x \in (\mathfrak{a},\mathfrak{b}), </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (28)

|-

| style="text-align: center;" | <math> u(x,t) = g(x,t),  \quad (x,t) \in \{ \mathfrak{a},\mathfrak{b}\} \times (0,T], </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (29)

|-

| style="text-align: center;" | <math> \left[u\right] =  \mathfrak{p}(t),\quad t\in [0,T] </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (30)

|-

| style="text-align: center;" | <math> \left[u_{x}\right] = \mathfrak{q}(t), \quad t\in [0,T].  </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (31)

|}

|}

​

where the source <math display="inline">f</math> and initial value <math display="inline">u_0</math> may be discontinuous or singular across a fixed interface located at <math display="inline">x_\alpha </math>.

​

===4.1 IFD-IIM for the 1D heat conduction problem===

​

Since the interface is fixed and all the quantities are continuous with time, we can approximate the time derivative using the Crank-Nicolson method, as follows

​

<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>u^{n+1} = u^{n} + \dfrac{\Delta t}{2} \left\{\left(u_{xx}-f\right)^{n+1} + \left(u_{xx}-f\right)^{n}\right\}+ O(\Delta t^2).  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (32)

|}

​

Applying the fourth-order operator [[#eq-3|(3)]] to equation [[#eq-32|(32)]] yields

​

<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>\left.\mathfrak{D}^2u\right|_{i}^{n+1}  = \left.\mathfrak{D}^2u\right|_{i}^{n}   + \dfrac{\Delta t}{2}\left\{ \left.\mathfrak{D}^2u_{xx}\right|_{i}^{n+1}  +\left.\mathfrak{D}^2u_{xx}\right|_{i}^{n} -\left.\mathfrak{D}^2f\right|_{i}^{n+1}  -\left.\mathfrak{D}^2f\right|_{i}^{n} \right\} + O(\Delta t^2), \quad i=2, \dots , N.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (33)

|}

​

For regular points, using the IFD method [[#eq-2|(2)]] and approximation [[#eq-12|(12)]], equation [[#eq-33|(33)]] can be approximated 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.\mathfrak{d}^2u\right|_{i}^{n+1}  = \left.\mathfrak{d}^2u\right|_{i}^{n}   + \dfrac{\Delta t}{2}\left\{ \left.\delta ^2u\right|_{i}^{n+1}  +\left.\delta ^2u\right|_{i}^{n} -\left.\mathfrak{d}^2f\right|_{i}^{n+1}  -\left.\mathfrak{d}^2f\right|_{i}^{n} \right\} +  O(\Delta t^2) + O(h^4), \quad i\neq I,I+1.   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (34)

|}

​

For irregular points, using the IFD-IIM [[#eq-6|(6)]] and approximation [[#eq-13|(13)]], the implicit scheme is given by

​

<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>\left.\mathfrak{d}^2u\right|_{i}^{n+1}  = \left.\mathfrak{d}^2u\right|_{i}^{n}  + \dfrac{\Delta t}{2}\left\{ \left.\delta ^2u\right|_{i}^{n+1} +\left.\delta ^2u\right|_{i}^{n} -\left.\mathfrak{d}^2f\right|_{i}^{n+1} -\left.\mathfrak{d}^2f\right|_{i}^{n} \right\}+ \mathfrak{C}_i   +  O(\Delta t^2) + O(h^3), \quad i= I,I+1,  </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>\mathfrak{C}_i = - (c_u)_{i}^{n+1} + (c_u)_{i}^{n} + \dfrac{\Delta t}{2}\left\{ C_{i}^{n+1} + C_{i}^{n} - (c_f)_{i}^{n+1} - (c_f)_{i}^{n} \right\}.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (36)

|}

​

Here, <math display="inline">(c_u)_i</math> and <math display="inline">(c_f)_i</math> are defined as [[#eq-14|(14)]], and <math display="inline">C_i</math> is given by [[#eq-7|(7)]]. Thus, for 1D heat conduction equation [[#eq-27|(27)]], the IFD-IIM reads

​

<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>(b-r)U_{i-1}^{n+1}+(1-2b+2r)U_i^{n+1}+(b-r)U_{i+1}^{n+1}  = (b+r)U_{i-1}^{n}+(1-2b-2r)U_i^{n}+(b+r)U_{i+1}^{n} </math>

|-

| style="text-align: center;" | <math>  - \dfrac{\Delta t}{2}\left\{bf_{i-1}^{n+1}+(1-2b)f_i^{n+1}+bf_{i+1}^{n+1} \right\} </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (37)

|-

| style="text-align: center;" | <math>  - \dfrac{\Delta t}{2}\left\{bf_{i-1}^{n}+(1-2b)f_i^{n}+bf_{i+1}^{n} \right\}+ \mathfrak{C}_i, </math>

|}

|}

​

where <math display="inline">r=\frac{1}{2}\Delta t/h^2</math> and <math display="inline">\mathfrak{C}_i</math> is given by [[#eq-36|(36)]]. Here, we have that <math display="inline">\mathfrak{C}_i=0</math>  for regular points.

​

Remark 4: The IFD-IIM [[#eq-37|(37)]] is unconditionally stable and the local truncation error is of <math display="inline">O(\Delta t^2 + h^4)</math> and <math display="inline">O(\Delta t^2 + h^3)</math> for regular and irregular grid points, respectively. Thus a global fourth-order method is expected by taking <math display="inline">\Delta t = O(h^2)</math>.

​

Remark 5: As the Crank-Nicolson method is implicit in time, the IFD-IIM solves a linear system at each time step. To address this efficiently, we employed Thomas' algorithm, given that the resulting matrix is tridiagonal. Furthermore, the implicit method in space preserves the original structure of this system of equations, thus yielding a higher-order method without compromising the efficiency of the standard scheme.

​

Remark 6: In addition to accounting for the contributions of the source term <math display="inline">f</math> and its derivatives, the finite-difference scheme presented in equation [[#eq-37|(37)]] requires additional knowledge of the jump conditions for the solution given by <math display="inline">\left[u_{xx}\right]</math>, <math display="inline">\left[u_{xxx}\right]</math>, and <math display="inline">\left[u_{xxxx}\right]</math>. It is worth noting that, although not presented here, there are techniques available for deriving all the necessary jump conditions from the principal jump conditions <span id='citeF-14'></span>[[#cite-14|[14]]].

​

==5 Numerical results for the Poisson equation==

​

In this section, we test the IFD-IIM for different examples of the Poisson equation. In the following simulations, we numerically solve the equation for a given right-hand side function and compare it with its analytic solution. In all cases the resulting linear system is solved using the Thomas' algorithm.

​

The numerical method is tested using three different examples. Example 1 considers a smooth solution to verify the fourth-order implicit method for smooth solutions. Example 2 studies a Poisson equation with a discontinuous solution in a single interface point. The Matlab code for this example can be found in Appendix A. Finally, Example 3 presents a discontinuous problem with multiple interface points.

​

For the all these examples, 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> sub-intervals. The errors are reported using the <math display="inline">L_\infty </math>-norm and the discrete <math display="inline">L_2</math>-norm calculated as

​

<span id="eq-38"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left\|e \right\|_\infty = \max _{i,j} \left|u_{i}-U_{i} \right|,   \quad \mathrm{and} \quad  \left\|e\right\|_2 = \left(\sum _{i=1}^{N_x} (u_{i}-U_{i})^2\Delta x\right)^{1/2}, </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (38)

|}

​

respectively, where <math display="inline">U_{i}</math> and <math display="inline">u_{i}</math> corresponds to the numerical and exact solution at <math display="inline">x_i</math>, respectively. The estimated order of accuracy is computed as

​

<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>\textrm{Order} := \dfrac{\log \left(\left\|e \right\|_{N_1}/ \left\|e \right\|_{N_2}\right)}{\log \left(N_1/N_2 \right)}, </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (39)

|}

​

where <math display="inline">N_1</math> and <math display="inline">N_2</math> indicates the different number of sub-intervals.

​

===5.1 Example 1. Poison equation with smooth solution===

​

Example 1 considers an infinitely smooth and differentiable exact solution of the one-dimensional Poisson problem [[#eq-15|(15)]] given by the following expression

​

<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>u(x) = e^{-4\pi (x-1/4)^2}.  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (40)

|}

​

The right-hand side function, <math display="inline">f</math>, is obtained directly from [[#eq-40|(40)]]. We impose Dirichlet boundary conditions according to the function <math display="inline">u</math>.  Due to the regularity of the solution, the jump contributions <math display="inline">\mathcal{C}_I</math> and <math display="inline">\mathcal{C}_{I+1}</math> in equation [[#eq-26|(26)]] are equal to zero.

​

Table [[#table-1|1]] presents the convergence analysis of Example 1 for different grid resolutions. If <math display="inline">b = 0</math>, then we recover the standard central finite-difference method of second-order accuracy. On the other hand, the fourth-order implicit scheme is recovered when <math display="inline">b=1/12</math>. Last row of table [[#table-1|1]] shows the numerical order calculated by the regression-line slope based on a least squares method (LSM) of the <math display="inline">L_{\infty }</math>- and <math display="inline">L_2</math>-norm error. A complete analysis of the IFD method for smooth solutions can be found in <span id='citeF-12'></span>[[#cite-12|[12]]].

​

​

{|  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 testing a Poisson equation with smooth solution.

|-  

| style="text-align: right;" |      

|    

|-<math display="inline"> N </math>      

| <math display="inline"> L_\infty </math>-norm     

|  Order     

| <math display="inline"> L_2</math>-norm     

|  Order    

| <math display="inline"> L_\infty </math>-norm     

|  Order     

| <math display="inline"> L_2</math>-norm     

|  Order    

|-

|  10 

|  7.52e-02 

|  ---   

|  2.84e-02 

|  ---  

|  9.30e-03 

|  ---   

|  3.56e-03 

|  ---   

|-

|  20 

|  1.70e-02 

|  2.15  

|  6.52e-03 

|  2.12 

|  5.05e-04 

|  4.20  

|  1.93e-04 

|  4.21  

|-

|  40 

|  4.15e-03 

|  2.03  

|  1.60e-03 

|  2.03 

|  3.06e-05 

|  4.04  

|  1.17e-05 

|  4.04  

|-

|  80 

|  1.03e-03 

|  2.01  

|  3.98e-04 

|  2.01 

|  1.90e-06 

|  4.01  

|  7.27e-07 

|  4.01  

|-

| 160 

|  2.57e-04 

|  2.00  

|  9.95e-05 

|  2.00 

|  1.18e-07 

|  4.00  

|  4.54e-08 

|  4.00  

|-

| LSM 

|    

|  '''2.04''' 

| 

|  '''2.04''' 

| 

|  '''4.06''' 

| 

|  '''4.06'''