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

Latest revision as of 01:37, 2 December 2023

1 Introduction

High-order numerical solutions to differential equations arising from discontinuous solutions find extensive utility across various research domains [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 [7,8,9,10,11], including those one based on the implicit finite-difference (IFD) formulation [12,13,14,15].

On the other hand, although several methods have been proposed to address discontinuous problems [2,16,17,18,19,20,21], the immersed interface method (IIM) [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 [26,27,28,29,30,31,32,33].

This paper focuses on the basic ideas of combining the IFD and 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 [34,35], while the IIM handles discontinuities through minimal adjustments made exclusively at grid points where the stencil intersects the interface [36,37], yielding additional terms known as jump contributions. The resulting method will be named as implicit finite-difference immersed interface method (IFD-IIM).

We illustrate the implementation of a global fourth-order 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 IFD 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 provide 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 IFD formulation with discontinuities

In this section, we outline the key attributes of the IFD formulation to achieve a global fourth-order method, demonstrating how the scheme can be adapted for addressing discontinuous problems through the utilization of the IIM.

We approximate the numerical solution on the domain ${\textstyle [{\mathfrak {a}},\,{\mathfrak {b}}]}$ that is divided into ${\textstyle N}$ sub-intervals using the points ${\textstyle x_{i}}$ , as follows

(1)

where the grid size is given by ${\textstyle h=({\mathfrak {b}}-{\mathfrak {a}})/N}$ . We employ ${\textstyle U_{i}}$ and ${\textstyle u_{i}=u(x_{i})}$ to denote the approximate and exact solution at the ${\textstyle i}$ th-point of the grid, respectively. Here, the interface is at ${\textstyle x=x_{\alpha }}$ located between grid points ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ as ${\textstyle x_{I}\leq x_{\alpha }{<}x_{I+1}}$ , see Fig. 1. The distances of the closest grid points to the interface are defined as ${\textstyle h_{L}=x_{I}-x_{\alpha }}$ and ${\textstyle h_{R}=x_{I+1}-x_{\alpha }}$ . Note that ${\textstyle h_{R}}$ and ${\textstyle h_{L}}$ are positive and negative values, respectively.

Figure 1: Example of a discontinuous function

with an interface located between the points

and

.

On the other hand, the jump for ${\textstyle u}$ at ${\textstyle x_{\alpha }}$ is defined as

We employ a similar definition for the jumps such as the ones of the right-hand side and the derivatives of ${\textstyle u}$ .

In this paper, we designate ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ 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.

Theorem 1: Regular points [34,35]. Let us consider a real-valued function ${\textstyle u}$ with an interface ${\textstyle x_{\alpha }}$ such that ${\textstyle x_{I}\leq x_{\alpha }{<}x_{I+1}}$ . Then ${\textstyle u_{xx}}$ can be approximated by a fourth-order IFD scheme

(2)

where

(3)

${\textstyle b=1/12}$ , and the central finite-difference formula ${\textstyle \delta ^{2}}$ is given by

(4)

Proof.

From Taylor series expansions and under some simplifications, the second-order derivative at any regular point ${\textstyle x_{i}}$ can be written in terms of the centered finite-difference operator, as follows

Previous equation holds due to the solution is smooth on a neighborhood around ${\textstyle x_{i}}$ . 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

Finally, the proof is completed by using the definition of operator ${\textstyle {\mathfrak {D}}^{2}}$ , ${\textstyle \delta ^{2}}$ , and ${\textstyle b}$ . This completes the proof. ♦

It is important to remark that formula (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 IFD 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 [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 ${\textstyle O(h^{4})}$ even if the local truncation error at ${\textstyle i=I}$ and ${\textstyle i=I+1}$ is ${\textstyle O(h^{3})}$ .

Theorem 2: Irregular points [14]. Let us consider the known jump conditions

(5)

at ${\textstyle x_{\alpha }}$ such that ${\textstyle x_{I}\leq x_{\alpha }{<}x_{I+1}}$ . Then ${\textstyle u_{xx}}$ can be approximated at ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ by the IFD-IIM given by

(6)

where ${\textstyle {\mathfrak {D}}^{2}}$ and ${\textstyle \delta ^{2}}$ are defined in (3) and (4), respectively, and

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}}

(7)

and where ${\textstyle b=1/12}$ .

Proof.

We obtain a third-order scheme for at ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ following similar ideas as the ones developed for the generalized Taylor expansion proposed by Xu & Wang [36] and the IIM for elliptic interface problems with straight interfaces proposed by Feng & Li [37]. The idea is to consider extended smooth solutions of ${\textstyle u}$ such that we can apply the standard central scheme to ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ . For instance, a function based on the original left solution is defined as

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}}

Using Taylor series expansions of ${\textstyle u_{I+1}}$ around ${\textstyle x_{\alpha }}$ , the definition of jumps (5), and some simplification yield

\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}).

Thus,

Finally, we get (6) which complete the proof. The same procedure can be applied for the proof at ${\textstyle x_{I+1}}$ . We refer to the reader to the work of Itza Balam and Uh Zapata [14] for more details. ♦

Remark 1: If ${\textstyle b=0}$ , then Theorem 1 yields the standard second-order finite-difference method for regular points given by

(8)

and Theorem 1 results in an IIM of first-order for irregular points ${\textstyle i=I,I+1}$ as follows

(9)

where

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}}

(10)

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

2.2 Implicit approximation for a real-valued function

Before applying previous results to approximate differential equations, it would be useful to express the operator ${\textstyle {\mathfrak {D}}^{2}u}$ using finite differences for a real-valued function ${\textstyle u}$ rather than its second-order derivative ${\textstyle u_{xx}}$ . The finite-difference formula is obtained by approximating ${\textstyle \left.u_{xx}\right|}$ with the central finite differences, as presented in (1)-(10) for regular and irregular points, respectively.

For regular points ( ${\textstyle i\neq I,I+1}$ ), it follows from equation (1):

Thus, if we define as ${\textstyle {\mathfrak {d}}^{2}}$ the resulting finite-difference

(11)

then, we have that

(12)

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

(13)

where ${\textstyle {\mathfrak {d}}^{2}}$ is given by finite difference (11) and

(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}}

(14)

Finally, for regular and irregular points, we have high-order finite-difference approximations of the implicit operator ${\textstyle {\mathfrak {D}}^{2}}$ applied to a real-valued function ${\textstyle u}$ as in (12) and (13), and to its second-order derivative ${\textstyle u_{xx}}$ as in (2) and (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 develop a fourth-order finite difference scheme for the Poisson equation. Let us consider ${\textstyle u}$ and ${\textstyle f}$ as the solution of the problem and known right-hand side function, respectively. Thus the interface problem is given by

We divide ${\textstyle \Omega }$ in two regions, ${\textstyle \Omega ^{+}}$ and ${\textstyle \Omega ^{-}}$ , separated by an immersed interface ${\textstyle \Gamma }$ . Dirichlet boundary conditions are defined on ${\textstyle \partial \Omega }$ through function ${\textstyle g}$ . We assume that ${\textstyle u}$ and ${\textstyle f}$ may have discontinuities at the interface ${\textstyle \Gamma }$ . Thus, we require additional conditions ${\textstyle \left[u\right]_{\Gamma }}$ and ${\textstyle \left[u_{\boldsymbol {n}}\right]_{\Gamma }}$ known as principal jump conditions. Here, ${\textstyle u_{\boldsymbol {n}}}$ is the derivative in the normal direction. Note that the principal jump conditions, ${\textstyle {\mathfrak {p}}}$ and ${\textstyle {\mathfrak {q}}}$ , are known functions defined on ${\textstyle \Gamma }$ .

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) interface problem as defined by

$u_{xx}(x)=f(x),\quad x\in ({\mathfrak {a}},x_{\alpha })\cup (x_{\alpha },{\mathfrak {b}}),$	(15)
$u(x)=g(x),\quad x\in \{{\mathfrak {a}},{\mathfrak {b}}\},$	(16)
$\left[u\right]={\mathfrak {p}},$	(17)
$\left[u_{x}\right]={\mathfrak {q}},$	(18)

Here, ${\textstyle u}$ and ${\textstyle f}$ can be discontinuous functions at a given fixed point ${\textstyle x=x_{\alpha }}$ , and principal jump conditions ${\textstyle \left[u\right]=\left[u\right]_{x_{\alpha }}}$ and ${\textstyle \left[u_{x}\right]=\left[u_{x}\right]_{x_{\alpha }}}$ are known values at ${\textstyle x_{\alpha }}$ .

3.1 IFD-IIM for the 1D Poisson problem

Let us consider that ${\textstyle x_{\alpha }}$ is located between the adjacent grid points ${\textstyle x_{I}}$ and ${\textstyle x_{I+1}}$ as ${\textstyle x_{I}\leq x_{\alpha }<x_{I+1}}$ . If we apply operator ${\textstyle {\mathfrak {D}}^{2}(\cdot )=(\cdot )+bh^{2}(\cdot )_{xx}}$ at both sides of (15), then we get

(19)

For ${\textstyle i=1}$ and ${\textstyle i=N+1}$ , the grid points are at the boundary, thus the Dirichlet boundary condition can be directly applied. Thus, using the IFD scheme (2) and approximation (12) in (19) at ${\textstyle x_{i}}$ , we have that the finite-difference scheme for regular points is given by

(20)

Similarly, using formulas (6) and (13), the scheme for irregular points is given by

(21)

where ${\textstyle C_{i}}$ and ${\textstyle (c_{f})_{i}}$ are given by (7) and (14), respectively. Thus, combining the results for all grid points, the IFD-IIM for the 1D Poisson equation (15) at ${\textstyle x_{i}}$ is given by

(22)

where total contribution is ${\textstyle {\mathfrak {C}}_{i}=(c_{f})_{i}-C_{i}}$ with

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,\\[4mm]{\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,\\[2mm]0,&{\textrm {elsewhere,}}\end{cases}}

(23)

and

(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 {elsewhere.}}\end{cases}}

(24)

We remark that scheme (22)-(24) results in an approximation with local truncation error of fourth- and third-order for regular and irregular grid points, respectively. Thus, a global ${\textstyle O(h^{4})}$ method is expected.

3.2 IFD-IIM and principal jump conditions

Note that ${\textstyle \left[u\right]}$ , ${\textstyle \left[u_{x}\right]}$ , ${\textstyle \left[u_{xx}\right]}$ , ${\textstyle \left[u_{xxx}\right]}$ , and ${\textstyle \left[u_{xxxx}\right]}$ must be known to apply the proposal fourth-order IFD-IIM. Thus, it seems that more jump conditions of ${\textstyle u}$ rather than the principal jump conditions (17) and (18) are required to have a fourth-order accurate method. However, we can use the Poisson equation (15) to obtain them as follows

(25)

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

{\mathfrak {C}}_{i}={\begin{cases}{\dfrac {1}{h^{2}}}\left([u]+h_{R}\left[u_{x}\right]+{\dfrac {h_{R}^{2}}{2}}\left[f\right]\right)&\\\qquad -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}}\left[f\right]\right)&\\\qquad +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,\\[4mm]0,&{\textrm {elsewhere}}.\end{cases}}

(26)

Thus contribution ${\textstyle {\mathfrak {C}}_{i}}$ depends only on the principal jump conditions, ${\textstyle [u]}$ and ${\textstyle [u_{x}]}$ , and right-hand side jumps: ${\textstyle [f]}$ , ${\textstyle \left[f_{x}\right]}$ , and ${\textstyle \left[f_{xx}\right]}$ . The additional jumps derivatives from the right-hand side can be approximated using the current values of ${\textstyle f}$ . In this paper, we will assume that we know them.

Remark 2: For the 1D Poisson problem, a global second-order IIM ( ${\textstyle b=0}$ ) only requires to know the principal jump conditions ${\textstyle [u]}$ , ${\textstyle \left[u_{x}\right]}$ and ${\textstyle \left[f\right]}$ .

Remark 3: If ${\textstyle h_{R}/h=1}$ , then ${\textstyle h_{L}=0}$ and both weight terms next to second-order derivative jump of ${\textstyle f}$ are equal to zero in (26). Thus, we do not require to know jump condition ${\textstyle \left[f_{xx}\right]}$ 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

where the source ${\textstyle f}$ and initial value ${\textstyle u_{0}}$ may be discontinuous or singular across a fixed interface ${\textstyle \Gamma }$ . The interface ${\textstyle \Gamma }$ is immersed in the domain ${\textstyle \Omega }$ and divides into two parts, ${\textstyle \Omega ^{+}}$ and ${\textstyle \Omega ^{-}}$ . As the Poisson equation, ${\textstyle g}$ , ${\textstyle {\mathfrak {p}}}$ and ${\textstyle {\mathfrak {q}}}$ are known functions. This paper only focuses on the one-dimensional problem given by

$u_{t}=u_{xx}-f,\quad (x,t)\in ({\mathfrak {a}},x_{\alpha })\cup (x_{\alpha },{\mathfrak {b}})\times (0,T],$	(27)
$u(x,0)=u_{0}(x),\quad x\in ({\mathfrak {a}},{\mathfrak {b}}),$	(28)
$u(x,t)=g(x,t),\quad (x,t)\in \{{\mathfrak {a}},{\mathfrak {b}}\}\times [0,T],$	(29)
$\left[u\right]={\mathfrak {p}}(t),\quad t\in [0,T],$	(30)
$\left[u_{x}\right]={\mathfrak {q}}(t),\quad t\in [0,T].$	(31)

where the source ${\textstyle f}$ and initial value ${\textstyle u_{0}}$ may be discontinuous or singular across a fixed interface located at ${\textstyle x_{\alpha }}$ .

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

(32)

Applying the fourth-order operator (3) to equation (32) for ${\textstyle i=2,\dots ,N}$ , yields

\left.{\mathfrak {D}}^{2}u\right|_{i}^{n+1}=\left.{\mathfrak {D}}^{2}u\right|_{i}^{n}+{\dfrac {\Delta t}{2}}\left\{\left.{\mathfrak {D}}^{2}u_{xx}\right|_{i}^{n+1}+\left.{\mathfrak {D}}^{2}u_{xx}\right|_{i}^{n}-\left.{\mathfrak {D}}^{2}f\right|_{i}^{n+1}-\left.{\mathfrak {D}}^{2}f\right|_{i}^{n}\right\}+O(\Delta t^{2}).

(33)

By employing the IFD method (2) and approximation (12), equation (33) can be approximated for regular points as follows

\left.{\mathfrak {d}}^{2}u\right|_{i}^{n+1}=\left.{\mathfrak {d}}^{2}u\right|_{i}^{n}+{\dfrac {\Delta t}{2}}\left\{\left.\delta ^{2}u\right|_{i}^{n+1}+\left.\delta ^{2}u\right|_{i}^{n}-\left.{\mathfrak {d}}^{2}f\right|_{i}^{n+1}-\left.{\mathfrak {d}}^{2}f\right|_{i}^{n}\right\}

(34)

By using the IFD-IIM (6) and approximation (13), the implicit scheme for irregular points is given by

\left.{\mathfrak {d}}^{2}u\right|_{i}^{n+1}=\left.{\mathfrak {d}}^{2}u\right|_{i}^{n}+{\dfrac {\Delta t}{2}}\left\{\left.\delta ^{2}u\right|_{i}^{n+1}+\left.\delta ^{2}u\right|_{i}^{n}-\left.{\mathfrak {d}}^{2}f\right|_{i}^{n+1}-\left.{\mathfrak {d}}^{2}f\right|_{i}^{n}\right\}+{\mathfrak {C}}_{i}

(35)

where

(36)

Here, ${\textstyle (c_{u})_{i}}$ and ${\textstyle (c_{f})_{i}}$ are defined as (14), and ${\textstyle C_{i}}$ is given by (7). Thus, for 1D heat conduction equation (27), the IFD-IIM reads

$(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}$
$-{\dfrac {\Delta t}{2}}\left\{bf_{i-1}^{n+1}+(1-2b)f_{i}^{n+1}+bf_{i+1}^{n+1}\right\}$	(37)
$-{\dfrac {\Delta t}{2}}\left\{bf_{i-1}^{n}+(1-2b)f_{i}^{n}+bf_{i+1}^{n}\right\}+{\mathfrak {C}}_{i},$

where ${\textstyle r={\frac {1}{2}}\Delta t/h^{2}}$ and ${\textstyle {\mathfrak {C}}_{i}}$ is given by (36). Here, we have that ${\textstyle {\mathfrak {C}}_{i}=0}$ for regular points.

Remark 4: The IFD-IIM (37) is unconditionally stable and the local truncation error is of ${\textstyle O(\Delta t^{2}+h^{4})}$ and ${\textstyle O(\Delta t^{2}+h^{3})}$ for regular and irregular grid points, respectively. Thus a global fourth-order method is expected by taking ${\textstyle \Delta t=O(h^{2})}$ .

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 ${\textstyle f}$ and its derivatives, the finite-difference scheme presented in equation (37) requires additional knowledge of the jump conditions for the solution given by ${\textstyle \left[u_{xx}\right]}$ , ${\textstyle \left[u_{xxx}\right]}$ , and ${\textstyle \left[u_{xxxx}\right]}$ . 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 [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. Finally, Example 3 presents a discontinuous problem with multiple interface points. A Matlab code of these examples can be download at https://github.com/CIMATMerida/IFD-IIM

For the all these examples, the computational domain is the interval ${\textstyle [0,1]}$ , and the grid spacing is ${\textstyle h=1/N}$ for different ${\textstyle N}$ sub-intervals. The errors are reported using the ${\textstyle L_{\infty }}$ -norm and the discrete ${\textstyle L_{2}}$ -norm calculated as

(38)

respectively, where ${\textstyle U_{i}}$ and ${\textstyle u_{i}}$ corresponds to the numerical and exact solution at ${\textstyle x_{i}}$ , respectively. The estimated order of accuracy is computed as

(39)

where ${\textstyle N_{1}}$ and ${\textstyle N_{2}}$ indicates the different number of sub-intervals of the corresponding norm.

5.1 Example 1. Poison equation with smooth solution

Example 1 considers the one-dimensional Poisson problem

(40)

where the right-hand side in (40) is given by

The exact solution of problem (40) is an infinitely smooth and differentiable function given by the following expression

(41)

Note that the function ${\textstyle f}$ and Dirichlet boundary conditions are obtained directly from (41). Due to the regularity of the solution, the jump contributions ${\textstyle {\mathcal {C}}_{I}}$ and ${\textstyle {\mathcal {C}}_{I+1}}$ are equal to zero.

Table 1 presents the convergence analysis of Example 1 for different grid resolutions. If ${\textstyle b=0}$ , 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 ${\textstyle b=1/12}$ . Last row of table 1 shows the numerical order calculated by the regression-line slope based on a least squares method (LSM) of the ${\textstyle L_{\infty }}$ - and ${\textstyle L_{2}}$ -norm error. A complete analysis of the IFD method for smooth solutions can be found in [12].

**Table. 1** Convergence analysis of Example 1 testing a Poisson equation with smooth solution.
	${\textstyle b=0}$				${\textstyle b=1/12}$
${\textstyle N}$	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -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

5.2 Example 2. Poison equation with a single interface

For Example 2, we show the method's capability by solving a single interface problem located at ${\textstyle x=x_{\alpha }}$ . Here, we consider the one-dimensional Poisson problem

{\begin{array}{rcl}u_{xx}(x)&=&f(x),\quad x\in (0,x_{\alpha })\cup (x_{\alpha },1),\\[1mm]u(0)&=&0\\[1mm]u(1)&=&-1,\\[1mm]\left[u\right]&=&\cos(\pi x_{\alpha })-\sin(\pi x_{\alpha }),\\[1mm]\left[u_{x}\right]&=&-\pi (\sin(\pi x_{\alpha })+\cos(\pi x_{\alpha })),\end{array}}

(42)

where the right-hand side in (42) is

The exact solution of problem (42) is a discontinuous function given by the expression

(43)

Note that the right-hand side, Dirichlet boundary conditions, and principal jump conditions are obtained directly from equation (43).

We test two different points: ${\textstyle x_{\alpha }=0.4}$ and ${\textstyle x_{\alpha }=0.63}$ . For the case ${\textstyle x_{\alpha }=0.4}$ , we always have ${\textstyle h_{R}/h=1}$ for ${\textstyle N=10\times {2}^{n}}$ ( ${\textstyle n=0,1,2,\dots }$ ); thus the interface is always located at one grid point of that resolution. In general, for ${\textstyle x_{\alpha }=0.63}$ , we have different ${\textstyle h_{R}/h}$ values for different ${\textstyle N}$ numbers. Fig. 2 shows the numerical and exact solution when the interface is located at these two values using ${\textstyle N=40}$ . As expected, the exact solution is accurately recovered for both cases.

Figure 2: Numerical and exact solution of Example 2 using

using (a)

, and (b)

. Here, the interface is located at one grid point using

and

.

Table 2 shows the convergence analysis for Example 2 for the two ${\textstyle x_{\alpha }}$ values. Observe that a second-order method is recovered for ${\textstyle b=0}$ and a fourth-order method for ${\textstyle b=1/12}$ . As expected, the IFD-IIM numerical order does not depend on the location of the interface. However, the magnitude of the error may present minor variations due to the interface position. Fig. 3 shows the error analysis corresponding to interface locations ${\textstyle x_{\alpha }=0.40}$ and ${\textstyle x_{\alpha }=0.63}$ for ${\textstyle N=10,11,12,\dots ,320}$ . The solid black line for this and the following figures shows the numerical order calculated by the regression-line slope based on a least squares method of the ${\textstyle L_{\infty }}$ -norm.

**Table. 2** Convergence analysis of Example 2 using the IFD-IIM.
	${\textstyle b=0}$				${\textstyle b=1/12}$
${\textstyle N}$	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order
10	1.69e-02	–-	5.19e-03	–-	5.75e-05	–-	3.42e-05	–-
20	4.21e-03	2.00	1.18e-03	2.13	3.58e-06	4.00	1.97e-06	4.12
40	1.05e-03	2.00	3.07e-04	1.95	2.24e-07	4.00	1.39e-07	3.82
80	2.63e-04	2.00	7.53e-05	2.03	1.40e-08	4.00	7.75e-09	4.16
160	6.57e-05	2.00	1.86e-05	2.01	8.74e-10	4.00	5.37e-10	3.85
LSM		2.04		2.02		4.00		3.99

$Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) x_α= 0.40, and (b) x_α= 0.63. The solid black line shows the numerical order 4.00 calculated by the regression-line slope based on a least squares method of the L∞-norm.$

Figure 3: Convergence analysis of Example 2 for

using (a)

, and (b)

. The solid black line shows the numerical order

calculated by the regression-line slope based on a least squares method of the

-norm.

The contribution formula includes jumps ${\textstyle [u]}$ , ${\textstyle [u_{x}]}$ , ${\textstyle [u_{xx}]}$ , ${\textstyle [u_{xxx}]}$ , and ${\textstyle [u_{xxxx}]}$ to obtain a fourth-order accurate method. Fig. 4 shows that if we add additional jumps of high-order derivatives into ${\textstyle {\mathcal {C}}}$ , such as ${\textstyle [u_{xxxxx}]=[f_{xxx}]}$ , we observe that the error oscillation decreases in comparison with Fig. 3 results. It is expected because the method is ${\textstyle O(h^{4})}$ for the whole computational domain, including the irregular points. Thus, we can mitigate error oscillations due to interface position by adding high-order jumps.

$Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) α= 0.40, and (b) α= 0.63. The contribution includes jumps up to fifth-order ([uₓₓₓₓₓ]=[fₓₓₓ]).$

Figure 4: Convergence analysis of Example 2 for

using (a)

, and (b)

. The contribution includes jumps up to fifth-order (

).

Now we study the effect of removing jumps in ${\textstyle {\mathcal {C}}}$ . Let us consider a contribution term up to the third derivative jump by dropping ${\textstyle [u_{xxxx}]=[f_{xx}]}$ . Thus, the numerical approximation is only third-order accurate, as shown in Fig. 5 for different interface values. Note that error oscillations may have a clear pattern or a scattered distribution depending on the interface location. In general, the error magnitude perturbations are related to the variations coming from ${\textstyle h_{R}/h}$ and ${\textstyle h_{L}/h}$ values in the ${\textstyle {\mathcal {C}}}$ formula (26).

$Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) α= 0.40, and (b) α= 0.63. The contribution includes jumps up to third-order ([uₓₓₓ]=[fₓ]).$

Figure 5: Convergence analysis of Example 2 for

using (a)

, and (b)

. The contribution includes jumps up to third-order (

).

Note that there are some points in Fig. 5, marked with circles, that are close to the fourth-order line. Those circled markers correspond to the values with ${\textstyle h_{R}/h=1}$ : (a) ${\textstyle N=5,10,15,\dots ,315}$ for ${\textstyle x_{\alpha }=0.40}$ , and (b) ${\textstyle N=100,200,300}$ for ${\textstyle x_{\alpha }=0.63}$ . Both set of points satisfy that the interface ${\textstyle x_{\alpha }}$ is located at a grid point ${\textstyle x_{I}}$ . In the case of ${\textstyle x_{\alpha }=0.4}$ , we observe that the global order (black points) is close to ${\textstyle 3.21}$ ; meanwhile, the circled points order is four. Similar behavior is obtained for ${\textstyle x_{\alpha }=0.63}$ . Thus, for a given mesh resolution ${\textstyle N}$ , the IFD-IIM is still fourth-order accurate for ${\textstyle h_{R}/h=1}$ , as proven theoretically in Section 4.2.

5.3 Example 3. Poison equation with multiple interface points

Example 3 investigates the numerical scheme capability to solve a multiple interface problem. We only focus on two interface points located at ${\textstyle x=x_{\alpha _{1}}}$ and ${\textstyle x=x_{\alpha _{2}}}$ . However, the methodology could be applied for multiple interfaces by doing minor modifications in the implementation. For this problem, the exact solution of the Poisson problem is given by

(44)

The right-hand function, Dirichlet boundary conditions, and principal jump conditions are obtained directly from the exact solution (44). We consider the same computational domain and grid resolution as previous 1D examples. Fig. 6 presents the analytical and numerical solution when the interface is located at ${\textstyle x_{\alpha _{1}}=0.3}$ and ${\textstyle x_{\alpha _{2}}=0.7}$ using ${\textstyle N=40}$ . This figure also shows the error analysis. As expected, the IFD-IIM is a fourth-order accurate method.

Figure 6: (a) Numerical and exact solution of Example 3 with multiple interfaces using

, and (b) convergence error analysis using grid resolutions of

, 20, 40, 80 and 160.

6 Numerical results for the Heat equation

This section tests the IFD-IIM for the Heat conduction equation in different scenarios. Example 4 verifies the fourth-order implicit method for smooth solutions in space. Example 5 studies a Heat equation with a discontinuous solution in a single interface point and no source term. Finally, Example 6 presents a general discontinuous problem. For the all these examples, the computational domain is the same interval ${\textstyle [0,1]}$ and space step, ${\textstyle h}$ , used for the Poisson examples. Here, the final time is at ${\textstyle T=0.5}$ and the time step is ${\textstyle \Delta t={\frac {1}{4}}h^{2}}$ . The error and estimated order of accuracy are reported using (38) and (39), respectively, at the final step approximation.

6.1 Example 4. Heat equation with smooth solution

This example considers the heat conduction problem

{\begin{array}{rcll}u_{t}(x,t)&=&u_{xx}(x,t)-f(x,t),&(x,t)\in (0,1)\times (0,{\frac {1}{2}}],\\[2mm]u(x,0)&=&\sin \left({\frac {3}{2}}\pi x\right),&x\in (0,1),\\[2mm]u(0,t)&=&0,&t\in [0,{\frac {1}{2}}],\\[2mm]u(1,t)&=&-e^{-t},&t\in [0,{\frac {1}{2}}],\end{array}}

(45)

where ${\textstyle f}$ in (45) is given by

This example is constructed so that the exact solution is

(46)

Note that the source term, ${\textstyle f}$ , initial condition and boundary conditions are obtained from the exact solution (46). The convergence analysis of this example is presented in Table 3 with the final time being ${\textstyle T=0.5}$ . As expected, the high-order methods reach their corresponding order of accuracy.

**Table. 3** Convergence analysis of Example 4 testing a heat equation with smooth solution.
	${\textstyle b=0}$				${\textstyle b=1/12}$
${\textstyle N}$	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order
10	1.66e-02	–-	1.07e-02	–-	1.84e-04	–-	1.18e-04	–-
20	4.12e-03	2.01	2.65e-03	2.01	1.14e-05	4.01	7.34e-06	4.01
40	1.03e-03	2.00	6.60e-04	2.00	7.15e-07	4.00	4.58e-07	4.00
80	2.58e-04	2.00	1.65e-04	2.00	4.47e-08	4.00	2.86e-08	4.00
160	6.44e-05	2.00	4.12e-05	2.00	2.79e-09	4.00	1.79e-09	4.00
LSM		2.00		2.00		4.00		4.00

6.2 Example 5. Heat equation with discontinuous solution

In this example, the heat equation with discontinuous solution is given by

{\begin{array}{rcll}u_{t}(x,t)&=&u_{xx}(x,t),&(x,t)\in (0,x_{\alpha })\cup (x_{\alpha },1)\times (0,{\frac {1}{2}}],\\[2mm]u(x,0)&=&{\begin{cases}\sin(x),&x<x_{\alpha },\\\cos(x),&x>x_{\alpha },\\\end{cases}}&x\in (0,x_{\alpha })\cup (x_{\alpha },1),\\\\[-3mm]u(0,t)&=&0,&t\in [0,{\frac {1}{2}}],\\[2mm]u(1,t)&=&\cos(1)e^{-t},&t\in [0,{\frac {1}{2}}],\\[2mm]\left[u\right](t)&=&\left(\cos(x_{\alpha })-\sin(x_{\alpha })\right)e^{-t},&t\in [0,{\frac {1}{2}}],\\[2mm]\left[u_{x}\right](t)&=&-\left(\sin(x_{\alpha })+\cos(x_{\alpha })\right)e^{-t},&t\in [0,{\frac {1}{2}}].\end{array}}

(47)

The exact solution ${\textstyle u}$ is defined as

(48)

Similar to previous examples, we use ${\textstyle x_{\alpha }=0.4}$ and ${\textstyle x_{\alpha }=0.63}$ . Note that the source term ${\textstyle f}$ of this problem is ${\textstyle f\equiv 0}$ . In this example, both the function ${\textstyle u}$ and their derivatives have jumps, and these jumps vary in time. The boundary condition, initial condition, and all jumps are specified by ${\textstyle u}$ .

The figure of numerical solution and absolute error plot using ${\textstyle N=40}$ for different time stages are shown in Fig. 7. On the other hand, the one-dimensional results at ${\textstyle T=0.5}$ are presented in Figs. 8 and 9. As expected, we can accurately recover the exact solution using the proposed method for ${\textstyle N=40}$ and considering different interface values. Additionally, we provide the absolute error for both cases, noting that the maximum error is found in close proximity to the point ${\textstyle x_{\alpha }}$ . This behavior is also expected since the local truncation error for irregular points is less accurate than for regular points. Finally, the convergence analysis using ${\textstyle N=10}$ , 20, 40, 80 and 160 confirms the fourth-order global convergence of the proposed method.

More details of the grid refinement analysis at ${\textstyle T=0.5}$ is presented in Table 4 for two different interface point locations. As expected, a fourth-order method is obtained for both interfaces. We remark that the variation of the errors using ${\textstyle x_{\alpha }=0.63}$ is because the different ${\textstyle h_{R}/h}$ values resulting from the discretization.

Figure 7: Numerical solution and absolute error of Example 5 using

,

for

.

Figure 8: Numerical results at final time simulation

and

of Example 5.

Figure 9: Numerical results at final time simulation

and

of Example 5.

**Table. 4** Convergence analysis of Example 5 using the IFD-IIM for the heat equation with a discontinuous solution.
	${\textstyle x_{\alpha }=0.4}$				${\textstyle x_{\alpha }=0.63}$
${\textstyle N}$	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order
10	1.11e-07	–-	6.52e-08	–-	7.57e-08	–-	4.49e-08	–-
20	6.90e-09	4.01	4.01e-09	4.02	3.75e-09	4.34	2.25e-09	4.32
40	4.29e-10	4.00	2.49e-10	4.01	3.58e-10	3.39	2.09e-10	4.43
80	2.68e-11	4.00	1.55e-11	4.00	1.53e-11	4.55	8.93e-12	4.55
160	1.63e-12	4.03	9.39e-13	4.04	1.35e-12	3.50	7.82e-13	3.51
LSM		4.01		4.02		3.95		3.96

From Table 5, we can find the results when the time step is chosen as ${\textstyle \Delta t=h}$ . As expected, the proposed IFD-IIM is stable and has two-order convergence either we take a second- or fourth-order approximation in space. We remark that this is a property of the selected time integration method. For instance, similar findings as Table 5 can be obtained for smooth solutions.

**Table. 5** Convergence analysis of Example 5 testing the heat equation with $\Delta t=h$ .
	${\textstyle b=0}$				${\textstyle b=1/12}$
${\textstyle N}$	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order
10	2.94e-04	–-	1.87e-04	–-	3.78e-05	–-	2.67e-05	–-
20	8.53e-05	1.78	4.85e-05	1.95	1.03e-05	1.87	7.22e-06	1.89
40	2.11e-05	2.01	1.26e-05	1.94	2.74e-06	1.91	1.92e-06	1.91
80	5.31e-06	1.99	3.15e-06	2.00	6.91e-07	1.99	4.83e-07	1.99
160	1.34e-06	1.99	7.83e-07	2.00	1.72e-07	2.00	1.21e-07	2.00
LSM		1.96		1.97		1.95		1.95

6.3 Example 6. Heat equation with a general discontinuous solution

Finally, Example 6 investigates the capability to solve a heat interface problem with a general discontinuous solution. For this problem, we slightly modified the previous example such that the exact solution of the heat conduction problem is given by

(49)

Here, the source term ${\textstyle f}$ is a general function that varies both in time and space. It is directly derived from equation (27) as well as the exact solution (49). Furthermore, the exact solution is utilized to specify the boundary condition, initial condition, and all jumps contributions.

In Fig. 10, we depict the temporal evolution of the numerical solution and absolute errors using parameters ${\textstyle N=40}$ and ${\textstyle x_{\alpha }=0.4}$ . It is noteworthy that the behavior of the solution differs from the previous example. More in-depth analysis of the results at ${\textstyle T=0.5}$ are illustrated in Figs. 11 and 12, corresponding to values of ${\textstyle x_{\alpha }=0.4}$ and ${\textstyle x_{\alpha }=0.63}$ , respectively. First, the numerical solution is accurately approximated using the proposed implicit method for ${\textstyle N=40}$ , regardless of the interface locations. It is worth noting that the maximum error for ${\textstyle x_{\alpha }=0.63}$ is concentrated near to this point; however, this is not observed for ${\textstyle x_{\alpha }=0.4}$ . This difference in behavior can be attributed to the interplay between the interface location ( ${\textstyle h_{R}/h=1}$ and ${\textstyle x_{\alpha }=0.4}$ ) and the non-zero right-hand side values, which contribute to the local truncation errors. Grid refinement analyses are also provided in these figures using ${\textstyle N=10}$ , 20, 40, 80 and 160. As anticipated, the IFD-IIM method demonstrates fourth-order accuracy, a validation that is further detailed in Table 6.

Figure 10: Numerical solution and absolute error of Example 6 using

,

for

.

Figure 11: Numerical results at final time simulation

and

of Example 6.

Figure 12: Numerical results at final time simulation

and

of Example 6.

**Table. 6** Convergence analysis of Example 6 using the IFD-IIM for the heat equation with a general discontinuous solution.
	${\textstyle x_{\alpha }=0.4}$				${\textstyle x_{\alpha }=0.63}$
${\textstyle N}$	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order	${\textstyle L_{\infty }}$ -norm	Order	${\textstyle L_{2}}$ -norm	Order
10	3.91e-04	–-	2.18e-04	–-	4.32e-04	–-	2.43e-04	–-
20	2.50e-05	3.97	1.34e-05	4.02	2.48e-05	4.12	1.45e-05	4.07
40	1.56e-06	4.00	8.37e-07	4.01	2.38e-06	3.38	1.12e-06	3.69
80	9.72e-08	4.00	5.22e-08	4.00	9.48e-08	4.65	5.56e-08	4.34
160	6.08e-09	4.03	3.26e-09	4.00	9.29e-09	3.35	4.37e-09	3.69
LSM		4.00		4.01		3.90		3.95

7 Conclusions

This work serves as the general foundation for addressing high-order IIMs to solve interface problems, facilitating comprehensive investigations into this theoretical framework. Here, we present a fourth-order finite-difference scheme for approximating the second-order derivative of real-valued continuous and discontinuous functions. This method combines an implicit formulation with the IIM. Our proposed scheme employs a three-point stencil, achieving different accuracy at regular and irregular points. To illustrate the effectiveness of the proposed IFD-IIM approach, we applied it to solve the one-dimensional Poisson equation and the heat conduction equation. The global accuracy of the fourth order was demonstrated using several numerical examples for both equations. Hence, this study establishes a general strategy for high-order IIMs, enabling their application to elliptic and time-evolving problems in several dimensions. Additionally, the implicit procedure lends itself to developing of higher-order numerical schemes based on the IIM, including sixth-order methods.

Acknowledgments

This work was partially supported by CONAHCYT under the program Investigadoras e Investigadores por México. We also thanks the facilities provided by the Facultad de Matemáticas (UADY) for the development of this work.

BIBLIOGRAPHY

[1] Diersch, H-J. (1990) "Fletcher, CAJ, Computational Techniques for Fluid Dynamics. Vol. I: Fundamental and General Techniques. Vol. II: Specific Techniques for Different Flow Categories. Berlin etc., Springer-Verlag 1988. XIV, 409 pp., 183 figs./XI, 484 pp., 183 figs., DM 198, 00 as a Set. ISBN 3-540-18151-2/3-540-18759-6 (Springer Series in Computational Physics)". Wiley Online Library

[2] Sethian, James A and others. (1999) "Level set methods and fast marching methods", Volume 98. Cambridge Cambridge UP 2

[3] Li, Zhilin and Ito, Kazufumi. (2006) "The immersed interface method: numerical solutions of PDEs involving interfaces and irregular domains". SIAM

[4] Javierre, E and Vuik, C and Vermolen, FJ and Van der Zwaag, S. (2006) "A comparison of numerical models for one-dimensional Stefan problems", Volume 192. Elsevier. Journal of Computational and Applied Mathematics 2 445–459

[5] Shi, Yu-e and Ray, Rajendra K and Nguyen, Kim Dan. (2013) "A projection method-based model with the exact C-property for shallow-water flows over dry and irregular bottom using unstructured finite-volume technique", Volume 76. Elsevier. Computers & Fluids 178–195

[6] Uh, Miguel and Xu, Sheng. (2014) "The immersed interface method for simulating two-fluid flows", Volume 7. Cambridge University Press. Numerical Mathematics: Theory, Methods and Applications 4 447–472

[7] Li, Ming and Tang, Tao and Fornberg, Bengt. (1995) "A compact fourth-order finite difference scheme for the steady incompressible Navier-Stokes equations", Volume 20. Wiley Online Library. International Journal for Numerical Methods in Fluids 10 1137–1151

[8] Zhang, Jun. (2002) "Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization", Volume 179. Elsevier. Journal of Computational Physics 1 170–179

[9] Nabavi, Majid and Siddiqui, MH Kamran and Dargahi, Javad. (2007) "A new 9-point sixth-order accurate compact finite-difference method for the Helmholtz equation", Volume 307. Elsevier. Journal of Sound and Vibration 3-5 972–982

[10] Wang, Yin and Zhang, Jun. (2009) "Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation", Volume 228. Elsevier. Journal of Computational Physics 1 137–146

[11] Zhai, Shuying and Feng, Xinlong and He, Yinnian. (2014) "A new method to deduce high-order compact difference schemes for two-dimensional Poisson equation", Volume 230. Elsevier. Applied Mathematics and Computation 9–26

[12] Uh Zapata, Miguel and Itzá Balam, Reymundo. (2017) "High-order implicit finite difference schemes for the two-dimensional Poisson equation", Volume 309. Elsevier. Applied Mathematics and Computation 222–244

[13] Itzá Balam, Reymundo and Uh Zapata, Miguel. (2020) "A new eighth-order implicit finite difference method to solve the three-dimensional Helmholtz equation", Volume 80. Elsevier. Computers & Mathematics with Applications 5 1176–1200

[14] Itza Balam, Reymundo and Uh Zapata, Miguel. (2022) "A fourth-order compact implicit immersed interface method for 2D Poisson interface problems", Volume 119. Elsevier. Computers & Mathematics with Applications 257–277

[15] Zapata, M Uh and Balam, R Itza and Montalvo-Urquizo, J. (2023) "A compact sixth-order implicit immersed interface method to solve 2D Poisson equations with discontinuities", Volume 210. Elsevier. Mathematics and Computers in Simulation 384–407

[16] Peskin, Charles S. (2002) "The immersed boundary method", Volume 11. Cambridge University Press. Acta numerica 479–517

[17] Liu, Xu-Dong and Sideris, Thomas. (2003) "Convergence of the ghost fluid method for elliptic equations with interfaces", Volume 72. Mathematics of computation 244 1731–1746

[18] Pan, Kejia and Tan, Yongji and Hu, Hongling. (2010) "An interpolation matched interface and boundary method for elliptic interface problems", Volume 234. Elsevier. Journal of computational and applied mathematics 1 73–94

[19] Mu, Lin and Wang, Junping and Ye, Xiu and Zhao, Shan. (2016) "A new weak Galerkin finite element method for elliptic interface problems", Volume 325. Elsevier. Journal of Computational Physics 157–173

[20] Cho, Hyuntae and Han, Heejae and Lee, Byungjoon and Ha, Youngsoo and Kang, Myungjoo. (2019) "A second-order boundary condition capturing method for solving the elliptic interface problems on irregular domains", Volume 81. Springer. Journal of Scientific Computing 217–251

[21] Itzá Balam, Reymundo and Hernandez-Lopez, Francisco and Trejo-Sánchez, Joel and Uh Zapata, Miguel. (2021) "An immersed boundary neural network for solving elliptic equations with singular forces on arbitrary domains", Volume 18. Math. Biosci. Eng 1 22–56

[22] LeVeque, Randall J and Li, Zhilin. (1994) "The immersed interface method for elliptic equations with discontinuous coefficients and singular sources", Volume 31. SIAM. SIAM Journal on Numerical Analysis 4 1019–1044

[23] Wiegmann, Andreas and Bube, Kenneth P. (2000) "The explicit-jump immersed interface method: finite difference methods for PDEs with piecewise smooth solutions", Volume 37. SIAM. SIAM Journal on Numerical Analysis 3 827–862

[24] Berthelsen, Petter Andreas. (2004) "A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions", Volume 197. Elsevier. Journal of Computational Physics 1 364–386

[25] Seo, Jung Hee and Mittal, Rajat. (2011) "A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries", Volume 230. Elsevier. Journal of computational physics 4 1000–1019

[26] Ito, Kazufumi and Li, Zhilin and Kyei, Yaw. (2005) "Higher-order, Cartesian grid based finite difference schemes for elliptic equations on irregular domains", Volume 27. SIAM. SIAM Journal on Scientific Computing 1 346–367

[27] Gibou, Frédéric and Fedkiw, Ronald. (2005) "A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem", Volume 202. Elsevier. Journal of Computational Physics 2 577–601

[28] Linnick, Mark N and Fasel, Hermann F. (2005) "A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains", Volume 204. Elsevier. Journal of Computational Physics 1 157–192

[29] Zhou, YC and Zhao, Shan and Feig, Michael and Wei, Guo-Wei. (2006) "High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources", Volume 213. Elsevier. Journal of Computational Physics 1 1–30

[30] Zhong, Xiaolin. (2007) "A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity", Volume 225. Elsevier. Journal of Computational Physics 1 1066–1099

[31] Feng, Xiufang and Li, Zhilin and Qiao, Zhonghua. (2011) "High order compact finite difference schemes for the Helmholtz equation with discontinuous coefficients". JSTOR. Journal of Computational Mathematics 324–340

[32] Pan, Kejia and He, Dongdong and Li, Zhilin. (2021) "A high order compact FD framework for elliptic BVPs involving singular sources, interfaces, and irregular domains", Volume 88. Springer. Journal of Scientific Computing 3 67

[33] Colnago, Marilaine and Casaca, Wallace and de Souza, Leandro Franco. (2020) "A high-order immersed interface method free of derivative jump conditions for Poisson equations on irregular domains", Volume 423. Elsevier. Journal of Computational Physics 109791

[34] Claerbout, Jon F. (1985) "The craft of wave-field extrapolation, in Imaging the earth's interior", Volume 1. Blackwell scientific publications Oxford 260–265

[35] Liu, Yang and Sen, Mrinal K. (2009) "A practical implicit finite-difference method: examples from seismic modelling", Volume 6. Oxford University Press. Journal of Geophysics and Engineering 3 231–249

[36] Xu, Sheng and Wang, Z Jane. (2006) "Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation", Volume 27. SIAM. SIAM Journal on Scientific Computing 6 1948–1980

[37] Feng, Xiufang and Li, Zhilin. (2012) "Simplified immersed interface methods for elliptic interface problems with straight interfaces", Volume 28. Wiley Online Library. Numerical Methods for Partial Differential Equations 1 188–203

Latest revision as of 01:37, 2 December 2023

1 Introduction

2 IFD formulation with discontinuities

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

2.2 Implicit approximation for a real-valued function

3 Poisson equation

3.1 IFD-IIM for the 1D Poisson problem

3.2 IFD-IIM and principal jump conditions

4 Heat equation with a fixed interface

4.1 IFD-IIM for the 1D heat conduction problem

5 Numerical results for the Poisson equation

5.1 Example 1. Poison equation with smooth solution

5.2 Example 2. Poison equation with a single interface

5.3 Example 3. Poison equation with multiple interface points

6 Numerical results for the Heat equation

6.1 Example 4. Heat equation with smooth solution

6.2 Example 5. Heat equation with discontinuous solution

6.3 Example 6. Heat equation with a general discontinuous solution

7 Conclusions

Acknowledgments

BIBLIOGRAPHY

Document information

Document Score

Share this document

Keywords

claim authorship

@@ Line 2: / Line 2: @@
 ==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>'''
+'''Heidy Escamilla Puc<sup>a</sup>, Reymundo Itzá Balam<sup>b,c</sup>, Miguel Uh Zapataangeluh@cimat.mx<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==
@@ Line 32: / Line 9: @@
 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]]].
+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.
+This paper focuses on the basic ideas of combining the IFD and 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. The resulting method will be named as implicit finite-difference immersed interface method (IFD-IIM).
-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 illustrate the implementation of a global fourth-order 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.
+We have organized our study as follows. In Section 2, we introduce a fourth-order IFD 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 provide 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==
+==2 IFD 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.
+In this section, we outline the key attributes of the IFD formulation to achieve a global fourth-order method, demonstrating how the scheme can be adapted for addressing discontinuous problems through the utilization of the IIM.
-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
+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 using the points <math display="inline">x_i</math>, as follows
 <span id="eq-1"></span>
@@ Line 62: / Line 39: @@
 {| 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>.]]
+|[[Image:Review_499745022132-Fig01_Scheme.png|600px|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>.
@@ Line 87: / Line 64: @@
 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='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 a fourth-order IFD scheme
 <span id="eq-2"></span>
@@ Line 113: / Line 90: @@
 |}
-<math display="inline">b=1/12</math>, and central finite-difference formula <math display="inline">\delta ^2</math> is given by
+<math display="inline">b=1/12</math>, and the central finite-difference formula <math display="inline">\delta ^2</math> is given by
 <span id="eq-4"></span>
@@ Line 125: / Line 102: @@
 | style="width: 5px;text-align: right;white-space: nowrap;" | (4)
 |}
+'''Proof.'''
 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
@@ Line 149: / Line 128: @@
 |}
-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>.
+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>. This completes the proof. ♦
-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>.
+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 IFD 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='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>
@@ Line 166: / Line 145: @@
 |}
-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
+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 IFD-IIM given by
 <span id="eq-6"></span>
@@ Line 194: / Line 173: @@
 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
+'''Proof.'''
+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>
@@ Line 224: / Line 205: @@
 {| 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>
+| 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 </math>
+|-
+| style="text-align: center;" | <math>  +\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.
+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
+<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;"
@@ Line 263: / Line 246: @@
 {| 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="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)
@@ Line 272: / Line 255: @@
 ===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.
+Before applying previous results to approximate differential equations, it would be useful to express the operator <math display="inline">\mathfrak{D}^2u</math> using finite differences  for a real-valued function <math display="inline">u</math> rather than 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|</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)]]:
@@ Line 306: / Line 289: @@
 {| 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="text-align: center;" | <math>\left.\mathfrak{D}^2 u \right|_i = \left.\mathfrak{d}^2 u \right|_i + O(h^4), \quad i\neq I,I+1.  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (12)
@@ Line 341: / Line 324: @@
 ==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
+In this section, we develop 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>
@@ Line 352: / Line 335: @@
 {| 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>\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>
@@ Line 362: / Line 345: @@
 |}
-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.
+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> through function <math display="inline">g</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 <math display="inline">\left[u\right]_\Gamma </math> and <math display="inline">\left[u_{\boldsymbol{n}} \right]_\Gamma </math> known as principal jump conditions. Here, <math display="inline">u_{\boldsymbol{n}}</math> is the derivative in the normal direction. Note that the principal jump conditions, <math display="inline">\mathfrak{p}</math> and <math display="inline">\mathfrak{q}</math>, are known functions defined on <math display="inline">\Gamma </math>.
-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
+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) interface problem as defined by
 <span id="eq-15"></span>
@@ Line 378: / Line 361: @@
 | 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="text-align: center;" | <math> u(x) = g(x),  \quad x\in \{ \mathfrak{a},\mathfrak{b}\} , </math>
 | style="width: 5px;text-align: right;white-space: nowrap;" | (16)
 |-
@@ Line 389: / Line 372: @@
 |}
-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>.
+Here, <math display="inline">u</math> and <math display="inline">f</math> can be discontinuous functions at a given fixed 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===
@@ Line 453: / Line 436: @@
 {| 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="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(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}\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,\\[2mm]    0, & \textrm{else where,}  \end{cases}  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (23)
@@ Line 466: / Line 449: @@
 {| 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="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)
@@ Line 495: / Line 478: @@
 {| 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="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}\left[f\right] \right)& \\ \qquad  - 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}\left[f\right] \right)& \\ \qquad  + 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,\\[4mm]    0, & \textrm{else where}. \end{cases}  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (26)
@@ Line 502: / Line 485: @@
 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 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.
+'''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==
@@ Line 519: / Line 502: @@
 {| 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_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> u(\boldsymbol{x},t) = g(\boldsymbol{x},t),  \quad (\boldsymbol{x},t) \in \partial \Omega \times [0,T], </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>
+| 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
+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 <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, <math display="inline">g</math>, <math display="inline">\mathfrak{p}</math> and <math display="inline">\mathfrak{q}</math> are known functions. This paper only focuses on the one-dimensional problem given by
 <span id="eq-27"></span>
@@ Line 549: / Line 532: @@
 | 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="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="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)
 |-
@@ Line 577: / Line 560: @@
 |}
-Applying the fourth-order operator [[#eq-3|(3)]] to equation [[#eq-32|(32)]] yields
+Applying the fourth-order operator [[#eq-3|(3)]] to equation [[#eq-32|(32)]] for <math display="inline">i=2, \dots , N</math>, yields
 <span id="eq-33"></span>
@@ Line 585: / Line 568: @@
 {| 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="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).  </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
+By employing the IFD method [[#eq-2|(2)]] and approximation [[#eq-12|(12)]], equation [[#eq-33|(33)]] can be approximated for regular points as follows
 <span id="eq-34"></span>
@@ Line 598: / Line 581: @@
 {| 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="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\}</math>
+|-
+| style="text-align: center;" | <math>  +  O(\Delta t^2) + O(h^4), \qquad 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
+By using the IFD-IIM [[#eq-6|(6)]] and approximation [[#eq-13|(13)]], the implicit scheme for irregular points is given by
 <span id="eq-35"></span>
@@ Line 611: / Line 596: @@
 {| 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="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 </math>
+|-
+| style="text-align: center;" | <math>  +  O(\Delta t^2) + O(h^3), \qquad i= I,I+1,  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (35)
@@ Line 637: / Line 624: @@
 {| 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>(b-r)U_{i-1}^{n+1}+(1-2b+2r)U_i^{n+1}+(b-r)U_{i+1}^{n+1} </math>
+|-
+| style="text-align: center;" | <math> = (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>
@@ Line 648: / Line 637: @@
 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 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 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]]].
+'''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==
@@ Line 658: / Line 647: @@
 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.
+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. Finally, Example 3 presents a discontinuous problem with multiple interface points. A Matlab code of these examples can be download at https://github.com/CIMATMerida/IFD-IIM
 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
@@ Line 668: / Line 657: @@
 {| 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="text-align: center;" | <math>\left\|e \right\|_\infty = \max _{1\leq i\leq N+1} \left|u_{i}-U_{i} \right|,   \quad \mathrm{and} \quad  \left\|e\right\|_2 = \left(\sum _{i=1}^{N+1} (u_{i}-U_{i})^2\Delta x\right)^{1/2}, </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (38)
@@ Line 686: / Line 675: @@
 |}
-where <math display="inline">N_1</math> and <math display="inline">N_2</math> indicates the different number of sub-intervals.
+where <math display="inline">N_1</math> and <math display="inline">N_2</math> indicates the different number of sub-intervals of the corresponding norm.
 ===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
+Example 1 considers the one-dimensional Poisson problem
 <span id="eq-40"></span>
@@ Line 698: / Line 687: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>u(x) = e^{-4\pi (x-1/4)^2}.  </math>
+| style="text-align: center;" | <math>\begin{array}{rcl}u_{xx}(x) &=& f(x),\quad x\in (0, 1), \\[1mm] u(0) &=& e^{-\pi /4}, \\  u(1) &=& e^{-9\pi /4}, \end{array}  </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.
+where the right-hand side in [[#eq-40|(40)]] is given by
-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]]].
+<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>f(x)= 4\pi (16\pi x^2-8\pi x+\pi{-2)}e^{-4\pi (x-1/4)^2}.  </math>
+|}
+|}
+The exact solution of problem [[#eq-40|(40)]] is an infinitely smooth and differentiable function given by the following expression
-{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
+<span id="eq-41"></span>
-|+ style="font-size: 75%;" |<span id='table-1'></span>'''Table. 1''' Convergence analysis of Example 1 testing a Poisson equation with smooth solution.
+{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|      0.95 @c cccc cccc
+|
-| style="text-align: right;" |
+{| style="text-align: left; margin:auto;width: 100%;"
-|
 |-
-| (rr)2-5  (rr)6-9 <math display="inline">N</math>
+| style="text-align: center;" | <math>u(x) = e^{-4\pi (x-1/4)^2}.  </math>
-| <math display="inline">L_\infty </math>-norm
+|}
-|  Order
+| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
-| <math display="inline">L_2</math>-norm
+|}
-|  Order
-| <math display="inline">L_\infty </math>-norm
+Note that the function <math display="inline">f</math> and Dirichlet boundary conditions are obtained directly from [[#eq-41|(41)]]. 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> are equal to zero.
-|  Order
-| <math display="inline">L_2</math>-norm
+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]]].
-|  Order
+{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:99%;"
+|+ style="font-size: 75%;" |<span id='table-1'></span>'''Table. 1''' Convergence analysis of Example 1 testing a Poisson equation with smooth solution.
+!
+! colspan="4" |<math display="inline">b=0</math>
+! colspan="4" |<math display="inline">b=1/12</math>
+|-
+! <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
@@ Line 785: / Line 798: @@
 |  '''4.06'''
 |-
-|
 |}
 ===5.2 Example 2. Poison equation with a single interface===
-For Example 2, we show the method's capability by solving a single interface problem located at <math display="inline">x = x_\alpha </math>. The exact solution is given by the function
+For Example 2, we show the method's capability by solving a single interface problem located at <math display="inline">x = x_\alpha </math>. Here, we consider the one-dimensional Poisson problem
-<span id="eq-41"></span>
+<span id="eq-42"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 798: / Line 810: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>u(x) = \begin{cases}\sin (\pi x)& x\leq x_\alpha ,\\ \cos (\pi x)& x   > x_\alpha{.} \end{cases}  </math>
+| style="text-align: center;" | <math>\begin{array}{rcl}u_{xx}(x) &=& f(x),\quad x\in (0,x_\alpha ) \cup (x_\alpha , 1), \\[1mm] u(0) &=& 0 \\[1mm] u(1) &=& -1, \\[1mm] \left[u\right]  & = & \cos (\pi x_\alpha ) - \sin (\pi x_\alpha ),\\[1mm] \left[u_x\right]& = &-\pi (\sin (\pi x_\alpha ) + \cos (\pi x_\alpha )),  \end{array}  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
 |}
-The right-hand side, <math display="inline">f</math>, is obtained directly from equation [[#eq-41|(41)]]. We test two different points: <math display="inline">x_\alpha=0.40</math> and <math display="inline">x_\alpha=0.63</math>.  For the first case, we always have <math display="inline">h_R/h=1</math> for <math display="inline">N=10\times{2}^{n}</math> (<math display="inline">n=0,1,2,\dots </math>); thus the interface is always located at one grid point of that resolution. In general, for <math display="inline">x_\alpha=0.63</math>, we have different <math display="inline">h_R/h</math> values for different <math display="inline">N</math> numbers. Figure [[#img-2|2]] shows the numerical and exact solution when the interface is located at these two values using <math display="inline">N = 40</math>. As expected, the exact solution is accurately recovered for both cases.
+where the right-hand side in [[#eq-42|(42)]] is
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>f(x) = \begin{cases}-\pi ^2\sin (\pi x), & x < x_\alpha ,\\ -\pi ^2\cos (\pi x), & x > x_\alpha{.} \end{cases} </math>
+|}
+|}
+The exact solution of problem [[#eq-42|(42)]] is a discontinuous function given by the expression
+<span id="eq-43"></span>
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>u(x) = \begin{cases}\sin (\pi x), & x < x_\alpha ,\\ \cos (\pi x), & x > x_\alpha{.} \end{cases}  </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
+|}
+Note that the right-hand side, Dirichlet boundary conditions, and principal jump conditions are obtained directly from equation [[#eq-43|(43)]].
+We test two different points: <math display="inline">x_\alpha=0.4</math> and <math display="inline">x_\alpha=0.63</math>.  For the case <math display="inline">x_\alpha=0.4</math>, we always have <math display="inline">h_R/h=1</math> for <math display="inline">N=10\times{2}^{n}</math> (<math display="inline">n=0,1,2,\dots </math>); thus the interface is always located at one grid point of that resolution. In general, for <math display="inline">x_\alpha=0.63</math>, we have different <math display="inline">h_R/h</math> values for different <math display="inline">N</math> numbers. Fig.&nbsp;[[#img-2|2]] shows the numerical and exact solution when the interface is located at these two values using <math display="inline">N = 40</math>. As expected, the exact solution is accurately recovered for both cases.
 <div id='img-2'></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_Poisson_Ex2_Sol.png|594px|Numerical and exact solution of Example 2 using N = 40 using (a) x_α= 0.4, and (b) x_α= 0.63.]]
+|[[Image:Review_499745022132-Fig02_Poisson_Ex2_Sol.png|594px|Numerical and exact solution of Example 2 using N = 40 using (a) x_α= 0.4, and (b) x_α= 0.63. Here, the interface is located at one grid point using N=40 and x_α=0.4.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 2:''' Numerical and exact solution of Example 2 using <math>N = 40</math> using (a) <math>x_\alpha = 0.4</math>, and (b) <math>x_\alpha = 0.63</math>.
+| colspan="1" | '''Figure 2:''' Numerical and exact solution of Example 2 using <math>N = 40</math> using (a) <math>x_\alpha = 0.4</math>, and (b) <math>x_\alpha = 0.63</math>. Here, the interface is located at one grid point using <math>N=40</math> and <math>x_\alpha=0.4</math>.
 |}
-Table [[#table-2|2]] shows the convergence analysis for Example 2 for the two <math display="inline">x_\alpha </math> values. Observe that a second-order method is recovered for <math display="inline">b = 0</math> and a fourth-order method for <math display="inline">b=1/12</math>. As expected, the IFD-IIM numerical order does not depend on the location of the interface. However, the magnitude of the error may present minor variations due to the interface position. Figure [[#img-3|3]] 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 ,320</math>.
+Table [[#table-2|2]] shows the convergence analysis for Example 2 for the two <math display="inline">x_\alpha </math> values. Observe that a second-order method is recovered for <math display="inline">b = 0</math> and a fourth-order method for <math display="inline">b=1/12</math>. As expected, the IFD-IIM numerical order does not depend on the location of the interface. However, the magnitude of the error may present minor variations due to the interface position. Fig.&nbsp;[[#img-3|3]] 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 ,320</math>. The solid black line for this and the following figures shows the numerical order calculated by the regression-line slope based on a least squares method of the <math display="inline">L_{\infty }</math>-norm.
-{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
+{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:99%;"
 |+ style="font-size: 75%;" |<span id='table-2'></span>'''Table. 2''' Convergence analysis of Example 2 using the IFD-IIM.
+!
+! colspan="4" |<math display="inline">b=0</math>
+! colspan="4" |<math display="inline">b=1/12</math>
 |-
-|     0.95 @c cccc cccc
+! <math display="inline">N</math>
-| style="text-align: right;" |
+! <math display="inline">L_\infty </math>-norm
-|
+!  Order
-|-
+! <math display="inline">L_2</math>-norm
-| (rrrr)2-5  (rrrr)6-9
+!  Order
-|
+! <math display="inline">L_\infty </math>-norm
-|
+!  Order
-|
+! <math display="inline">L_2</math>-norm
-|
+!  Order
-|-
-| (rr)2-3  (rr)4-5 (rr)6-7 (rr)8-9 <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
-| <math display="inline">L_\infty </math>-norm
-|  Order
 |-
 | 10
@@ Line 898: / Line 929: @@
 |  '''3.99'''
 |-
-|
 |}
 <div id='img-3'></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_Poisson_Ex2_OrderDispersion.png|594px|Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) x_α= 0.40, and (b) x_α= 0.63.]]
+|[[Image:Review_499745022132-Fig03_Poisson_Ex2_Order4th.png|594px|Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) x_α= 0.40, and (b) x_α= 0.63. The solid black line shows the numerical order 4.00 calculated by the regression-line slope based on a least squares method of the L<sub>∞</sub>-norm.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 3:''' Convergence analysis of Example 2 for <math>N = 10,\dots ,320</math> using (a) <math>x_\alpha = 0.40</math>, and (b) <math>x_\alpha = 0.63</math>.
+| colspan="1" | '''Figure 3:''' Convergence analysis of Example 2 for <math>N = 10,\dots ,320</math> using (a) <math>x_\alpha = 0.40</math>, and (b) <math>x_\alpha = 0.63</math>. The solid black line shows the numerical order <math>4.00</math> calculated by the regression-line slope based on a least squares method of the <math>L_{\infty }</math>-norm.
 |}
-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-4|4]] 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}]=[f_{xxx}]</math>, we observe that the error oscillation decreases in comparison with Fig.&nbsp;[[#img-3|3]] results. It is expected because 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.
+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-4|4]] shows that if we add additional jumps of high-order derivatives into <math display="inline">\mathcal{C}</math>, such as <math display="inline">[u_{xxxxx}]=[f_{xxx}]</math>, we observe that the error oscillation decreases in comparison with Fig.&nbsp;[[#img-3|3]] results. It is expected because 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.
 <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_Uh Zapata_643844164-Fig_Poisson_Ex2_OrderCloud5th.png|594px|Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) α= 0.40, and (b) α= 0.63. The contribution term includes jumps up to fifth-order ([uₓₓₓₓₓ]=[fₓₓₓ]).]]
+|[[Image:Review_499745022132-Fig04_Poisson_Ex2_Order5th.png|594px|Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) α= 0.40, and (b) α= 0.63. The contribution includes jumps up to fifth-order ([uₓₓₓₓₓ]=[fₓₓₓ]).]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 4:''' Convergence analysis of Example 2 for <math>N = 10,\dots ,320</math> using (a) <math>\alpha = 0.40</math>, and (b) <math>\alpha = 0.63</math>. The contribution term includes jumps up to fifth-order (<math>[u_{xxxxx}]=[f_{xxx}]</math>).
+| colspan="1" | '''Figure 4:''' Convergence analysis of Example 2 for <math>N = 10,\dots ,320</math> using (a) <math>\alpha = 0.40</math>, and (b) <math>\alpha = 0.63</math>. The contribution includes jumps up to fifth-order (<math>[u_{xxxxx}]=[f_{xxx}]</math>).
 |}
-Now we study the effect of removing jumps in <math display="inline">\mathcal{C}</math>. Let us consider a contribution term up to the third derivative jump by dropping <math display="inline">[u_{xxxx}]=[f_{xx}]</math>. Thus, the numerical approximation is only third-order accurate, as shown in Fig.&nbsp;[[#img-5|5]] for different interface values.  Note that error oscillations may have a clear pattern or a random distribution depending on the interface location. In general, the error magnitude perturbations are related to the variations coming from <math display="inline">h_R/h</math> and <math display="inline">h_L/h</math> values in the <math display="inline">\mathcal{C}</math> formula [[#eq-26|(26)]].
+Now we study the effect of removing jumps in <math display="inline">\mathcal{C}</math>. Let us consider a contribution term up to the third derivative jump by dropping <math display="inline">[u_{xxxx}]=[f_{xx}]</math>. Thus, the numerical approximation is only third-order accurate, as shown in Fig.&nbsp;[[#img-5|5]] for different interface values.  Note that error oscillations may have a clear pattern or a scattered distribution depending on the interface location. In general, the error magnitude perturbations are related to the variations coming from <math display="inline">h_R/h</math> and <math display="inline">h_L/h</math> values in the <math display="inline">\mathcal{C}</math> formula [[#eq-26|(26)]].
 <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_Uh Zapata_643844164-Fig_Poisson_Ex2_OrderCloud3rd.png|594px|Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) α= 0.40, and (b) α= 0.63. The contribution term includes jumps up to third-order ([uₓₓₓ]=[fₓ]).]]
+|[[Image:Review_499745022132-Fig05_Poisson_Ex2_Order3rd.png|594px|Convergence analysis of Example 2 for N = 10,\dots ,320 using (a) α= 0.40, and (b) α= 0.63. The contribution includes jumps up to third-order ([uₓₓₓ]=[fₓ]).]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 5:''' Convergence analysis of Example 2 for <math>N = 10,\dots ,320</math> using (a) <math>\alpha = 0.40</math>, and (b) <math>\alpha = 0.63</math>. The contribution term includes jumps up to third-order (<math>[u_{xxx}]=[f_{x}]</math>).
+| colspan="1" | '''Figure 5:''' Convergence analysis of Example 2 for <math>N = 10,\dots ,320</math> using (a) <math>\alpha = 0.40</math>, and (b) <math>\alpha = 0.63</math>. The contribution includes jumps up to third-order (<math>[u_{xxx}]=[f_{x}]</math>).
 |}
-Note that there are some points in Fig.&nbsp;[[#img-5|5]], marked with circles, that are close to the fourth-order line. Those circled markers correspond to the values with <math display="inline">h_R/h=1</math>: <math display="inline">N = 5,10,15,\dots ,315</math> for <math display="inline">x_\alpha=0.40</math>, and <math display="inline">N = 100, 200, 300</math> for <math display="inline">x_\alpha=0.63</math>.  Both set of points satisfy that the interface <math display="inline">x_\alpha </math> is located at a grid point <math display="inline">x_I</math>. In the case of <math display="inline">x_\alpha=0.4</math>, we observe that the global order (black points) is close to <math display="inline">3.21</math>; meanwhile, the circled points order is four. Similar behavior is obtained for <math display="inline">x_\alpha=0.63</math>. Thus, for a given mesh resolution <math display="inline">N</math>, the IFD-IIM is still fourth-order accurate for <math display="inline">h_R/h = 1</math>, as  proven theoretically in Section [[#3.2 IFD-IIM and principal jump conditions|3.2]].
+Note that there are some points in Fig.&nbsp;[[#img-5|5]], marked with circles, that are close to the fourth-order line. Those circled markers correspond to the values with <math display="inline">h_R/h=1</math>: (a) <math display="inline">N = 5,10,15,\dots ,315</math> for <math display="inline">x_\alpha=0.40</math>, and (b) <math display="inline">N = 100, 200, 300</math> for <math display="inline">x_\alpha=0.63</math>.  Both set of points satisfy that the interface <math display="inline">x_\alpha </math> is located at a grid point <math display="inline">x_I</math>. In the case of <math display="inline">x_\alpha=0.4</math>, we observe that the global order (black points) is close to <math display="inline">3.21</math>; meanwhile, the circled points order is four. Similar behavior is obtained for <math display="inline">x_\alpha=0.63</math>. Thus, for a given mesh resolution <math display="inline">N</math>, the IFD-IIM is still fourth-order accurate for <math display="inline">h_R/h = 1</math>, as  proven theoretically in Section [[#4.2 IFD-IIM and principal jump conditions|4.2]].
 ===5.3 Example 3. Poison equation with multiple interface points===
@@ Line 931: / Line 961: @@
 Example 3 investigates the numerical scheme capability to solve a multiple interface problem. We only focus on two interface points located at <math display="inline">x = x_{\alpha _1}</math> and <math display="inline">x = x_{\alpha _2}</math>. However, the methodology could be applied for multiple interfaces by doing minor modifications in the implementation. For this problem, the exact solution of the Poisson problem is given by
-<span id="eq-42"></span>
+<span id="eq-44"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 937: / Line 967: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>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="text-align: center;" | <math>u(x) = \begin{cases}\sin (2\pi x), & x < x_{\alpha _1},\\ \cos (2\pi x), & x_{\alpha _1}< x < x_{\alpha _2}, \\ \sin (2\pi x), & x_{\alpha _2} < x. \end{cases}  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
 |}
-The right-hand function is obtained directly from the exact solution [[#eq-42|(42)]].  We consider the same computational domain and grid resolution as previous 1D examples. Fig.&nbsp;[[#img-6|6]] presents the analytical and numerical solution when the interface is located at <math display="inline">x_{\alpha _1}= 0.3</math> and <math display="inline">x_{\alpha _2} = 0.7</math> using <math display="inline">N=40</math>. This figure also shows the error analysis. As expected, the IFD-IIM is a fourth-order accurate method.
+The right-hand function, Dirichlet boundary conditions, and principal jump conditions are obtained directly from the exact solution [[#eq-44|(44)]].  We consider the same computational domain and grid resolution as previous 1D examples. Fig.&nbsp;[[#img-6|6]] presents the analytical and numerical solution when the interface is located at <math display="inline">x_{\alpha _1}= 0.3</math> and <math display="inline">x_{\alpha _2} = 0.7</math> using <math display="inline">N=40</math>. This figure also shows the error analysis. As expected, the IFD-IIM is a fourth-order accurate method.
 <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_Uh Zapata_643844164-Fig_Poisson_Ex3.png|594px|(a) Numerical and exact solution of Example 3 with multiple interfaces using N = 40, and (b) convergence error analysis using different grid resolutions.]]
+|[[Image:Review_499745022132-Fig06_Poisson_Ex3.png|594px|(a) Numerical and exact solution of Example 3 with multiple interfaces using N = 40, and (b) convergence error analysis using  grid resolutions of N = 10, 20, 40, 80 and 160.]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 6:''' (a) Numerical and exact solution of Example 3 with multiple interfaces using <math>N = 40</math>, and (b) convergence error analysis using different grid resolutions.
+| colspan="1" | '''Figure 6:''' (a) Numerical and exact solution of Example 3 with multiple interfaces using <math>N = 40</math>, and (b) convergence error analysis using  grid resolutions of <math>N = 10</math>, 20, 40, 80 and 160.
 |}
 ==6 Numerical results for the Heat equation==
-This section tests the IFD-IIM for the Heat conduction equation in different scenarios. Example 4 verifies the fourth-order implicit method for smooth solutions. Example 5 studies a Heat equation with a discontinuous solution in a single interface point and no source term. Finally, Example 6 presents a general discontinuous problem. For the all these examples, the computational domain is the same interval <math display="inline">[0,1]</math> and space step, <math display="inline">h</math>, used for the Poisson examples. Here, the final time is at <math display="inline">T=0.5</math> and the time step is <math display="inline">\Delta t = \frac{1}{4}h^2</math>. The error and estimated order of accuracy are reported using  [[#eq-38|(38)]] and [[#eq-39|(39)]], respectively, at the final step.
+This section tests the IFD-IIM for the Heat conduction equation in different scenarios. Example 4 verifies the fourth-order implicit method for smooth solutions in space. Example 5 studies a Heat equation with a discontinuous solution in a single interface point and no source term. Finally, Example 6 presents a general discontinuous problem. For the all these examples, the computational domain is the same interval <math display="inline">[0,1]</math> and space step, <math display="inline">h</math>, used for the Poisson examples. Here, the final time is at <math display="inline">T=0.5</math> and the time step is <math display="inline">\Delta t = \frac{1}{4}h^2</math>. The error and estimated order of accuracy are reported using  [[#eq-38|(38)]] and [[#eq-39|(39)]], respectively, at the final step approximation.
 ===6.1 Example 4. Heat equation with smooth solution===
+This example considers the heat conduction problem
+<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>\begin{array}{rcll}u_t(x,t) &=& u_{xx}(x,t)-f(x,t), & (x,t)\in (0, 1)\times (0,\frac{1}{2}], \\[2mm] u(x,0) &=& \sin \left(\frac{3}{2}\pi x\right),&  x \in (0, 1), \\[2mm] u(0,t) &=& 0, &  t\in [0,\frac{1}{2}],\\[2mm] u(1,t) &=& -e^{-t}, &  t\in [0,\frac{1}{2}], \end{array}  </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
+|}
+where <math display="inline">f</math> in [[#eq-45|(45)]] is given by
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>f(x,t)= \left(1-\frac{9}{4}\pi ^2\right)\sin \left(\frac{3}{2}\pi x\right)e^{-t}. </math>
+|}
+|}
 This example is constructed so that the exact solution is
-<span id="eq-43"></span>
+<span id="eq-46"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 966: / Line 1,020: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>u(x,t) = \sin \left(\frac{3}{2}\pi x\right)\exp (-t),  </math>
+| style="text-align: center;" | <math>u(x,t) = \sin \left(\frac{3}{2}\pi x\right)e^{-t}.  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
 |}
-where the source term <math display="inline">f</math> is derived from [[#eq-27|(27)]] and [[#eq-43|(43)]]. The initial and boundary conditions are also obtained from the exact solution. The convergence analysis of this example is presented in Table [[#table-3|3]] with the final time being <math display="inline">T = 0.5</math>. As expected, the high-order methods reach their corresponding order of accuracy.
+Note that the source term, <math display="inline">f</math>, initial condition and boundary conditions are obtained from the exact solution [[#eq-46|(46)]]. The convergence analysis of this example is presented in Table [[#table-3|3]] with the final time being <math display="inline">T = 0.5</math>. As expected, the high-order methods reach their corresponding order of accuracy.
-{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
+{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:99%;"
 |+ style="font-size: 75%;" |<span id='table-3'></span>'''Table. 3''' Convergence analysis of Example 4 testing a heat equation with smooth solution.
+!
+! colspan="4" |<math display="inline">b=0</math>
+! colspan="4" |<math display="inline">b=1/12</math>
 |-
-|      0.95 @c cccc cccc
+! <math display="inline">N</math>
-| style="text-align: right;" |
+! <math display="inline">L_\infty </math>-norm
-|
+!  Order
-|-
+! <math display="inline">L_2</math>-norm
-| (rr)2-5  (rr)6-9 <math display="inline">N</math>
+!  Order
-| <math display="inline">L_\infty </math>-norm
+! <math display="inline">L_\infty </math>-norm
-|  Order
+!  Order
-| <math display="inline">L_2</math>-norm
+! <math display="inline">L_2</math>-norm
-|  Order
+!  Order
-| <math display="inline">L_\infty </math>-norm
-|  Order
-| <math display="inline">L_2</math>-norm
-|  Order
 |-
 | 10
@@ Line 1,051: / Line 1,104: @@
 |  '''4.00'''
 |-
-|
 |}
 ===6.2 Example 5. Heat equation with discontinuous solution===
-In this example, the exact solution <math display="inline">u</math> is defined as
+In this example, the heat equation with discontinuous solution is given by
-<span id="eq-44"></span>
+<span id="eq-47"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 1,064: / Line 1,116: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>u(x) = \begin{cases}\sin (x)\exp (-t) & x\leq x_\alpha ,\\ \cos (x)\exp (-t) & x   > x_\alpha{.} \end{cases}  </math>
+| style="text-align: center;" | <math>\begin{array}{rcll}u_t(x,t) &=& u_{xx}(x,t), & (x,t)\in (0,x_\alpha ) \cup (x_\alpha , 1) \times (0,\frac{1}{2}], \\[2mm] u(x,0) &=&\begin{cases}\sin (x), & x < x_\alpha ,\\ \cos (x), & x > x_\alpha ,\\ \end{cases}  & x \in (0,x_\alpha ) \cup (x_\alpha , 1), \\ \\[-3mm] u(0,t) &=& 0, & t\in [0,\frac{1}{2}],\\[2mm] u(1,t) &=& \cos (1)e^{-t}, & t\in [0,\frac{1}{2}],\\[2mm] \left[u\right](t)  &=&     \left(\cos (x_\alpha )-\sin (x_\alpha )\right)e^{-t}, &  t\in [0,\frac{1}{2}],\\[2mm] \left[u_{x}\right](t)  &=& -\left(\sin (x_\alpha )+\cos (x_\alpha )\right)e^{-t}, & t\in [0,\frac{1}{2}].  \end{array}  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
 |}
-The source term <math display="inline">f</math> of this problem is <math display="inline">f\equiv 0</math>. In this example, both the function <math display="inline">u</math> and their derivatives have jumps, and these jumps vary in time. The boundary condition, initial condition, and all jumps are specified by <math display="inline">u</math>.
+The exact solution <math display="inline">u</math> is defined as
+<span id="eq-48"></span>
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>u(x,t) = \begin{cases}\sin (x)e^{-t}, & x < x_\alpha ,\\ \cos (x)e^{-t}, & x > x_\alpha{.} \end{cases}  </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
+|}
+Similar to previous examples, we use <math display="inline">x_\alpha=0.4</math> and <math display="inline">x_\alpha=0.63</math>. Note that the source term <math display="inline">f</math> of this problem is <math display="inline">f\equiv 0</math>. In this example, both the function <math display="inline">u</math> and their derivatives have jumps, and these jumps vary in time. The boundary condition, initial condition, and all jumps are specified by <math display="inline">u</math>.
+The figure of numerical solution and absolute error plot using <math display="inline">N=40</math> for different time stages are shown in Fig.&nbsp;[[#img-7|7]]. On the other hand, the one-dimensional results at <math display="inline">T=0.5</math> are presented in Figs.&nbsp;[[#img-8|8]] and [[#img-9|9]]. As expected, we can accurately recover the exact solution using the proposed method for <math display="inline">N=40</math> and considering different interface values.  Additionally, we provide the absolute error for both cases, noting that the maximum error is found in close proximity to the point <math display="inline">x_\alpha </math>.  This behavior is also expected since the local truncation error for irregular points is less accurate than for regular points. Finally, the convergence analysis using <math display="inline">N=10</math>, 20, 40, 80 and 160 confirms the fourth-order global convergence of the proposed method.
-The figure of numerical solution and absolute error plot using <math display="inline">N=40</math> for different time stages are shown in Fig.&nbsp;[[#img-7|7]]. On the other hand, the one-dimensional results at <math display="inline">t=0.5</math> are presented in Figs.&nbsp;[[#img-8|8]] and [[#img-9|9]]. More details of the grid refinement analysis at <math display="inline">t=0.5</math> is presented in Table [[#table-4|4]] for two different interface point locations. As expected, a fourth-order method is obtained for both interfaces. We remark that the variation of the errors using <math display="inline">x_\alpha=0.63</math> is because the different <math display="inline">h_R/h</math> values resulting from the discretization.
+More details of the grid refinement analysis at <math display="inline">T=0.5</math> is presented in Table [[#table-4|4]] for two different interface point locations. As expected, a fourth-order method is obtained for both interfaces. We remark that the variation of the errors using <math display="inline">x_\alpha=0.63</math> is because the different <math display="inline">h_R/h</math> values resulting from the discretization.
 <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_Uh Zapata_643844164-Fig_Heat_Ex5_3D.png|594px|Numerical solution and absolute error of Example 5 using N = 40, x_α=0.4 for tퟄ[0,0.5].]]
+|[[Image:Review_499745022132-Fig07_Heat_Ex5_3D.png|594px|Numerical solution and absolute error of Example 5 using N = 40, x_α=0.4 for tퟄ[0,0.5].]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 7:''' Numerical solution and absolute error of Example 5 using <math>N = 40</math>, <math>x_\alpha=0.4</math> for <math>t\in [0,0.5]</math>.
@@ Line 1,083: / Line 1,150: @@
 {| 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_Heat_SolEx5_0p40.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.4 of Example 5.]]
+|[[Image:Review_499745022132-Fig08_Heat_Ex5_0p40.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.4 of Example 5.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 8:''' Numerical results at final time simulation <math>T=0.5</math> and <math>x_\alpha=0.4</math> of Example 5.
@@ Line 1,090: / Line 1,157: @@
 {| 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_Heat_SolEx5_0p63.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.63 of Example 5.]]
+|[[Image:Review_499745022132-Fig09_Heat_Ex5_0p63.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.63 of Example 5.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 9:''' Numerical results at final time simulation <math>T=0.5</math> and <math>x_\alpha=0.63</math> of Example 5.
 |}
-{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
+{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:99%;"
 |+ style="font-size: 75%;" |<span id='table-4'></span>'''Table. 4''' Convergence analysis of Example 5 using the IFD-IIM for the heat equation with a discontinuous solution.
+!
+! colspan="4" |<math display="inline">x_\alpha=0.4</math>
+! colspan="4" |<math display="inline">x_\alpha=0.63</math>
 |-
-|     0.95 @c cccc cccc
+! <math display="inline">N</math>
-| style="text-align: right;" |
+! <math display="inline">L_\infty </math>-norm
-|
+!  Order
-|-
+! <math display="inline">L_2</math>-norm
-| (rr)2-5  (rr)6-9 <math display="inline">N</math>
+!  Order
-| <math display="inline">L_\infty </math>-norm
+! <math display="inline">L_\infty </math>-norm
-|  Order
+!  Order
-| <math display="inline">L_2</math>-norm
+! <math display="inline">L_2</math>-norm
-|  Order
+!  Order
-| <math display="inline">L_\infty </math>-norm
-|  Order
-| <math display="inline">L_2</math>-norm
-|  Order
 |-
 | 10
@@ Line 1,172: / Line 1,238: @@
 |  '''3.96'''
 |-
-|
 |}
 From Table&nbsp;[[#table-5|5]], we can find the results when the time step is chosen as <math display="inline">\Delta t = h</math>. As expected, the proposed IFD-IIM is stable and has two-order convergence either we take a second- or fourth-order approximation in space. We remark that this is a property of the selected time integration method. For instance, similar findings as Table&nbsp;[[#table-5|5]] can be obtained for smooth solutions.
-{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
+{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:99%;"
 |+ style="font-size: 75%;" |<span id='table-5'></span>'''Table. 5''' Convergence analysis of Example 5 testing the heat equation with <math>\Delta t=h</math>.
+!
+! colspan="4" |<math display="inline">b=0</math>
+! colspan="4" |<math display="inline">b=1/12</math>
 |-
-|      0.95 @c cccc cccc
+! <math display="inline">N</math>
-|
+! <math display="inline">L_\infty </math>-norm
-|
+!  Order
-|-
+! <math display="inline">L_2</math>-norm
-| (rr)2-5  (rr)6-9 <math display="inline">N</math>
+!  Order
-| <math display="inline">L_\infty </math>-norm
+! <math display="inline">L_\infty </math>-norm
-|  Order
+!  Order
-| <math display="inline">L_2</math>-norm
+! <math display="inline">L_2</math>-norm
-|  Order
+!  Order
-| <math display="inline">L_\infty </math>-norm
-|  Order
-| <math display="inline">L_2</math>-norm
-|  Order
 |-
 | 10
@@ Line 1,254: / Line 1,318: @@
 |  '''1.95'''
 |-
-|
 |}
@@ Line 1,261: / Line 1,324: @@
 Finally, Example 6 investigates the capability to solve a heat interface problem with a general discontinuous solution. For this problem, we slightly modified the previous example such that the exact solution of the heat conduction problem is given by
-<span id="eq-45"></span>
+<span id="eq-49"></span>
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 1,267: / Line 1,330: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>u(x) = \begin{cases}\sin (2\pi x)\exp (-t) & x\leq x_{\alpha },\\ \cos (2\pi x)\exp (-t) & x > x_{\alpha }. \end{cases}  </math>
+| style="text-align: center;" | <math>u(x) = \begin{cases}\sin (2\pi x)e^{-t}, & x < x_{\alpha },\\ \cos (2\pi x)e^{-t}, & x > x_{\alpha }. \end{cases}  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
 |}
-Here, the source term <math display="inline">f</math> is a general function that varies both in time and space. It is directly derived from equation [[#eq-27|(27)]] as well as the exact solution [[#eq-45|(45)]]. Furthermore, the exact solution is utilized to specify the boundary condition, initial condition, and all jumps contributions.
+Here, the source term <math display="inline">f</math> is a general function that varies both in time and space. It is directly derived from equation [[#eq-27|(27)]] as well as the exact solution [[#eq-49|(49)]]. Furthermore, the exact solution is utilized to specify the boundary condition, initial condition, and all jumps contributions.
-In Fig.&nbsp;[[#img-10|10]], we depict the temporal evolution of the numerical solution and absolute errors using parameters <math display="inline">N=40</math> and <math display="inline">x_\alpha=0.4</math>. It is noteworthy that the behavior of the solution differs from the previous example. More in-depth analysis of the results at <math display="inline">T = 0.5</math> are illustrated in Figs.&nbsp;[[#img-11|11]] and [[#img-12|12]], corresponding to values of <math display="inline">x_\alpha=0.4</math> and <math display="inline">x_\alpha=0.63</math>, respectively. Grid refinement analyses are also provided in these figures. As anticipated, the IFD-IIM method demonstrates fourth-order accuracy, a validation that is further detailed in Table [[#table-6|6]].
+In Fig.&nbsp;[[#img-10|10]], we depict the temporal evolution of the numerical solution and absolute errors using parameters <math display="inline">N=40</math> and <math display="inline">x_\alpha=0.4</math>. It is noteworthy that the behavior of the solution differs from the previous example. More in-depth analysis of the results at <math display="inline">T = 0.5</math> are illustrated in Figs.&nbsp;[[#img-11|11]] and [[#img-12|12]], corresponding to values of <math display="inline">x_\alpha=0.4</math> and <math display="inline">x_\alpha=0.63</math>, respectively. First, the numerical solution is accurately approximated using the proposed implicit method for <math display="inline">N=40</math>, regardless of the interface locations. It is worth noting that the maximum error for <math display="inline">x_\alpha=0.63</math> is concentrated near to this point; however, this is not observed for <math display="inline">x_\alpha=0.4</math>. This difference in behavior can be attributed to the interplay between the interface location (<math display="inline">h_R/h=1</math> and <math display="inline">x_\alpha=0.4</math>) and the non-zero right-hand side values, which contribute to the local truncation errors. Grid refinement analyses are also provided in these figures using <math display="inline">N=10</math>, 20, 40, 80 and 160. As anticipated, the IFD-IIM method demonstrates fourth-order accuracy, a validation that is further detailed in Table [[#table-6|6]].
 <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_Uh Zapata_643844164-Fig_Heat_Ex6_3D.png|594px|Numerical solution and absolute error of Example 6 using N = 40, x_α=0.4 for tퟄ[0,0.5].]]
+|[[Image:Review_499745022132-Fig10_Heat_Ex6_3D.png|594px|Numerical solution and absolute error of Example 6 using N = 40, x_α=0.4 for tퟄ[0,0.5].]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 10:''' Numerical solution and absolute error of Example 6 using <math>N = 40</math>, <math>x_\alpha=0.4</math> for <math>t\in [0,0.5]</math>.
@@ Line 1,286: / Line 1,349: @@
 {| 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_Heat_SolEx6_0p40.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.4 of Example 6. ]]
+|[[Image:Review_499745022132-Fig11_Heat_Ex6_0p40.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.4 of Example 6. ]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 11:''' Numerical results at final time simulation <math>T=0.5</math> and <math>x_\alpha=0.4</math> of Example 6.
@@ Line 1,293: / Line 1,356: @@
 {| 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_Heat_SolEx6_0p63.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.63 of Example 6.]]
+|[[Image:Review_499745022132-Fig12_Heat_Ex6_0p63.png|594px|Numerical results at final time simulation T=0.5 and x_α=0.63 of Example 6.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 12:''' Numerical results at final time simulation <math>T=0.5</math> and <math>x_\alpha=0.63</math> of Example 6.
 |}
-{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;"
+{|  class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:99%;"
 |+ style="font-size: 75%;" |<span id='table-6'></span>'''Table. 6''' Convergence analysis of Example 6 using the IFD-IIM for the heat equation with a general discontinuous solution.
+!
+! colspan="4" |<math display="inline">x_\alpha=0.4</math>
+! colspan="4" |<math display="inline">x_\alpha=0.63</math>
 |-
-|     0.95 @c cccc cccc
+! <math display="inline">N</math>
-| style="text-align: right;" |
+! <math display="inline">L_\infty </math>-norm
-|
+!  Order
-|-
+! <math display="inline">L_2</math>-norm
-| (rr)2-5  (rr)6-9 <math display="inline">N</math>
+!  Order
-| <math display="inline">L_\infty </math>-norm
+! <math display="inline">L_\infty </math>-norm
-|  Order
+!  Order
-| <math display="inline">L_2</math>-norm
+! <math display="inline">L_2</math>-norm
-|  Order
+!  Order
-| <math display="inline">L_\infty </math>-norm
-|  Order
-| <math display="inline">L_2</math>-norm
-|  Order
 |-
 | 10
@@ Line 1,375: / Line 1,437: @@
 |  '''3.95'''
 |-
-|
 |}
 ==7 Conclusions==
-This work serves as the general foundation for addressing high-order IIMs to solve interface problems, facilitating comprehensive investigations into this theoretical framework. Here, we present a fourth-order finite-difference scheme for approximating the second-order derivative of real-valued continuous and discontinuous functions. This method combines an implicit formulation with the immersed interface method. Our proposed scheme employs a three-point stencil, achieving different accuracy at regular and irregular points. To illustrate the effectiveness of the proposed IFD-IIM approach, we applied it to solve the one-dimensional Poisson equation and the heat conduction equation. The global accuracy of the fourth order was demonstrated using several numerical examples for both equations. Hence, this study establishes a general strategy for high-order immersed interface methods, enabling their application to elliptic and time-evolving problems in several dimensions. Additionally, the implicit procedure lends itself to developing of higher-order numerical schemes based on the IIM, including sixth-order methods.
+This work serves as the general foundation for addressing high-order IIMs to solve interface problems, facilitating comprehensive investigations into this theoretical framework. Here, we present a fourth-order finite-difference scheme for approximating the second-order derivative of real-valued continuous and discontinuous functions. This method combines an implicit formulation with the IIM. Our proposed scheme employs a three-point stencil, achieving different accuracy at regular and irregular points. To illustrate the effectiveness of the proposed IFD-IIM approach, we applied it to solve the one-dimensional Poisson equation and the heat conduction equation. The global accuracy of the fourth order was demonstrated using several numerical examples for both equations. Hence, this study establishes a general strategy for high-order IIMs, enabling their application to elliptic and time-evolving problems in several dimensions. Additionally, the implicit procedure lends itself to developing of higher-order numerical schemes based on the IIM, including sixth-order methods.
 ==Acknowledgments==
-This work was partially supported by CONAHCYT under the program ''Investigadoras e Investigadores por México''.
+This work was partially supported by CONAHCYT under the program ''Investigadoras e Investigadores por México''. We also thanks the facilities provided by the ''Facultad de Matemáticas (UADY)'' for the development of this work.
-==Appendix A: Matlab code to solve the 1D Poisson equation==
+===BIBLIOGRAPHY===
-This Appendix is focused on the Matlab code to solve Example 2, corresponding to a one-dimensional Poison equation with Dirichlet boundary conditions. For better exposition, the code was divided in four sections. In the first part, we provide the computational domain, the location of the interface, and a vector of different sub-divisions <code>Mvec</code>. Next, we present the main loop corresponding to the IFD-IIM implementation. The third part complement this loop solving the the linear system by the Thomas' Algorithm and the estimation of the order of accuracy using the exact solution. In the fourth section, we display the norm errors and order of accuracy. Finally, the program plots the solution for the last entree of <code>Mvec</code>. This program can be also download at https://github.com/CIMATMerida/IFD-IIM.
+<div id="cite-1"></div>
+'''[[#citeF-1|[1]]]''' Diersch, H-J. (1990) "Fletcher, CAJ, Computational Techniques for Fluid Dynamics. Vol. I: Fundamental and General Techniques. Vol. II: Specific Techniques for Different Flow Categories. Berlin etc., Springer-Verlag 1988. XIV, 409 pp., 183 figs./XI, 484 pp., 183 figs., DM 198, 00 as a Set. ISBN 3-540-18151-2/3-540-18759-6 (Springer Series in Computational Physics)". Wiley Online Library
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+<div id="cite-2"></div>
-|-
+'''[[#citeF-2|[2]]]''' Sethian, James A and others. (1999) "Level set methods and fast marching methods", Volume 98. Cambridge Cambridge UP 2
-<pre>
+<div id="cite-3"></div>
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+'''[[#citeF-3|[3]]]''' Li, Zhilin and Ito, Kazufumi. (2006) "The immersed interface method: numerical solutions of PDEs involving interfaces and irregular domains". SIAM
-%                                                                         %
-%                  HIGH-ORDER IMMERSED INTERFACE METHOD                   %
-%            H. EScamilla Puc,  R. Itza Balam & M. Uh Zapata              %
-%                                Sept 2023                                %
-%  It solves the one-dimensional Poisson equation:                        %
-%                           u_xx  = f                                     %
-%  knowing jump conditions: [u], [u_x], [f], [fx], and [fxx] at the       %
-%  interface and Dirichlet boundary conditions.                           %
-%                                                                         %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+<div id="cite-4"></div>
-% PROBLEM PARAMETERS & FUNCTIONS
+'''[[#citeF-4|[4]]]''' Javierre, E and Vuik, C and Vermolen, FJ and Van der Zwaag, S. (2006) "A comparison of numerical models for one-dimensional Stefan problems", Volume 192. Elsevier. Journal of Computational and Applied Mathematics 2 445&#8211;459
-clear
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
-% Discretization
-xI   = 0.0;               % Initial of the domain: (xI,xF)
-xF   = 1.0;               % Final   of the domain: (xI,xF)
-alf  = 0.4;               % Location of the interface
-Mvec = [10,20,40,80,160]; % Number of sub-divitions (vector)
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
-% Method
-b = 1/12;                 % b=1/12 (4th-order), b=0 (2nd-order)
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
-% Functions
-fun_uL    = @(x)         sin(pi*x);
-fun_uxL   = @(x)      pi*cos(pi*x);
-fun_fL    = @(x) -(pi^2)*sin(pi*x);
-fun_fxL   = @(x) -(pi^3)*cos(pi*x);
-fun_fxxL  = @(x)  (pi^4)*sin(pi*x);
-fun_uR    = @(x)         cos(pi*x);
+<div id="cite-5"></div>
-fun_uxR   = @(x)     -pi*sin(pi*x);
+'''[[#citeF-5|[5]]]''' Shi, Yu-e and Ray, Rajendra K and Nguyen, Kim Dan. (2013) "A projection method-based model with the exact C-property for shallow-water flows over dry and irregular bottom using unstructured finite-volume technique", Volume 76. Elsevier. Computers & Fluids 178&#8211;195
-fun_fR    = @(x) -(pi^2)*cos(pi*x);
-fun_fxR   = @(x)  (pi^3)*sin(pi*x);
-fun_fxxR  = @(x)  (pi^4)*cos(pi*x);
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+<div id="cite-6"></div>
-% WORKSPACE
+'''[[#citeF-6|[6]]]''' Uh, Miguel and Xu, Sheng. (2014) "The immersed interface method for simulating two-fluid flows", Volume 7. Cambridge University Press. Numerical Mathematics: Theory, Methods and Applications 4 447&#8211;472
-M = length(Mvec);
+<div id="cite-7"></div>
-hvec     = zeros(M,1);
+'''[[#citeF-7|[7]]]''' Li, Ming and Tang, Tao and Fornberg, Bengt. (1995) "A compact fourth-order finite difference scheme for the steady incompressible Navier-Stokes equations", Volume 20. Wiley Online Library. International Journal for Numerical Methods in Fluids 10 1137&#8211;1151
-NormE1   = zeros(M,1);
-NormE2   = zeros(M,1);
-OrderOr1 = zeros(M,1);
-OrderOr2 = zeros(M,1);
-for s=1:M
+<div id="cite-8"></div>
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
+'''[[#citeF-8|[8]]]''' Zhang, Jun. (2002) "Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization", Volume 179. Elsevier. Journal of Computational Physics 1 170&#8211;179
-% DISCRETIZATION
-n = Mvec(s);
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Points and step size
-x  = linspace(xI,xF,n+1)';
-h  = x(2)-x(1);
-h2 = h*h;
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Find I: the interval where alpha is
-for i=1:n
-if  (x(i)<=alf) && (alf<x(i+1))
-I = i;
-break;
-end
-end
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
-% MATRIX & RIGHT-HAND SIDE
-A1 = zeros(n,1);
-A2 = zeros(n+1,1);
-A3 = zeros(n,1);
-rhs= zeros(n+1,1);
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-f(1:I)     = fun_fL(x(1:I));
-f(I+1:n+1) = fun_fR(x(I+1:n+1));
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Boundary
-A2(1)    = 1/h2;
-rhs(1)   = fun_uL(x(1))/h2;
-A2(n+1)  = 1/h2;
-rhs(n+1) = fun_uR(x(n+1))/h2;
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Regular points: Equation (22) in Escamilla et al. 2023
-A1(1:n-1) =  1/h2;
-A3(2:n)   =  1/h2;
-A2(2:n)   = -2/h2;
-rhs(2:n)  =  b*f(3:n+1)+(1-2*b)*f(2:n)+b*f(1:n-1);
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Irregular points: Equation (26) in Escamilla et al. 2023
-hL   = x(I)   - alf;
-hR   = x(I+1) - alf;
-%&#8211;&#8211;-
-uJ   = fun_uR(alf)   - fun_uL(alf);    % [u]
-uxJ  = fun_uxR(alf)  - fun_uxL(alf);   % [ux]
-fJ   = fun_fR(alf)   - fun_fL(alf);    % [f]
-fxJ  = fun_fxR(alf)  - fun_fxL(alf);   % [fx]
-fxxJ = fun_fxxR(alf) - fun_fxxL(alf);  % [fxx]
-%&#8211;&#8211;-
-CI   =  (1/h2)*(uJ + hR*uxJ + 0.5*(hR^2)*fJ) ...
-- b*(fJ + hR*(1-2*hR^2/h2)*fxJ + 0.5*hR^2*(1-hR^2/h2)*fxxJ);
-CIp1 = -(1/h2)*(uJ + hL*uxJ + 0.5*(hL^2)*fJ) ...
-+ b*(fJ + hL*(1-2*hL^2/h2)*fxJ + 0.5*hL^2*(1-hL^2/h2)*fxxJ);
-%&#8211;&#8211;-
-rhs(I)   = rhs(I)   + CI;
-rhs(I+1) = rhs(I+1) + CIp1;
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
-% LINEAR SYSTEM SOLUTION BY THOMAS
-U = zeros(n+1,1);
-for i=1:n
-A2(i+1)=A2(i+1)-A3(i)*A1(i)/A2(i);
-rhs(i+1)=rhs(i+1)-rhs(i)*A1(i)/A2(i);
-end
-U(n+1)=rhs(n+1)/A2(n+1);
-for i=n:-1:1
-U(i)=(rhs(i)-U(i+1)*A3(i))/A2(i);
-end
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;
-% ERRORS & ORDER
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Exact solution
-Uexact(1:I)     = fun_uL(x(1:I));
-Uexact(I+1:n+1) = fun_uR(x(I+1:n+1));
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Norm errors
-Err       = abs(Uexact'-U);
-NormE1(s) = norm(Err,inf);
-NormE2(s) = sqrt(h*sum(Err.^2));
-%&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;-
-% Estimated order
-hvec(s) = h;
-if s&nbsp;=1
-div = log((hvec(s-1))/(hvec(s)));
-OrderOr1(s) = log(NormE1(s-1)/NormE1(s))/div;
-OrderOr2(s) = log(NormE2(s-1)/NormE2(s))/div;
-end
-end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+<div id="cite-9"></div>
-% DISPLAY & PLOT
+'''[[#citeF-9|[9]]]''' Nabavi, Majid and Siddiqui, MH Kamran and Dargahi, Javad. (2007) "A new 9-point sixth-order accurate compact finite-difference method for the Helmholtz equation", Volume 307. Elsevier. Journal of Sound and Vibration 3-5 972&#8211;982
-disp('   N   — Max-norm Order — L2-norm  Order')
+<div id="cite-10"></div>
-disp('&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;&#8211;')
+'''[[#citeF-10|[10]]]''' Wang, Yin and Zhang, Jun. (2009) "Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation", Volume 228. Elsevier. Journal of Computational Physics 1 137&#8211;146
-fprintf(' %5i — %.2e  &#8211;&#8211; — %.2e  &#8211;&#8211; \n',...
-Mvec(1),NormE1(1),NormE2(1))
-for i=2:M
-fprintf(' %5i — %.2e  %.2f — %.2e  %.2f \n', ...
-Mvec(i),NormE1(i),OrderOr1(i),NormE2(i),OrderOr2(i))
-end
-figure
+<div id="cite-11"></div>
-hold on
+'''[[#citeF-11|[11]]]''' Zhai, Shuying and Feng, Xinlong and He, Yinnian. (2014) "A new method to deduce high-order compact difference schemes for two-dimensional Poisson equation", Volume 230. Elsevier. Applied Mathematics and Computation 9&#8211;26
-plot(x,U,'o')
-plot([x(1:I);alf],[Uexact(1:I)';fun_uL(alf)],'-k','LineWidth',2)
-plot([alf;x(I+1:end)],[fun_uR(alf);Uexact(I+1:end)'],'-k','LineWidth',2)
-plot([alf,alf],[min(U),max(U)],'&#8211;r')
-legend('Numerical','Analytical')
-xlabel('x')
-ylabel('u')
-title(sprintf('Solution (N=%d)',n))
-figure
+<div id="cite-12"></div>
-hold on
+'''[[#citeF-12|[12]]]''' Uh Zapata, Miguel and Itzá Balam, Reymundo. (2017) "High-order implicit finite difference schemes for the two-dimensional Poisson equation", Volume 309. Elsevier. Applied Mathematics and Computation 222&#8211;244
-plot(x,Err,'-k','LineWidth',2)
-plot([alf,alf],[min(Err),max(Err)],'&#8211;r')
-xlabel('x')
-ylabel('—U-u—')
-title(sprintf('Absolute error (N=%d)',n))
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+<div id="cite-13"></div>
-</pre>
+'''[[#citeF-13|[13]]]''' Itzá Balam, Reymundo and Uh Zapata, Miguel. (2020) "A new eighth-order implicit finite difference method to solve the three-dimensional Helmholtz equation", Volume 80. Elsevier. Computers & Mathematics with Applications 5 1176&#8211;1200
+<div id="cite-14"></div>
+'''[[#citeF-14|[14]]]''' Itza Balam, Reymundo and Uh Zapata, Miguel. (2022) "A fourth-order compact implicit immersed interface method for 2D Poisson interface problems", Volume 119. Elsevier. Computers & Mathematics with Applications 257&#8211;277
-|}
+<div id="cite-15"></div>
+'''[[#citeF-15|[15]]]''' Zapata, M Uh and Balam, R Itza and Montalvo-Urquizo, J. (2023) "A compact sixth-order implicit immersed interface method to solve 2D Poisson equations with discontinuities", Volume 210. Elsevier. Mathematics and Computers in Simulation 384&#8211;407
-===BIBLIOGRAPHY===
+<div id="cite-16"></div>
+'''[[#citeF-16|[16]]]''' Peskin, Charles S. (2002) "The immersed boundary method", Volume 11. Cambridge University Press. Acta numerica 479&#8211;517
-<div id="cite-1"></div>
+<div id="cite-17"></div>
-'''[[#citeF-1|[1]]]'''    H. J. Diersch, Fletcher, C.A.J. (1988). Computational Techniques for Fluid Dynamics. Vol. I: Fundamental and General Techniques. Vol. II: Specific Techniques for Different Flow Categories. Springer-Verlag. <div id="cite-2"></div>
+'''[[#citeF-17|[17]]]''' Liu, Xu-Dong and Sideris, Thomas. (2003) "Convergence of the ghost fluid method for elliptic equations with interfaces", Volume 72. Mathematics of computation 244 1731&#8211;1746
-'''[[#citeF-2|[2]]]'''   Sethian, J. A. (1999). Level Set Methods and Fast Marching Methods: Evolving Interfaces in Geometry, Fluid Mechanics, Computer Vision, and Materials Sciences, Cambridge University Press. <div id="cite-3"></div>
-'''[[#citeF-3|[3]]]'''   Li, Z. & Ito, K. (2006). The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains, SIAM: Frontiers in Applied Mathematics. <div id="cite-4"></div>
+<div id="cite-18"></div>
-'''[[#citeF-4|[4]]]'''  Javierre E., Vuik C., Vermolen, F. J. & Van der Zwaag, S. (2006). A comparison of numerical models for one-dimensional Stefan problems. Journal of Computational and Applied Mathematics, 192(2), 445&#8211;459. <div id="cite-5"></div>
+'''[[#citeF-18|[18]]]''' Pan, Kejia and Tan, Yongji and Hu, Hongling. (2010) "An interpolation matched interface and boundary method for elliptic interface problems", Volume 234. Elsevier. Journal of computational and applied mathematics 1 73&#8211;94
-'''[[#citeF-5|[5]]]'''    Shi, Y. E., Ray, R. K., & Nguyen, K. D. (2013). A projection method-based model with the exact C-property for shallow-water flows over dry and irregular bottom using unstructured finite-volume technique. Computers & Fluids, 76, 178&#8211;195. <div id="cite-6"></div>
-'''[[#citeF-6|[6]]]'''   Uh, M. & Xu, S. (2014). The immersed interface method for simulating two-fluid flows, Numerical Mathematics: Theory, Methods and Applications, 7(4), 447&#8211;472. <div id="cite-7"></div>
+<div id="cite-19"></div>
-'''[[#citeF-7|[7]]]'''    Li M., Fornberg B. & Tang T. (1995). A compact fourth order finite difference scheme for the steady incompressible Navier-Stokes equations, Int. J. Numer. Methods Fluids 20, 1137&#8211;1151. <div id="cite-8"></div>
+'''[[#citeF-19|[19]]]''' Mu, Lin and Wang, Junping and Ye, Xiu and Zhao, Shan. (2016) "A new weak Galerkin finite element method for elliptic interface problems", Volume 325. Elsevier. Journal of Computational Physics 157&#8211;173
-'''[[#citeF-8|[8]]]'''      Zhang, J. (2002). Multigrid method and fourth-order compact scheme for 2D Poisson equation with unequal mesh-size discretization. J. Comput. Phys. 179(1), 170&#8211;179. <div id="cite-9"></div>
-'''[[#citeF-9|[9]]]'''    Nabavi, M., Siddiqui, M.H.K., & Dargahi J. (2007). A new 9-point sixth-order accurate compact finite difference method for the Helmholtz equation, J. Sound Vib. 307, 972&#8211;982. <div id="cite-10"></div>
+<div id="cite-20"></div>
-'''[[#citeF-10|[10]]]'''    Wang Y., &  Zhang J. (2009). Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D poisson equation, J. Comput. Phys. 228, 137&#8211;146. <div id="cite-11"></div>
+'''[[#citeF-20|[20]]]''' Cho, Hyuntae and Han, Heejae and Lee, Byungjoon and Ha, Youngsoo and Kang, Myungjoo. (2019) "A second-order boundary condition capturing method for solving the elliptic interface problems on irregular domains", Volume 81. Springer. Journal of Scientific Computing 217&#8211;251
-'''[[#citeF-11|[11]]]'''      Zhai, S., Feng, X., & He, Y. (2014). A new method to deduce high-order compact difference schemes for two-dimensional Poisson equation. Applied Mathematics and Computation, 230, 9&#8211;26. <div id="cite-12"></div>
-'''[[#citeF-12|[12]]]'''  Uh Zapata, M. , & Itzá Balam, R. (2017). High-order implicit finite difference schemes for the two-dimensional Poisson equation. Applied Mathematics and Computation, 309, 222&#8211;244. <div id="cite-13"></div>
+<div id="cite-21"></div>
-'''[[#citeF-13|[13]]]'''  Itzá Balam, R., & Uh Zapata, M (2020). A new eighth-order implicit finite difference method to solve the three-dimensional Helmholtz equation. Computers & Mathematics with Applications, 80(5), 1176&#8211;1200. <div id="cite-14"></div>
+'''[[#citeF-21|[21]]]''' Itzá Balam, Reymundo and Hernandez-Lopez, Francisco and Trejo-Sánchez, Joel and Uh Zapata, Miguel. (2021) "An immersed boundary neural network for solving elliptic equations with singular forces on arbitrary domains", Volume 18. Math. Biosci. Eng 1 22&#8211;56
-'''[[#citeF-14|[14]]]'''  Itza Balam, R., & UhZapata, M. (2022). A fourth-order compact implicit immersed interface method for 2D Poisson interface problems. Computers & Mathematics with Applications, 119, 257&#8211;277. <div id="cite-15"></div>
-'''[[#citeF-15|[15]]]'''  Uh Zapata, M., Itza Balam, R., & Montalvo-Urquizo, J. (2023). A compact sixth-order implicit immersed interface method to solve 2D Poisson equations with discontinuities. Mathematics and Computers in Simulation, 210, 384&#8211;407. <div id="cite-16"></div>
+<div id="cite-22"></div>
-'''[[#citeF-16|[16]]]'''   Peskin, C. S. (2002). The immersed boundary method, Acta Numer. 11, 479&#8211;517. <div id="cite-17"></div>
+'''[[#citeF-22|[22]]]''' LeVeque, Randall J and Li, Zhilin. (1994) "The immersed interface method for elliptic equations with discontinuous coefficients and singular sources", Volume 31. SIAM. SIAM Journal on Numerical Analysis 4 1019&#8211;1044
-'''[[#citeF-17|[17]]]'''  Liu, X. D., & Soderis, T. C. (2003). Convergence of the ghost fluid method for elliptic equations with interfaces, J. Math. Comp. 72, 1731&#8211;1746. <div id="cite-18"></div>
-'''[[#citeF-18|[18]]]'''  Hu, H., Pan, K., & Tan, Y. (2010). An interpolation matched interface and boundary method for elliptic interface problems, J. Comput. Appl. Math. 234, 73&#8211;94. <div id="cite-19"></div>
+<div id="cite-23"></div>
-'''[[#citeF-19|[19]]]'''   Mu, L., Wang, J., Ye, X., & Zhao, S. (2016). A new weak Galerkin finite element method for elliptic interface problems, J. Comput. Phys., 325, 157&#8211;173. <div id="cite-20"></div>
+'''[[#citeF-23|[23]]]''' Wiegmann, Andreas and Bube, Kenneth P. (2000) "The explicit-jump immersed interface method: finite difference methods for PDEs with piecewise smooth solutions", Volume 37. SIAM. SIAM Journal on Numerical Analysis 3 827&#8211;862
-'''[[#citeF-20|[20]]]'''  Cho, H., Han, H., Lee, B., Ha, Y.,  & Kang, M. (2019). A second-order boundary condition capturing method for solving the elliptic interface problems on irregular domains, Journal of Scientific Computing, 81(3), 217&#8211;251. <div id="cite-21"></div>
-'''[[#citeF-21|[21]]]'''  Itzá Balam, R. , Hernandez-Lopez, F., Trejo-Sánchez, J., & Uh Zapata, M (2020). An immersed boundary neural network for solving elliptic equations with singular forces on arbitrary domains. Mathematical Biosciences and Engineering: MBE, 18(1), 22&#8211;56. <div id="cite-22"></div>
+<div id="cite-24"></div>
-'''[[#citeF-22|[22]]]'''   Leveque R.J., & Li Z.  (1994). The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (4), 1019-1044. <div id="cite-23"></div>
+'''[[#citeF-24|[24]]]''' Berthelsen, Petter Andreas. (2004) "A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions", Volume 197. Elsevier. Journal of Computational Physics 1 364&#8211;386
-'''[[#citeF-23|[23]]]'''   Wiegmann A., &  Bube K.P. (2000) The explicit-jump immersed interface method: finite difference methods for PDEs with piecewise smooth solutions, SIAM J. Numer. Anal. 37 (3), 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>
+<div id="cite-25"></div>
-'''[[#citeF-25|[25]]]'''   Seo, J. H. & Mittal, R. (2011). A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries, J. Comput. Phys., 230, 1000&#8211;1019. <div id="cite-26"></div>
+'''[[#citeF-25|[25]]]''' Seo, Jung Hee and Mittal, Rajat. (2011) "A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries", Volume 230. Elsevier. Journal of computational physics 4 1000&#8211;1019
-'''[[#citeF-26|[26]]]'''    Ito, K., Li, Z., & Kyei, Y. (2005). Higher-order, Cartesian grid based finite difference schemes for elliptic equations on irregular domains. SIAM Journal on Scientific Computing, 27(1), 346-367. <div id="cite-27"></div>
-'''[[#citeF-27|[27]]]'''   Gibou, F., & Fedkiw, R. (2005). 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), 577-601. <div id="cite-28"></div>
+<div id="cite-26"></div>
-'''[[#citeF-28|[28]]]'''   Linnick, M. N., & Fasel, H. F. (2005). A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains. J. Comput. Phys., 204(1), 157-192. <div id="cite-29"></div>
+'''[[#citeF-26|[26]]]''' Ito, Kazufumi and Li, Zhilin and Kyei, Yaw. (2005) "Higher-order, Cartesian grid based finite difference schemes for elliptic equations on irregular domains", Volume 27. SIAM. SIAM Journal on Scientific Computing 1 346&#8211;367
-'''[[#citeF-29|[29]]]'''    Zhou, Y. C., Zhao, S., Feig, M., & Wei, G. W. (2006). High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources. J. Comput. Phys., 213(1), 1&#8211;30.  <div id="cite-30"></div>
-'''[[#citeF-30|[30]]]'''   Zhong, X. (2007). A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity. J. Comput. Phys., 225(1), 1066-1099. <div id="cite-31"></div>
+<div id="cite-27"></div>
-'''[[#citeF-31|[31]]]'''   Feng, X., Li, Z., & Qiao, Z. (2011). High order compact finite difference schemes for the Helmholtz equation with discontinuous coefficients. J. of Comput. Math., 324&#8211;340. <div id="cite-32"></div>
+'''[[#citeF-27|[27]]]''' Gibou, Frédéric and Fedkiw, Ronald. (2005) "A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem", Volume 202. Elsevier. Journal of Computational Physics 2 577&#8211;601
-'''[[#citeF-32|[32]]]'''    Pan, K., He, D., & Li, Z. (2021). A high order compact FD framework for elliptic BVPs involving singular sources, interfaces, and irregular domains. J. Sci. Comput. 88(3), 1&#8211;25. <div id="cite-33"></div>
-'''[[#citeF-33|[33]]]'''    Colnago, M., Casaca, W., & de Souza, L. F. (2020). A high-order immersed interface method free of derivative jump conditions for Poisson equations on irregular domains. J. Comput. Phys. 423, 109791. <div id="cite-34"></div>
+<div id="cite-28"></div>
-'''[[#citeF-34|[34]]]'''           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-35"></div>
+'''[[#citeF-28|[28]]]''' Linnick, Mark N and Fasel, Hermann F. (2005) "A high-order immersed interface method for simulating unsteady incompressible flows on irregular domains", Volume 204. Elsevier. Journal of Computational Physics 1 157&#8211;192
-'''[[#citeF-35|[35]]]'''    Liu, Y., & Sen, M. K. (2009). A practical implicit finite-difference method: examples from seismic modeling. Journal of Geophysics and Engineering, 6(3),  231. <div id="cite-36"></div>
-'''[[#citeF-36|[36]]]'''   Xu S., & Wang Z.J.,  (2006). Systematic Derivation of Jump Conditions for the Immersed Interface Method in Three-Dimensional Flow Simulation, J. Sci. Comput., 27(6), 1948-1980. <div id="cite-37"></div>
+<div id="cite-29"></div>
-'''[[#citeF-37|[37]]]'''  Feng, X., & Li, Z. (2012). Simplified immersed interface methods for elliptic interface problems with straight interfaces. Num. Meth. for Par. Diff. Eqs. 28(1), 188&#8211;203.
+'''[[#citeF-29|[29]]]''' Zhou, YC and Zhao, Shan and Feig, Michael and Wei, Guo-Wei. (2006) "High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources", Volume 213. Elsevier. Journal of Computational Physics 1 1&#8211;30
+<div id="cite-30"></div>
+'''[[#citeF-30|[30]]]''' Zhong, Xiaolin. (2007) "A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity", Volume 225. Elsevier. Journal of Computational Physics 1 1066&#8211;1099
+<div id="cite-31"></div>
+'''[[#citeF-31|[31]]]''' Feng, Xiufang and Li, Zhilin and Qiao, Zhonghua. (2011) "High order compact finite difference schemes for the Helmholtz equation with discontinuous coefficients". JSTOR. Journal of Computational Mathematics 324&#8211;340
+<div id="cite-32"></div>
+'''[[#citeF-32|[32]]]''' Pan, Kejia and He, Dongdong and Li, Zhilin. (2021) "A high order compact FD framework for elliptic BVPs involving singular sources, interfaces, and irregular domains", Volume 88. Springer. Journal of Scientific Computing 3 67
+<div id="cite-33"></div>
+'''[[#citeF-33|[33]]]''' Colnago, Marilaine and Casaca, Wallace and de Souza, Leandro Franco. (2020) "A high-order immersed interface method free of derivative jump conditions for Poisson equations on irregular domains", Volume 423. Elsevier. Journal of Computational Physics 109791
+<div id="cite-34"></div>
+'''[[#citeF-34|[34]]]''' Claerbout, Jon F. (1985) "The craft of wave-field extrapolation, in Imaging the earth's interior", Volume 1. Blackwell scientific publications Oxford 260&#8211;265
+<div id="cite-35"></div>
+'''[[#citeF-35|[35]]]''' Liu, Yang and Sen, Mrinal K. (2009) "A practical implicit finite-difference method: examples from seismic modelling", Volume 6. Oxford University Press. Journal of Geophysics and Engineering 3 231&#8211;249
+<div id="cite-36"></div>
+'''[[#citeF-36|[36]]]''' Xu, Sheng and Wang, Z Jane. (2006) "Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation", Volume 27. SIAM. SIAM Journal on Scientific Computing 6 1948&#8211;1980
+<div id="cite-37"></div>
+'''[[#citeF-37|[37]]]''' Feng, Xiufang and Li, Zhilin. (2012) "Simplified immersed interface methods for elliptic interface problems with straight interfaces", Volume 28. Wiley Online Library. Numerical Methods for Partial Differential Equations 1 188&#8211;203