m (Muhzapata moved page Draft Uh Zapata 643844164 to Review 499745022132) |
|||
Line 1: | Line 1: | ||
<!-- metadata commented in wiki content | <!-- metadata commented in wiki content | ||
− | ==General foundations of high-order immersed interface methods to solve interface problems== | + | ==General foundations of high-order immersed interface methods to solve interface problems== |
− | '''H. Escamilla | + | '''H. Escamilla Puc<sup>a</sup>, R. Itza Balam<sup>b,c</sup>, M. Uh Zapata<sup>b,c</sup>''' |
{| | {| | ||
|- | |- | ||
− | |<sup>a</sup> Facultad de Matemáticas; Anillo Periférico Norte, Tablaje Cat. 13615, Colonia Chuburná Hidalgo Inn, Mérida, México | + | |<sup>a</sup> Facultad de Matemáticas, Universidad Autónoma de Yucatán; Anillo Periférico Norte, Tablaje Cat. 13615, Colonia Chuburná Hidalgo Inn, Mérida, México |
|} | |} | ||
Line 36: | Line 36: | ||
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 implicit finite-difference and immersed interface method (IFD-IIM) to achieve high-order approximations for second-order derivatives of both continuous and discontinuous real-valued functions. The IFD scheme offers a highly accurate numerical method <span id='citeF-34'></span><span id='citeF-35'></span>[[#cite-34|[34,35]]], while the IIM handles discontinuities through minimal adjustments made exclusively at grid points where the stencil intersects the interface <span id='citeF-36'></span><span id='citeF-37'></span>[[#cite-36|[36,37]]], yielding additional terms known as jump contributions. | ||
− | We illustrate the implementation of | + | 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 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. | ||
Line 42: | Line 42: | ||
==2 Implicit finite-difference formulation with discontinuities== | ==2 Implicit finite-difference formulation with discontinuities== | ||
− | In this section, we outline the key attributes of the implicit finite difference formulation, demonstrating how the scheme can be adapted for addressing discontinuous problems through the utilization of the immersed interface method. | + | In this section, we outline the key attributes of the implicit finite difference formulation to achieve a global fourth-order method, demonstrating how the scheme can be adapted for addressing discontinuous problems through the utilization of the immersed interface method. |
We approximate the numerical solution on the domain <math display="inline">[\mathfrak{a},\,\mathfrak{b}]</math> that is divided into <math display="inline">N</math> sub-intervals, as follows | 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 | ||
Line 62: | Line 62: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-Fig_Scheme.png|594px|Example of a discontinuous function u with an interface located between the points x<sub>I</sub> and x<sub>I+1</sub>.]] |
|- style="text-align: center; font-size: 75%;" | |- 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>. | | 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 87: | ||
The following two Theorems state the main results to approximate the second-order derivative using high-order schemes for regular and irregular points. | The following two Theorems state the main results to approximate the second-order derivative using high-order schemes for regular and irregular points. | ||
− | <span id='theorem-teo:regular'></span>Theorem 1: '''Regular points''' <span id='citeF-34'></span><span id='citeF-35'></span>[[#cite-34|[34,35]]]. Let us consider a real-valued function <math display="inline">u</math> with an interface <math display="inline">x_\alpha </math> such that <math display="inline">x_I\leq x_\alpha{<}x_{I+1}</math>. Then <math display="inline">u_{xx}</math> can be approximated by | + | <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 implicit finite-difference (IFD) scheme |
<span id="eq-2"></span> | <span id="eq-2"></span> | ||
Line 263: | Line 263: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>c_i = \begin{cases}- \dfrac{1}{h^2}\left | + | | 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) | | style="width: 5px;text-align: right;white-space: nowrap;" | (10) | ||
Line 306: | Line 306: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>\left.\mathfrak{D}^2 u \right|_i = \left.\mathfrak{d}^2 | + | | 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) | | style="width: 5px;text-align: right;white-space: nowrap;" | (12) | ||
Line 341: | Line 341: | ||
==3 Poisson equation== | ==3 Poisson equation== | ||
− | In this section, we | + | 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> | <span id="eq-15"></span> | ||
Line 352: | Line 352: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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> | | style="text-align: center;" | <math> u(\boldsymbol{x}) = g(\boldsymbol{x}), \quad \boldsymbol{x}\in \partial \Omega , </math> | ||
Line 362: | Line 362: | ||
|} | |} | ||
− | 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 | + | 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) | + | 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> | <span id="eq-15"></span> | ||
Line 378: | Line 378: | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (15) | | style="width: 5px;text-align: right;white-space: nowrap;" | (15) | ||
|- | |- | ||
− | | style="text-align: center;" | <math> u(x) = g(x), \quad | + | | 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) | | style="width: 5px;text-align: right;white-space: nowrap;" | (16) | ||
|- | |- | ||
Line 389: | Line 389: | ||
|} | |} | ||
− | 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=== | ===3.1 IFD-IIM for the 1D Poisson problem=== | ||
Line 453: | Line 453: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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} + b \left(2 h_R^3\left[u_{xxx}\right]+\dfrac{h_R^4}{2}\left[u_{xxxx}\right]\right)\right), & i=I,\\[4mm] \dfrac{1}{h^2} \left([u]+h_L\left[u_{x}\right]+\dfrac{h_L^2}{2} + b \left(2h_L^3\left[u_{xxx}\right]+\dfrac{h_L^4}{2}\left[u_{xxxx}\right]\right)\right), & i=I+1,\\[2mm] 0, & \textrm{else where,} \end{cases} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (23) | | style="width: 5px;text-align: right;white-space: nowrap;" | (23) | ||
Line 466: | Line 466: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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) | | style="width: 5px;text-align: right;white-space: nowrap;" | (24) | ||
Line 519: | Line 519: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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 | + | | 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},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 | + | | 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\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] | + | | 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> |
|} | |} | ||
|} | |} | ||
− | + | 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> | <span id="eq-27"></span> | ||
Line 552: | Line 552: | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (29) | | 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) | | style="width: 5px;text-align: right;white-space: nowrap;" | (30) | ||
|- | |- | ||
Line 658: | Line 658: | ||
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. | 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 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 | 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 668: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| style="text-align: left; margin:auto;width: 100%;" | ||
|- | |- | ||
− | | style="text-align: center;" | <math>\left\|e \right\|_\infty = \max _{i | + | | 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) | | style="width: 5px;text-align: right;white-space: nowrap;" | (38) | ||
Line 686: | Line 686: | ||
|} | |} | ||
− | 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=== | ===5.1 Example 1. Poison equation with smooth solution=== | ||
Line 708: | Line 708: | ||
− | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto; border-top: 2px solid; border-bottom: 2px solid; min-width: | + | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;" |
− | |+ style="font-size: 75%;" |<span id='table-1'></span>'''Table. 1''' Convergence analysis of Example 1 testing a Poisson equation with smooth solution. | + | |+ style="font-size: 75%;" |<span id='table-1'></span>'''Table. 1''' Convergence analysis of Example 1 testing a Poisson equation with smooth solution. |
|- | |- | ||
− | | | + | | 0.95 @c cccc cccc |
− | | | + | | style="text-align: right;" | |
− | | | + | | |
− | | | + | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
|- | |- | ||
− | | 10 | + | | (rr)2-5 (rr)6-9 <math display="inline">N</math> |
+ | | <math display="inline">L_\infty </math>-norm | ||
+ | | Order | ||
+ | | <math display="inline">L_2</math>-norm | ||
+ | | Order | ||
+ | | <math display="inline">L_\infty </math>-norm | ||
+ | | Order | ||
+ | | <math display="inline">L_2</math>-norm | ||
+ | | Order | ||
+ | |- | ||
+ | | 10 | ||
| 7.52e-02 | | 7.52e-02 | ||
− | | | + | | –- |
− | | | + | | 2.84e-02 |
− | | | + | | –- |
| 9.30e-03 | | 9.30e-03 | ||
− | | | + | | –- |
− | | | + | | 3.56e-03 |
− | | | + | | –- |
|- | |- | ||
− | | | + | | 20 |
| 1.70e-02 | | 1.70e-02 | ||
| 2.15 | | 2.15 | ||
− | | | + | | 6.52e-03 |
| 2.12 | | 2.12 | ||
| 5.05e-04 | | 5.05e-04 | ||
| 4.20 | | 4.20 | ||
− | | | + | | 1.93e-04 |
| 4.21 | | 4.21 | ||
|- | |- | ||
− | | | + | | 40 |
| 4.15e-03 | | 4.15e-03 | ||
| 2.03 | | 2.03 | ||
− | | | + | | 1.60e-03 |
| 2.03 | | 2.03 | ||
| 3.06e-05 | | 3.06e-05 | ||
| 4.04 | | 4.04 | ||
− | | | + | | 1.17e-05 |
| 4.04 | | 4.04 | ||
|- | |- | ||
− | | | + | | 80 |
| 1.03e-03 | | 1.03e-03 | ||
| 2.01 | | 2.01 | ||
− | | | + | | 3.98e-04 |
| 2.01 | | 2.01 | ||
| 1.90e-06 | | 1.90e-06 | ||
| 4.01 | | 4.01 | ||
− | | | + | | 7.27e-07 |
| 4.01 | | 4.01 | ||
|- | |- | ||
Line 764: | Line 768: | ||
| 2.57e-04 | | 2.57e-04 | ||
| 2.00 | | 2.00 | ||
− | | | + | | 9.95e-05 |
| 2.00 | | 2.00 | ||
| 1.18e-07 | | 1.18e-07 | ||
| 4.00 | | 4.00 | ||
− | | | + | | 4.54e-08 |
− | | 4.00 | + | | 4.00 |
|- | |- | ||
| LSM | | LSM | ||
| | | | ||
− | | | + | | '''2.04''' |
| | | | ||
| '''2.04''' | | '''2.04''' | ||
Line 781: | Line 785: | ||
| '''4.06''' | | '''4.06''' | ||
|- | |- | ||
+ | | | ||
|} | |} | ||
Line 793: | Line 798: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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>u(x) = \begin{cases}\sin (\pi x), & x\leq x_\alpha ,\\ \cos (\pi x), & x > x_\alpha{.} \end{cases} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (41) | | style="width: 5px;text-align: right;white-space: nowrap;" | (41) | ||
Line 803: | Line 808: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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.]] |
|- style="text-align: center; font-size: 75%;" | |- 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>. | ||
Line 810: | Line 815: | ||
− | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width: | + | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;" |
− | |+ style="font-size: 75%;" |<span id='table-2'></span>'''Table. 2''' Convergence analysis of Example 2 using the IFD-IIM. | + | |+ style="font-size: 75%;" |<span id='table-2'></span>'''Table. 2''' Convergence analysis of Example 2 using the IFD-IIM. |
|- | |- | ||
− | | <math display="inline">N</math> | + | | 0.95 @c cccc cccc |
+ | | style="text-align: right;" | | ||
+ | | | ||
+ | |- | ||
+ | | (rrrr)2-5 (rrrr)6-9 | ||
+ | | | ||
+ | | | ||
+ | | | ||
+ | | | ||
+ | |- | ||
+ | | (rr)2-3 (rr)4-5 (rr)6-7 (rr)8-9 <math display="inline">N</math> | ||
| <math display="inline">L_\infty </math>-norm | | <math display="inline">L_\infty </math>-norm | ||
| Order | | Order | ||
Line 823: | Line 838: | ||
| Order | | Order | ||
|- | |- | ||
− | | | + | | 10 |
| 1.69e-02 | | 1.69e-02 | ||
− | | | + | | –- |
| 5.19e-03 | | 5.19e-03 | ||
− | | | + | | –- |
| 5.75e-05 | | 5.75e-05 | ||
− | | | + | | –- |
| 3.42e-05 | | 3.42e-05 | ||
− | | | + | | –- |
|- | |- | ||
− | | | + | | 20 |
| 4.21e-03 | | 4.21e-03 | ||
| 2.00 | | 2.00 | ||
Line 843: | Line 858: | ||
| 4.12 | | 4.12 | ||
|- | |- | ||
− | | | + | | 40 |
| 1.05e-03 | | 1.05e-03 | ||
| 2.00 | | 2.00 | ||
Line 853: | Line 868: | ||
| 3.82 | | 3.82 | ||
|- | |- | ||
− | | | + | | 80 |
| 2.63e-04 | | 2.63e-04 | ||
| 2.00 | | 2.00 | ||
Line 863: | Line 878: | ||
| 4.16 | | 4.16 | ||
|- | |- | ||
− | | | + | | 160 |
| 6.57e-05 | | 6.57e-05 | ||
| 2.00 | | 2.00 | ||
Line 873: | Line 888: | ||
| 3.85 | | 3.85 | ||
|- | |- | ||
− | | | + | | LSM |
| | | | ||
− | | | + | | '''2.04''' |
| | | | ||
| '''2.02''' | | '''2.02''' | ||
Line 883: | Line 898: | ||
| '''3.99''' | | '''3.99''' | ||
|- | |- | ||
+ | | | ||
|} | |} | ||
<div id='img-3'></div> | <div id='img-3'></div> | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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.]] |
|- style="text-align: center; font-size: 75%;" | |- 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 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. [[#img-4|4]] shows that if we add additional jumps of high-order derivatives | + | 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. [[#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. [[#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> | <div id='img-4'></div> | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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ₓₓₓ]).]] |
|- style="text-align: center; font-size: 75%;" | |- 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 term includes jumps up to fifth-order (<math>[u_{xxxxx}]=[f_{xxx}]</math>). | ||
Line 905: | Line 921: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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ₓ]).]] |
|- style="text-align: center; font-size: 75%;" | |- 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 term includes jumps up to third-order (<math>[u_{xxx}]=[f_{x}]</math>). | ||
|} | |} | ||
− | Note that there are some points in Fig. [[#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. [[#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 [[#3.2 IFD-IIM and principal jump conditions|3.2]]. |
===5.3 Example 3. Poison equation with multiple interface points=== | ===5.3 Example 3. Poison equation with multiple interface points=== | ||
Line 921: | Line 937: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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\leq x_{\alpha _1},\\ \cos (2\pi x), & x_{\alpha _1}<x \leq x_{\alpha _2}, \\ \sin (2\pi x), & x_{\alpha _2} < x. \end{cases} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (42) | | style="width: 5px;text-align: right;white-space: nowrap;" | (42) | ||
Line 931: | Line 947: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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.]] |
|- style="text-align: center; font-size: 75%;" | |- 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 different grid resolutions. | ||
Line 938: | Line 954: | ||
==6 Numerical results for the Heat equation== | ==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=== | ===6.1 Example 4. Heat equation with smooth solution=== | ||
Line 958: | Line 974: | ||
− | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width: | + | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;" |
|+ style="font-size: 75%;" |<span id='table-3'></span>'''Table. 3''' Convergence analysis of Example 4 testing a heat equation with smooth solution. | |+ style="font-size: 75%;" |<span id='table-3'></span>'''Table. 3''' Convergence analysis of Example 4 testing a heat equation with smooth solution. | ||
|- | |- | ||
− | | <math display="inline">N</math> | + | | 0.95 @c cccc cccc |
+ | | style="text-align: right;" | | ||
+ | | | ||
+ | |- | ||
+ | | (rr)2-5 (rr)6-9 <math display="inline">N</math> | ||
| <math display="inline">L_\infty </math>-norm | | <math display="inline">L_\infty </math>-norm | ||
| Order | | Order | ||
Line 971: | Line 991: | ||
| Order | | Order | ||
|- | |- | ||
− | | | + | | 10 |
| 1.66e-02 | | 1.66e-02 | ||
− | | | + | | –- |
− | | | + | | 1.07e-02 |
− | | | + | | –- |
| 1.84e-04 | | 1.84e-04 | ||
− | | | + | | –- |
| 1.18e-04 | | 1.18e-04 | ||
− | | | + | | –- |
|- | |- | ||
− | | | + | | 20 |
| 4.12e-03 | | 4.12e-03 | ||
| 2.01 | | 2.01 | ||
− | | | + | | 2.65e-03 |
| 2.01 | | 2.01 | ||
| 1.14e-05 | | 1.14e-05 | ||
Line 991: | Line 1,011: | ||
| 4.01 | | 4.01 | ||
|- | |- | ||
− | | | + | | 40 |
| 1.03e-03 | | 1.03e-03 | ||
| 2.00 | | 2.00 | ||
− | | | + | | 6.60e-04 |
| 2.00 | | 2.00 | ||
| 7.15e-07 | | 7.15e-07 | ||
Line 1,001: | Line 1,021: | ||
| 4.00 | | 4.00 | ||
|- | |- | ||
− | | | + | | 80 |
| 2.58e-04 | | 2.58e-04 | ||
| 2.00 | | 2.00 | ||
− | | | + | | 1.65e-04 |
| 2.00 | | 2.00 | ||
| 4.47e-08 | | 4.47e-08 | ||
Line 1,011: | Line 1,031: | ||
| 4.00 | | 4.00 | ||
|- | |- | ||
− | | | + | | 160 |
| 6.44e-05 | | 6.44e-05 | ||
| 2.00 | | 2.00 | ||
− | | | + | | 4.12e-05 |
| 2.00 | | 2.00 | ||
| 2.79e-09 | | 2.79e-09 | ||
Line 1,021: | Line 1,041: | ||
| 4.00 | | 4.00 | ||
|- | |- | ||
− | | | + | | LSM |
| | | | ||
− | | | + | | '''2.00''' |
| | | | ||
| '''2.00''' | | '''2.00''' | ||
Line 1,031: | Line 1,051: | ||
| '''4.00''' | | '''4.00''' | ||
|- | |- | ||
+ | | | ||
|} | |} | ||
Line 1,043: | Line 1,064: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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>u(x) = \begin{cases}\sin (x)\exp (-t), & x\leq x_\alpha ,\\ \cos (x)\exp (-t), & x > x_\alpha{.} \end{cases} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (44) | | style="width: 5px;text-align: right;white-space: nowrap;" | (44) | ||
Line 1,050: | Line 1,071: | ||
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 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. [[#img-7|7]]. On the other hand, the one-dimensional results at <math display="inline"> | + | The figure of numerical solution and absolute error plot using <math display="inline">N=40</math> for different time stages are shown in Fig. [[#img-7|7]]. On the other hand, the one-dimensional results at <math display="inline">T=0.5</math> are presented in Figs. [[#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. |
<div id='img-7'></div> | <div id='img-7'></div> | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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].]] |
|- style="text-align: center; font-size: 75%;" | |- 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>. | | 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,062: | Line 1,083: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-Fig_Heat_SolEx5_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%;" | |- 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. | | 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,069: | Line 1,090: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-Fig_Heat_SolEx5_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%;" | |- 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. | | 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: | + | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;" |
|+ style="font-size: 75%;" |<span id='table-4'></span>'''Table. 4''' Convergence analysis of Example 5 using the IFD-IIM for the heat equation with a discontinuous solution. | |+ 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. | ||
|- | |- | ||
− | | <math display="inline">N</math> | + | | 0.95 @c cccc cccc |
+ | | style="text-align: right;" | | ||
+ | | | ||
+ | |- | ||
+ | | (rr)2-5 (rr)6-9 <math display="inline">N</math> | ||
| <math display="inline">L_\infty </math>-norm | | <math display="inline">L_\infty </math>-norm | ||
| Order | | Order | ||
Line 1,087: | Line 1,112: | ||
| Order | | Order | ||
|- | |- | ||
− | | | + | | 10 |
| 1.11e-07 | | 1.11e-07 | ||
− | | | + | | –- |
− | | | + | | 6.52e-08 |
− | | | + | | –- |
| 7.57e-08 | | 7.57e-08 | ||
− | | | + | | –- |
| 4.49e-08 | | 4.49e-08 | ||
− | | | + | | –- |
|- | |- | ||
| 20 | | 20 | ||
| 6.90e-09 | | 6.90e-09 | ||
| 4.01 | | 4.01 | ||
− | | | + | | 4.01e-09 |
| 4.02 | | 4.02 | ||
| 3.75e-09 | | 3.75e-09 | ||
Line 1,110: | Line 1,135: | ||
| 4.29e-10 | | 4.29e-10 | ||
| 4.00 | | 4.00 | ||
− | | | + | | 2.49e-10 |
| 4.01 | | 4.01 | ||
| 3.58e-10 | | 3.58e-10 | ||
Line 1,120: | Line 1,145: | ||
| 2.68e-11 | | 2.68e-11 | ||
| 4.00 | | 4.00 | ||
− | | | + | | 1.55e-11 |
| 4.00 | | 4.00 | ||
| 1.53e-11 | | 1.53e-11 | ||
Line 1,127: | Line 1,152: | ||
| 4.55 | | 4.55 | ||
|- | |- | ||
− | | | + | | 160 |
| 1.63e-12 | | 1.63e-12 | ||
| 4.03 | | 4.03 | ||
− | | | + | | 9.39e-13 |
| 4.04 | | 4.04 | ||
| 1.35e-12 | | 1.35e-12 | ||
Line 1,137: | Line 1,162: | ||
| 3.51 | | 3.51 | ||
|- | |- | ||
− | | | + | | LSM |
| | | | ||
| '''4.01''' | | '''4.01''' | ||
Line 1,147: | Line 1,172: | ||
| '''3.96''' | | '''3.96''' | ||
|- | |- | ||
+ | | | ||
|} | |} | ||
From Table [[#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 [[#table-5|5]] can be obtained for smooth solutions. | From Table [[#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 [[#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: | + | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;" |
− | |+ style="font-size: 75%;" |<span id='table-5'></span>'''Table. 5''' Convergence analysis of Example 5 testing the heat equation with <math>\Delta t=h</math>. | + | |+ 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>. |
|- | |- | ||
− | | <math display="inline">N</math> | + | | 0.95 @c cccc cccc |
+ | | | ||
+ | | | ||
+ | |- | ||
+ | | (rr)2-5 (rr)6-9 <math display="inline">N</math> | ||
| <math display="inline">L_\infty </math>-norm | | <math display="inline">L_\infty </math>-norm | ||
| Order | | Order | ||
Line 1,164: | Line 1,194: | ||
| Order | | Order | ||
|- | |- | ||
− | | | + | | 10 |
| 2.94e-04 | | 2.94e-04 | ||
− | | | + | | –- |
− | | | + | | 1.87e-04 |
− | | | + | | –- |
| 3.78e-05 | | 3.78e-05 | ||
− | | | + | | –- |
− | | | + | | 2.67e-05 |
− | | | + | | –- |
|- | |- | ||
− | | | + | | 20 |
| 8.53e-05 | | 8.53e-05 | ||
| 1.78 | | 1.78 | ||
− | | | + | | 4.85e-05 |
| 1.95 | | 1.95 | ||
| 1.03e-05 | | 1.03e-05 | ||
| 1.87 | | 1.87 | ||
− | | | + | | 7.22e-06 |
| 1.89 | | 1.89 | ||
|- | |- | ||
− | | | + | | 40 |
| 2.11e-05 | | 2.11e-05 | ||
| 2.01 | | 2.01 | ||
− | | | + | | 1.26e-05 |
| 1.94 | | 1.94 | ||
| 2.74e-06 | | 2.74e-06 | ||
| 1.91 | | 1.91 | ||
− | | | + | | 1.92e-06 |
| 1.91 | | 1.91 | ||
|- | |- | ||
− | | | + | | 80 |
| 5.31e-06 | | 5.31e-06 | ||
| 1.99 | | 1.99 | ||
− | | | + | | 3.15e-06 |
| 2.00 | | 2.00 | ||
| 6.91e-07 | | 6.91e-07 | ||
| 1.99 | | 1.99 | ||
− | | | + | | 4.83e-07 |
| 1.99 | | 1.99 | ||
|- | |- | ||
− | | | + | | 160 |
| 1.34e-06 | | 1.34e-06 | ||
| 1.99 | | 1.99 | ||
− | | | + | | 7.83e-07 |
| 2.00 | | 2.00 | ||
| 1.72e-07 | | 1.72e-07 | ||
| 2.00 | | 2.00 | ||
− | | | + | | 1.21e-07 |
| 2.00 | | 2.00 | ||
|- | |- | ||
− | | | + | | LSM |
| | | | ||
− | | | + | | '''1.96''' |
| | | | ||
| '''1.97''' | | '''1.97''' | ||
Line 1,224: | Line 1,254: | ||
| '''1.95''' | | '''1.95''' | ||
|- | |- | ||
+ | | | ||
|} | |} | ||
Line 1,236: | Line 1,267: | ||
{| style="text-align: left; margin:auto;width: 100%;" | {| 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)\exp (-t), & x\leq x_{\alpha },\\ \cos (2\pi x)\exp (-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;" | (45) | ||
Line 1,248: | Line 1,279: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-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].]] |
|- style="text-align: center; font-size: 75%;" | |- 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>. | | 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,255: | Line 1,286: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-Fig_Heat_SolEx6_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%;" | |- 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. | | 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,262: | Line 1,293: | ||
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
|- | |- | ||
− | |[[Image: | + | |[[Image:Review_499745022132-Fig_Heat_SolEx6_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%;" | |- 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. | | 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: | + | {| class="floating_tableSCP" style="text-align: left; margin: 1em auto;border-top: 2px solid;border-bottom: 2px solid;min-width:50%;" |
− | |+ style="font-size: 75%;" |<span id='table-6'></span>'''Table. 6''' Convergence analysis of Example 6 using the IFD-IIM for the heat equation with a general discontinuous solution. | + | |+ 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. |
|- | |- | ||
− | | <math display="inline">N</math> | + | | 0.95 @c cccc cccc |
+ | | style="text-align: right;" | | ||
+ | | | ||
+ | |- | ||
+ | | (rr)2-5 (rr)6-9 <math display="inline">N</math> | ||
| <math display="inline">L_\infty </math>-norm | | <math display="inline">L_\infty </math>-norm | ||
| Order | | Order | ||
Line 1,280: | Line 1,315: | ||
| Order | | Order | ||
|- | |- | ||
− | | | + | | 10 |
| 3.91e-04 | | 3.91e-04 | ||
− | | | + | | –- |
− | | | + | | 2.18e-04 |
− | | | + | | –- |
| 4.32e-04 | | 4.32e-04 | ||
− | | | + | | –- |
| 2.43e-04 | | 2.43e-04 | ||
− | | | + | | –- |
|- | |- | ||
− | | | + | | 20 |
| 2.50e-05 | | 2.50e-05 | ||
| 3.97 | | 3.97 | ||
− | | | + | | 1.34e-05 |
| 4.02 | | 4.02 | ||
| 2.48e-05 | | 2.48e-05 | ||
Line 1,300: | Line 1,335: | ||
| 4.07 | | 4.07 | ||
|- | |- | ||
− | | | + | | 40 |
| 1.56e-06 | | 1.56e-06 | ||
| 4.00 | | 4.00 | ||
− | | | + | | 8.37e-07 |
| 4.01 | | 4.01 | ||
| 2.38e-06 | | 2.38e-06 | ||
Line 1,310: | Line 1,345: | ||
| 3.69 | | 3.69 | ||
|- | |- | ||
− | | | + | | 80 |
| 9.72e-08 | | 9.72e-08 | ||
| 4.00 | | 4.00 | ||
− | | | + | | 5.22e-08 |
| 4.00 | | 4.00 | ||
| 9.48e-08 | | 9.48e-08 | ||
Line 1,320: | Line 1,355: | ||
| 4.34 | | 4.34 | ||
|- | |- | ||
− | | | + | | 160 |
| 6.08e-09 | | 6.08e-09 | ||
| 4.03 | | 4.03 | ||
− | | | + | | 3.26e-09 |
| 4.00 | | 4.00 | ||
| 9.29e-09 | | 9.29e-09 | ||
Line 1,330: | Line 1,365: | ||
| 3.69 | | 3.69 | ||
|- | |- | ||
− | | | + | | LSM |
| | | | ||
| '''4.00''' | | '''4.00''' | ||
Line 1,340: | Line 1,375: | ||
| '''3.95''' | | '''3.95''' | ||
|- | |- | ||
+ | | | ||
|} | |} | ||
Line 1,348: | Line 1,384: | ||
==Acknowledgments== | ==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. |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | for | + | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
===BIBLIOGRAPHY=== | ===BIBLIOGRAPHY=== |
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
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 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 [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.
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.
In this section, we outline the key attributes of the implicit finite difference formulation to achieve a global fourth-order method, demonstrating how the scheme can be adapted for addressing discontinuous problems through the utilization of the immersed interface method.
We approximate the numerical solution on the domain that is divided into sub-intervals, as follows
|
(1) |
where the grid size is given by . We employ and to denote the approximate and exact solution at the th-point of the grid, respectively. Here, the interface is at located between grid points and as , see Fig. 1. The distances of the closest grid points to the interface are defined as and . Note that and 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 at is defined as
|
We employ a similar definition for the jumps such as the ones of the right-hand side and the derivatives of .
In this paper, we designate and 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.
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 with an interface such that . Then can be approximated by a fourth-order implicit finite-difference (IFD) scheme
|
(2) |
where
|
(3) |
, and central finite-difference formula is given by
|
(4) |
From Taylor series expansions and under some simplifications, the second-order derivative at any regular point 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 . 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 , , and .
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 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 [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 even if the local truncation error at and is .
Theorem 2: Irregular points [14]. Let us consider the known jump conditions
|
(5) |
at such that . Then can be approximated at and by the implicit finite-difference immersed interface method (IFD-IIM) given by
|
(6) |
where and are defined in (3) and (4), respectively, and
|
(7) |
and where .
We obtain a third-order scheme for at and 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 such that we can apply the standard central scheme to and . For instance, a function based on the original left solution is defined as
|
Using Taylor series expansions of around , the definition of jumps (5), and some simplification yield
|
Thus,
|
Finally, we get (6) which complete the proof. The same procedure can be applied for the proof at . We refer to the reader to the work of Itza Balam and Uh Zapata [14] for more details.
Remark 1: If , 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 as follows
|
(9) |
where
|
(10) |
Note that, in this case, we only require to explicitly know jump conditions , and .
Before to apply the previous results for approximations to differential equations, it would be useful to express with finite differences the operator for a real-valued function and not its second-order derivative . The finite-difference formula is obtained by approximating with the central finite differences, as presented in (1)-(10) for regular and irregular points, respectively.
For regular points (), it follows from equation (1):
|
Thus, if we define as the resulting finite-difference
|
(11) |
then, we have that
|
(12) |
However, we still have the same issue to overcome for a discontinuous function . We remark that this second-order derivative of is already multiplied for . Consequently, we can use the IIM technique where the contribution term is only first-order accurate to ; thus, we still have a local to keep a global fourth-order accurate method. Then for irregular points, the IIM applied for this term follows
|
(13) |
where is given by finite difference (11) and
|
(14) |
Finally, for regular and irregular points, we have high-order finite-difference approximations of the implicit operator applied to a real-valued function as in (12) and (13), and to its second-order derivative as in (2) and (6). Now, we proceed to implement them to the solution of differential equations, as presented in the following two sections.
In this section, we develop a fourth-order finite difference scheme for the Poisson equation. Let us consider and as the solution of the problem and known right-hand side function, respectively. Thus the interface problem is given by
|
We divide in two regions, and , separated by an immersed interface . Dirichlet boundary conditions are defined on through function . We assume that and may have discontinuities at the interface . Thus, we require additional conditions and known as principal jump conditions. Here, is the derivative in the normal direction. Note that the principal jump conditions, and , are known functions defined on .
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
|
Here, and can be discontinuous functions at a given fixed point , and principal jump conditions and are known values at .
Let us consider that is located between the adjacent grid points and as . If we apply operator at both sides of (15), then we get
|
(19) |
For and , 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 , 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 and 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 is given by
|
(22) |
where total contribution is with
|
(23) |
and
|
(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 method is expected.
Note that , , , , and must be known to apply the proposal fourth-order IFD-IIM. Thus, it seems that more jump conditions of 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
|
(26) |
Thus contribution depends only on the principal jump conditions, and , and right-hand side jumps: , , and . The additional jumps derivatives from the right-hand side can be approximated using the current values of . In this paper, we will assume that we know them.
Remark 2: For the 1D Poisson problem, a global second-order IIM () only requires to know the principal jump conditions , and .
Remark 3: If , then and both weight terms next to second-order derivative jump of are equal to zero in (26). Thus, we do not require to know jump condition to obtain a fourth-order method when the interface is located at a grid point.
For the second differential equation, we consider the heat conduction equation with initial, boundary, and principal jump conditions, as follows
|
where the source and initial value may be discontinuous or singular across a fixed interface . The interface is immersed in the domain and divides into two parts, and . As the Poisson equation, , and are known functions. This paper only focuses on the one-dimensional problem given by
|
where the source and initial value may be discontinuous or singular across a fixed interface located at .
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) yields
|
(33) |
For regular points, using the IFD method (2) and approximation (12), equation (33) can be approximated as follows
|
(34) |
For irregular points, using the IFD-IIM (6) and approximation (13), the implicit scheme is given by
|
(35) |
where
|
(36) |
Here, and are defined as (14), and is given by (7). Thus, for 1D heat conduction equation (27), the IFD-IIM reads
|
where and is given by (36). Here, we have that for regular points.
Remark 4: The IFD-IIM (37) is unconditionally stable and the local truncation error is of and for regular and irregular grid points, respectively. Thus a global fourth-order method is expected by taking .
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 and its derivatives, the finite-difference scheme presented in equation (37) requires additional knowledge of the jump conditions for the solution given by , , and . 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].
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 , and the grid spacing is for different sub-intervals. The errors are reported using the -norm and the discrete -norm calculated as
|
(38) |
respectively, where and corresponds to the numerical and exact solution at , respectively. The estimated order of accuracy is computed as
|
(39) |
where and indicates the different number of sub-intervals of the corresponding norm.
Example 1 considers an infinitely smooth and differentiable exact solution of the one-dimensional Poisson problem (15) given by the following expression
|
(40) |
The right-hand side function, , is obtained directly from (40). We impose Dirichlet boundary conditions according to the function . Due to the regularity of the solution, the jump contributions and in equation (26) are equal to zero.
Table 1 presents the convergence analysis of Example 1 for different grid resolutions. If , 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 . Last row of table 1 shows the numerical order calculated by the regression-line slope based on a least squares method (LSM) of the - and -norm error. A complete analysis of the IFD method for smooth solutions can be found in [12].
0.95 @c cccc cccc | ||||||||
(rr)2-5 (rr)6-9 | -norm | Order | -norm | Order | -norm | Order | -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 | ||||
For Example 2, we show the method's capability by solving a single interface problem located at . The exact solution is given by the function
|
(41) |
The right-hand side, , is obtained directly from equation (41). We test two different points: and . For the first case, we always have for (); thus the interface is always located at one grid point of that resolution. In general, for , we have different values for different numbers. Figure 2 shows the numerical and exact solution when the interface is located at these two values using . 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) . |
Table 2 shows the convergence analysis for Example 2 for the two values. Observe that a second-order method is recovered for and a fourth-order method for . 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 3 shows the error analysis corresponding to interface locations and for .
0.95 @c cccc cccc | ||||||||
(rrrr)2-5 (rrrr)6-9 | ||||||||
(rr)2-3 (rr)4-5 (rr)6-7 (rr)8-9 | -norm | Order | -norm | Order | -norm | Order | -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 | ||||
![]() |
Figure 3: Convergence analysis of Example 2 for using (a) , and (b) . |
The contribution formula includes jumps , , , , and to obtain a fourth-order accurate method. Fig. 4 shows that if we add additional jumps of high-order derivatives into , such as , we observe that the error oscillation decreases in comparison with Fig. 3 results. It is expected because the method is for the whole computational domain, including the irregular points. Thus, we can mitigate error oscillations due to interface position by adding high-order jumps.
![]() |
Figure 4: Convergence analysis of Example 2 for using (a) , and (b) . The contribution term includes jumps up to fifth-order (). |
Now we study the effect of removing jumps in . Let us consider a contribution term up to the third derivative jump by dropping . 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 random distribution depending on the interface location. In general, the error magnitude perturbations are related to the variations coming from and values in the formula (26).
![]() |
Figure 5: Convergence analysis of Example 2 for using (a) , and (b) . The contribution term 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 : (a) for , and (b) for . Both set of points satisfy that the interface is located at a grid point . In the case of , we observe that the global order (black points) is close to ; meanwhile, the circled points order is four. Similar behavior is obtained for . Thus, for a given mesh resolution , the IFD-IIM is still fourth-order accurate for , as proven theoretically in Section 3.2.
Example 3 investigates the numerical scheme capability to solve a multiple interface problem. We only focus on two interface points located at and . 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
|
(42) |
The right-hand function is obtained directly from the exact solution (42). 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 and using . 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 different grid resolutions. |
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 and space step, , used for the Poisson examples. Here, the final time is at and the time step is . The error and estimated order of accuracy are reported using (38) and (39), respectively, at the final step approximation.
This example is constructed so that the exact solution is
|
(43) |
where the source term is derived from (27) and (43). The initial and boundary conditions are also obtained from the exact solution. The convergence analysis of this example is presented in Table 3 with the final time being . As expected, the high-order methods reach their corresponding order of accuracy.
0.95 @c cccc cccc | ||||||||
(rr)2-5 (rr)6-9 | -norm | Order | -norm | Order | -norm | Order | -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 | ||||
In this example, the exact solution is defined as
|
(44) |
The source term of this problem is . In this example, both the function and their derivatives have jumps, and these jumps vary in time. The boundary condition, initial condition, and all jumps are specified by .
The figure of numerical solution and absolute error plot using for different time stages are shown in Fig. 7. On the other hand, the one-dimensional results at are presented in Figs. 8 and 9. More details of the grid refinement analysis at 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 is because the different 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. |
0.95 @c cccc cccc | ||||||||
(rr)2-5 (rr)6-9 | -norm | Order | -norm | Order | -norm | Order | -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 . 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.
0.95 @c cccc cccc | ||||||||
(rr)2-5 (rr)6-9 | -norm | Order | -norm | Order | -norm | Order | -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 | ||||
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
|
(45) |
Here, the source term is a general function that varies both in time and space. It is directly derived from equation (27) as well as the exact solution (45). 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 and . It is noteworthy that the behavior of the solution differs from the previous example. More in-depth analysis of the results at are illustrated in Figs. 11 and 12, corresponding to values of and , 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 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. |
0.95 @c cccc cccc | ||||||||
(rr)2-5 (rr)6-9 | -norm | Order | -norm | Order | -norm | Order | -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 | ||||
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 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.
[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–203.
Published on 16/11/23
Submitted on 19/09/23
Licence: CC BY-NC-SA license
Are you one of the authors of this document?