The high accuracy conserved splitting domain decomposition scheme for solving the parabolic equations

The one dimensional parabolic equations with variable coefficients are considered as follows,

∂<!-- ? -->u∂<!-- ? -->t−<!-- - -->∂<!-- ? -->∂<!-- ? -->x(D(x)∂<!-- ? -->u∂<!-- ? -->x)=f(x,t),x∈<!-- ? -->[0,L],t∈<!-- ? -->[0,T],∂<!-- ? -->u∂<!-- ? -->x∣<!-- | -->x=0=0,∂<!-- ? -->u∂<!-- ? -->x|x=L=0,t∈<!-- ? -->[0,T],u(x,0)=u0(x),x∈<!-- ? -->[0,L].$$\begin{array}{} \displaystyle \left\{ \begin{array}{}\frac{\partial u}{\partial t}-\frac{\partial}{\partial x}( D(x)\frac{\partial u}{\partial x})=f(x,t), \qquad x\in [0, L], \quad t\in[0,\,T], \\ \frac{\partial u}{\partial x}\mid_{x=0}=0, \,\frac{\partial u}{\partial x}|_{x=L}=0, \qquad t\in[0,\,T],\\ u(x,0)=u_0(x),\qquad x\in[0, L]. \end{array} \right. \end{array}$$(1)

where D(x) is diffusion coefficients and 0 < D₀ ≤ D(x) ≤ D₁. Let h > 0 be the space step length, and {xi+12}$\begin{array}{} \big\{x_{i+\frac{1}{2}}\big\} \end{array}$ , {x_i} be the uniform staggered partition points as

x12=0,xI+12=L,xi+12=ih,i=1,2,⋯<!-- ? -->,I−<!-- - -->1,xi=(i−<!-- - -->12)h,i=1,2,⋯<!-- ? -->,I.$$\begin{array}{} \displaystyle x_{\frac{1}{2}}=0,\, x_{I+\frac{1}{2}}=L,\, x_{i+\frac{1}{2}}=ih,\, i=1, 2,\cdots, I-1,\, x_{i}=(i-\frac{1}{2})h,\quad i=1, 2,\cdots, I. \end{array}$$(2)

Let τ > 0 be the time step length and tⁿ be the uniform partition points of [0, T] as

t0=0,tM=T,tn=nτ<!-- t -->,n=1,2,⋯<!-- ? -->,M−<!-- - -->1.$$\begin{array}{} \displaystyle t^0=0,\quad t^M = T,\quad t^n=n\tau,\quad n=1, 2,\cdots, M-1. \end{array}$$(3)

For functions F(x, t), define the difference operators as follows,

∂<!-- ? -->τ<!-- t -->Fin=Fin−<!-- - -->Fin−<!-- - -->1τ<!-- t -->,Fin+12=Fin+Fin+12,δ<!-- d -->xFi+12n=Fi+1n−<!-- - -->Finh,Δ<!-- ? -->hFin=δ<!-- d -->xFi+12n−<!-- - -->δ<!-- d -->xFi−<!-- - -->12nh=Fi+1n−<!-- - -->2Fin+Fi−<!-- - -->1nh2,δ<!-- d -->x3Fi+12n=Δ<!-- ? -->h[δ<!-- d -->xFi+12n]=Fi+2n−<!-- - -->3Fi+1n+3Fin−<!-- - -->Fi−<!-- - -->1nh3.$$\begin{array}{} \displaystyle \qquad\qquad\qquad\qquad\partial_{\tau} F_{i}^n=\frac{F_{i}^{n}-F_{i}^{n-1}}{\tau},\quad F_{i}^{n+\frac{1}{2}}=\frac{F_{i}^n+F_{i}^{n+1}}{2},\quad\delta_{x} F_{i+\frac{1}{2}}^n=\frac{F_{i+1}^n-F_{i}^n}{h}, \\ \displaystyle\Delta_h F_i^n=\frac{\delta_{x} F_{i+\frac{1}{2}}^n-\delta_{x} F_{i-\frac{1}{2}}^n}{h}=\frac{F_{i+1}^n-2F_{i}^n+F_{i-1}^n}{h^2},\quad \delta^3_{x} F_{i+\frac{1}{2}}^n=\Delta_h [\delta_{x}F_{i+\frac{1}{2}}^n]=\frac{F_{i+2}^n-3F_{i+1}^n+3F_{i}^n-F_{i-1}^n}{h^3}. \end{array}$$(4)

For the stake of simplicity, the domain [0, L] is divided into two sub-domains and xk+12$\begin{array}{} x_{k+\frac{1}{2}} \end{array}$ is the interface point of the sub-domains. Let {Uin}and{Qi+12n}$\begin{array}{} \{U_i^n\}~~ \text{and}~~\{Q_{i+\frac{1}{2}}^n\} \end{array}$ be the numerical approximations of the exact solutions {uin}$\begin{array}{} \{u_i^n\} \end{array}$ and {Di+12∂<!-- ? -->un∂<!-- ? -->x∣<!-- | -->i+12}$\begin{array}{} \{{D_{i+\frac{1}{2}}\frac{\partial u^n}{\partial x}}\mid_{i+\frac{1}{2}}\} \end{array}$.

Define π<!-- p -->i+12n=δ<!-- d -->xUi+12n,qi+12n=Di+12π<!-- p -->i+12n$\begin{array}{} \pi^n_{i+\frac{1}{2}}=\delta_{x} U^{n}_{i+\frac{1}{2}}, q^n_{i+\frac{1}{2}}=D_{i+\frac{1}{2}}\pi^n_{i+\frac{1}{2}} \end{array}$ and

π<!-- p -->~<!-- ~ -->i+12n=1mH[54(∑<!-- ? -->l=i+1i+mUln−<!-- - -->∑<!-- ? -->l=i−<!-- - -->m+1iUln)−<!-- - -->112(∑<!-- ? -->l=i+m+1i+2mUln−<!-- - -->∑<!-- ? -->l=i−<!-- - -->2m+1i−<!-- - -->mUln)]+h212δ<!-- d -->x3Ui+12n,$$\begin{array}{} \displaystyle \tilde \pi^{n}_{i+\frac{1}{2}}=\frac{1}{mH}[\frac{5}{4}(\sum_{l=i+1}^{i+m}U_l^{n} -\sum_{l=i-m+1}^{i}U_l^{n})-\frac{1}{12}(\sum_{l=i+m+1}^{i+2m}U_l^{n} -\sum_{l=i-2m+1}^{i-m}U_l^{n})]+\frac{h^2}{12}\delta^3_{x} U_{i+\frac{1}{2}}^n, \end{array}$$(5)

and q~<!-- ~ -->i+12n=Di+12π<!-- p -->~<!-- ~ -->ni+12$\begin{array}{} \tilde q^{n}_{i+\frac{1}{2}}={D_{i+\frac{1}{2}}}{\tilde \pi^{n}}_{i+\frac{1}{2}} \end{array}$. Where the large space step H = mh and m is the integer.

Now, we propose the time second-order and space fourth-order conservative domain decomposition scheme of Eqns. (1) in two steps at every time [tⁿ⁻¹, tⁿ].

Step 1

The interface fluxes {Qk+12n+1}$\begin{array}{} \{{Q}_{k+\frac{1}{2}}^{{n+1}}\} \end{array}$ on the interface are firstly computed by

Qk+12n+1=qˇ<!-- ? -->k+12n+1−<!-- - -->h224Δ<!-- ? -->hqˇ<!-- ? -->k+12n+1−<!-- - -->h224Di+12Δ<!-- ? -->hπ<!-- p -->ˇ<!-- ? -->k+12n+1,$$\begin{array}{} \displaystyle Q^{n+1}_{k+\frac{1}{2}}=\check q^{n+1}_{k+\frac{1}{2}}-\frac{{h}^2}{24}\Delta_{h} \check q^{n+1}_{k+\frac{1}{2}}-\frac{{h}^2}{24} D_{i+\frac{1}{2}}\Delta_{h}\check\pi^{n+1}_{k+\frac{1}{2}}, \end{array}$$(6)

where qˇ<!-- ? -->k+12n+1andπ<!-- p -->ˇ<!-- ? -->k+12n+1$\begin{array}{} \check q^{n+1}_{k+\frac{1}{2}}~~ \text{and}~~ \check\pi^{n+1}_{k+\frac{1}{2}} \end{array}$ are computed by the time extrapolation and local multi-point weighted average scheme as follows,

qˇ<!-- ? -->k+12n+1=2q~<!-- ~ -->k+12n−<!-- - -->q~<!-- ~ -->k+12n−<!-- - -->1,π<!-- p -->ˇ<!-- ? -->k+12n+1=2π<!-- p -->~<!-- ~ -->k+12n−<!-- - -->π<!-- p -->~<!-- ~ -->k+12n−<!-- - -->1.$$\begin{array}{} \displaystyle \check q^{n+1}_{k+\frac{1}{2}}=2\tilde q^{n}_{k+\frac{1}{2}}-\tilde q^{{n-1}}_{k+\frac{1}{2}},\quad \check \pi^{n+1}_{k+\frac{1}{2}}=2\tilde \pi^{n}_{k+\frac{1}{2}}-\tilde \pi^{n-1}_{k+\frac{1}{2}}. \end{array}$$(7)

The two-dimensional parabolic equations are considered as

∂<!-- ? -->u∂<!-- ? -->t−<!-- - -->∂<!-- ? -->∂<!-- ? -->x(a1(x,y)∂<!-- ? -->u∂<!-- ? -->x)−<!-- - -->∂<!-- ? -->∂<!-- ? -->y(a2(x,y)∂<!-- ? -->u∂<!-- ? -->y)=f(x,y,t),(x,y,t)∈<!-- ? -->Ω<!-- O -->×<!-- ? -->(0,T],∇<!-- ? -->u⋅<!-- ? -->n→<!-- ? -->=0,(x,y,t)∈<!-- ? -->∂<!-- ? -->Ω<!-- O -->×<!-- ? -->(0,T],u(x,y,0)=u0(x,y),(x,y)∈<!-- ? -->Ω<!-- O -->,$$\begin{array}{} \left\{ \begin{array}{} \frac{\partial u}{\partial t}-\frac{\partial }{\partial x}(a^1(x,y)\frac{\partial u}{\partial x}) -\frac{ \partial }{\partial y} (a^2(x,y)\frac{\partial u}{\partial y})=f(x,y,t), \quad (x,y,t) \in \Omega\times(0,\,T],\\ \nabla u\cdot \vec{n}=0, \quad (x,y,t) \in \partial\,\Omega\times(0,\,T],\\ u{(x,y,0)}=u_0{(x,y)}, \quad (x,y)\in{\Omega}, \end{array} \right. \end{array}$$(25)

where Ω = [0, 1] × [0, 1], a¹(x, y) and a²(x, y) are the diffusion coefficients. Assume that 0 < a₀ ≤ {a¹(x, y), a²(x, y)} ≤ a₁ are the known smooth functions. Define h_x = 1I$\begin{array}{} \frac{1}{I} \end{array}$ and h_y = 1J$\begin{array}{} \frac{1}{J} \end{array}$ be spatial step size along x-directional and y-directional, respectively. I and J are the positive integers. Introducing the following staggered meshes as

xi+12=ihx,i=0,1,⋯<!-- ? -->,I,xi=(i+12)hx,i=0,1,⋯<!-- ? -->,I−<!-- - -->1,yj+12=jhy,j=0,1,⋯<!-- ? -->,J,yj=(j+12)hy,j=0,1,⋯<!-- ? -->,J−<!-- - -->1.$$\begin{array}{} x_{i+\frac{1}{2}}=i h_x,\, i=0,1,\cdots,I,\quad x_{i}=(i+\frac{1}{2}) h_x, \, i=0,1,\cdots, I-1, \\ y_{j+\frac{1}{2}}=j h_y, \, j=0,1,\cdots, J, \quad y_{j}=(j+\frac{1}{2})h_y, \, j=0,1,\cdots,J-1. \end{array}$$(26)

For simplicity, we assume that a¹(x, y) and a²(x, y) are constants and f ≡ 0 as below, when a¹(x, y) and a²(x, y) are variable coefficients, the schemes are modified as similar as 1-dimension problem.

For simplicity of description, we assume that the domain Ω be divided into 2 × 2 block sub-domains (see Fig. 1). Let {(xi+12,yj)}and{(xi,yj+12)}$\begin{array}{} \{(x_{i+\frac{1}{2}},y_j)\}~~ \text{and}~~ \{(x_i,y_{j+\frac{1}{2}})\} \end{array}$ be the nodes for the fluxes while {(x_i, y_j)} are the nodes used for the solution. The line Γ<!-- G -->1x:x=xi1+12$\begin{array}{} \Gamma^x_1: x=x_{i_1+\frac{1}{2}} \end{array}$ is the interface of Ω_1,q and Ω_2,q, q = 1, 2, where i₁ denotes the mesh point index of interface location of Γ<!-- G -->1x$\begin{array}{} \Gamma ^x_1 \end{array}$ along x-direction. The line Γ<!-- G -->1y:y=yj1+12$\begin{array}{} \Gamma^{y}_1: y=y_{j_1+\frac{1}{2}} \end{array}$ is the interface of sub-domains Ω_p,1 and Ω_p,2, p = 1, 2, and j₁ denotes the mesh point index of interface location of Γ<!-- G -->1y$\begin{array}{} \Gamma ^y_1 \end{array}$ along y-direction.

Fig. 1

Let qi+12,jx=a1δ<!-- d -->xUi+12,j,qi,j+12y=a2δ<!-- d -->yUi,j+12$\begin{array}{} q^{x}_{i+\frac{1}{2},j}=a^1 \delta_{x} U_{i+\frac{1}{2},j},~ q^{y}_{i,j+\frac{1}{2}}=a^2 \delta_{y} U_{i,j+\frac{1}{2}} \end{array}$ , and

q~<!-- ~ -->i+12,jx=a1m1H1[54(∑<!-- ? -->l=i+1i+mUl,j−<!-- - -->∑<!-- ? -->l=i−<!-- - -->m+1iUl,j)−<!-- - -->112(∑<!-- ? -->l=i+m+1i+2mUl,j−<!-- - -->∑<!-- ? -->l=i−<!-- - -->2m+1i−<!-- - -->mUl,j)]+hx212a1δ<!-- d -->x3Ui+12,j,q~<!-- ~ -->i,j+12y=a2m2H2[54(∑<!-- ? -->l=i+1i+mUi,l−<!-- - -->∑<!-- ? -->l=i−<!-- - -->m+1iUi,l)−<!-- - -->112(∑<!-- ? -->l=i+m+1i+2mUi,l−<!-- - -->∑<!-- ? -->l=i−<!-- - -->2m+1i−<!-- - -->mUi,l)]+hy212a2δ<!-- d -->y3Ui,j+12,$$\begin{array}{} \displaystyle \tilde q^{x}_{i+\frac{1}{2},j}=\frac{a^1}{m_1H_1}[\frac{5}{4}(\sum_{l=i+1}^{i+m}U_{l,j} -\sum_{l=i-m+1}^{i}U_{l,j})-\frac{1}{12}(\sum_{l=i+m+1}^{i+2m}U_{l,j} -\sum_{l=i-2m+1}^{i-m}U_{l,j})]+\frac{{h_x}^2}{12}a^1\delta^3_{x} U_{i+\frac{1}{2},j}, \\\displaystyle~~ \tilde q^{y}_{i,j+\frac{1}{2}}=\frac{a^2}{m_2H_2}[\frac{5}{4}(\sum_{l=i+1}^{i+m}U_{i,l} -\sum_{l=i-m+1}^{i}U_{i,l})-\frac{1}{12}(\sum_{l=i+m+1}^{i+2m}U_{i,l} -\sum_{l=i-2m+1}^{i-m}U_{i,l})]+\frac{{h_y}^2}{12}a^2\delta^3_{y} U_{i,j+\frac{1}{2}}, \end{array}$$(27)

where H₁ = m₁h_x, H₂ = m₂h_y.

Now we describe the algorithm of our time second-order and space fourth-order conserved splitting domain decomposition scheme on Ω_1,1 at each time [tⁿ⁻¹, tⁿ] in details as

Step 1

Along x-direction.

The intermediate interface fluxes {Qi1+12,jx,n∗<!-- * -->}$\begin{array}{} \{{Q}_{i_1+\frac{1}{2},j}^{x,{n^{*}}}\} \end{array}$ on the interface are firstly computed by

Qi1+12,jx,n∗<!-- * -->=qˇ<!-- ? -->i1+12,jx,n∗<!-- * -->−<!-- - -->hx212Δ<!-- ? -->hxqˇ<!-- ? -->i1+12,jx,n∗<!-- * -->,$$\begin{array}{} \displaystyle Q^{x,n^{*}}_{i_1+\frac{1}{2},j}=\check q^{x,n^{*}}_{i_1+\frac{1}{2},j}-\frac{{h_x}^2}{12}\Delta_{h_x} \check q^{x,n^{*}}_{i_1+\frac{1}{2},j}, \end{array}$$(28)

where qˇ<!-- ? -->i1+12,jx,n∗<!-- * -->$\begin{array}{} \check q^{x,n^{*}}_{i_1+\frac{1}{2},j} \end{array}$ are computed by the time extrapolation and local multi-point weighted average scheme as follows,

qˇ<!-- ? -->i1+12,jx,n∗<!-- * -->=14(5q~<!-- ~ -->i1+12,jx,1−<!-- - -->q~<!-- ~ -->i1+12,jx,0),n=1,2q~<!-- ~ -->i1+12,jx,n−<!-- - -->q~<!-- ~ -->i1+12,jx,n−<!-- - -->1∗<!-- * -->∗<!-- * -->,n≥<!-- = -->2.$$\begin{array}{} \displaystyle \check q^{x,n^{*}}_{i_1+\frac{1}{2},j}=\left\{ \begin{array}{ll} \frac{1}{4}(5\tilde q^{x,1}_{i_1+\frac{1}{2},j}-\tilde q^{x,0}_{i_1+\frac{1}{2},j}),\quad n=1,\\ 2\tilde q^{x,n}_{i_1+\frac{1}{2},j}-\tilde q^{x,{n-1}^{**}}_{i_1+\frac{1}{2},j},\quad n\geq2. \end{array} \right. \end{array}$$(29)

The boundary conditions are approximated by

Q12,jx=0,Qi,12y=0,{(x12,yj),(xi,y12)}∈<!-- ? -->∂<!-- ? -->Ω<!-- O -->h,$$\begin{array}{} \displaystyle Q_{\frac{1}{2},j}^{x}=0,\, Q_{i,\frac{1}{2}}^{y}=0,\, \{(x_{\frac{1}{2}},y_j), (x_i,y_{\frac{1}{2}}) \}\in\partial\Omega_h, \end{array}$$(35)

The initial values are computed by Ui,j0$\begin{array}{} U_{i,j}^0 \end{array}$ = u₀(x_i, y_j), and the first time level values {Ui,j1}$\begin{array}{} \{U^{1}_{i,j}\} \end{array}$ are need to compute by splitting scheme.