Boundary corrections on multi-dimensional PDEs

In this paper, we firstly consider multi-dimensional semi-linear parabolic partial differential equations on rectangular domains \(\Omega =(0,1)^N\subset \mathbb {R}^N\) and \(t\in [0,t^*]\),

with time-dependent Dirichlet boundary conditions

and an initial condition

Here and henceforth, we denote \(\partial _{t} u\equiv \frac{\partial u}{\partial t},\; \partial ^2_{x_i x_i} u\equiv \frac{\partial ^2 u}{\partial x_i^2},\) and \(\partial \Omega \) is the boundary of \(\Omega \) (\(\overline{\Omega }= \Omega \cup \partial \Omega \)). The term \(c(t,\textbf{x},u)\) stands for the reaction term, and we will assume that it is either a non-stiff term (\(\partial c/\partial u\) has a Lipschitz constant of moderate size) or it is simply a source term (it does not depend on u). The diffusion coefficients are supposed to be non-negative, i.e., \(d_i(\textbf{x})\ge 0, \forall \textbf{x}\in \overline{\Omega }\). Although we consider that the spatial domain is the N-dimensional unit square, the results extend straightforwardly to any N-dimensional box.

Our focus is on the time integration of the differential systems that come from some spatial discretization of (1), (2), (3) by using finite differences. In particular, we are interested in AMF-type W-methods (based on approximate matrix factorization) which avoid the solution of non-linear equations and only require the solution of linear systems with tridiagonal matrices, as it can be seen in [5,6,7]. They have also been successful in applications to the case of variable coefficients in [8] and multi-dimensional problems (\(N > 3\)) in [13].

These methods belong to the class of ADI methods, and they, as well as other splitting methods, suffer from order reduction when time-dependent boundary conditions are considered. Many authors have proposed correction techniques to restore the order [4, 12]. In [10, IV.2.4], it is stated that “the general idea is always to treat boundary values as far as possible in the same way as the interior points.”

Following this idea, we are going to propose two different ways of correcting the boundary values of the internal stages of the AMF-type W-methods in order to ensure that these stages are consistent with the boundary values. Through a number of numerical tests in several spatial dimensions with time-dependent Dirichlet boundary conditions, it will be illustrated that the new correction techniques to be proposed allow to retain the accuracy and the order of convergence that the AMF-type W-methods reach on PDEs with homogeneous Dirichlet boundary conditions.

The order reduction phenomenon associated with the time integration of initial boundary value problems with time-dependent boundary conditions has been addressed in a number of papers in the literature in the last decades. To mitigate such order reduction, several techniques have been proposed for linear and nonlinear problems with several numerical integrators, such as Runge–Kutta and Rosenbrock methods and exponential, Lie and Strang splitting methods, among others (see, e.g., [1,2,3] and the references therein). To the best of our knowledge, the correction techniques to be proposed in the current paper are new in the context of W-methods with the AMF approach.

The rest of the paper is divided into sections as follows. Firstly, in Sect. 2, we introduce the spatial discretization using a notation similar to the one introduced in [13] with the purpose of simplifying the treatment of finite differences on high-dimensional PDEs and to visualize the boundary values “in the same way as the interior points.” In fact, with this notation, the whole discretization can be seen as a system of differential-algebraic equations (DAE), so we explore this possibility in Sect. 3, and it will be seen that it does not avoid the order reduction of the splitting methods. Then, in Sect. 4, two different boundary corrections are proposed. Both of them expand the discretization to the borders, and, at the same time, they mitigate the order reduction when applying AMF-W-methods on PDEs (1), (2), and (3). This is experimentally seen in Sect. 5. In Sects. 6 and 7, these new techniques to avoid the order reduction are extended to the important case of parabolic PDEs with mixed derivatives.

Following the ideas given in [13], a space discretization is performed on a uniform spatial grid on \(\overline{\Omega }\).

Given N integers \((M_1,\dots ,M_{N-1},M_N)\), we consider on \(\overline{\Omega }\) the uniform spatial grid with \(M_i+2\) points at the \(x_i-\)direction

and discretize the derivatives in (1) at each spatial node \((x_{1,j_1},\dots ,x_{N,j_{N}})\in \Omega \) with finite differences. With this approach, we get approximations \(u(t,x_{1,j_1},\dots ,x_{N,j_{N}})\approx w_{j_1,\dots ,j_{N-1},j_N}(t),\) with \(w_{j_1,\dots ,j_{N-1},j_N}(t)\) satisfying the semi-discretized ODE system

where \(\textbf{x}_{r,j_r}\) denotes the mesh point \((x_{1,j_1},\dots ,x_{N,j_N}) \). Here, \(\Delta ^{(i)}w\) represents the finite differences approximation for the derivative \(\partial ^2_{x_i x_i}u \). It must be noticed that non-uniform spatial meshes can be also considered, but the notation would be more complicated and would not help to make this presentation clearer to read.

The handling of these differences becomes very cumbersome when the dimension is greater than 3, especially when we have to deal with the values on the boundary. In [13, Lemma 3.1], a map is defined that assigns an integer J to each multi-index \((j_1,\dots ,j_N)\) for all \(j_i\) running from an integer \(m_i\) to another \(M_i\ge m_i\). Moreover, the set of all possible integers J for a given spatial grid is split into two subsets. The first one, \(\textsf {Outer}\), groups the integers corresponding to the grid points on the boundaries where the Dirichlet conditions are imposed, and the second one, \(\textsf {Inner}\), stores the remaining integers, which are those where the solution of the PDE (1) is unknown.

Here, we suppose that time-dependent Dirichlet conditions are imposed in all the boundary of \(\Omega \). Then, we apply [13, Lemma 3.1] to define the two following maps, with \(M=\prod _{k=1}^N (M_k+2)\), \(\hat{M}=\prod _{k=1}^N M_k\),

Then, the set of integers \(\{0,1,\dots ,M-1\}\) can be separated as the union of the two disjoint sets

As a consequence, for each \(J\in \{0,1,\dots ,M-1\}\), \(\vartheta _0^{-1}(J)=(j_1,\dots ,j_N)\), we denote \(\textbf{x}_J:=(x_{1,j_1},\dots x_{N,j_N})\), that is, the corresponding mesh point, or \(\Omega _h=\{ \textbf{x}_J:\, 0\le J \le M-1 \}.\) If \(J\in \textsf {Outer}\), the solution \(u(t,\textbf{x}_J)\) is determined by the Dirichlet boundary conditions (2). The unknown in our problem is the value of \(u(t,\textbf{x}_J)\) for each \(J\in \textsf {Inner}\). Throughout this article, J will denote \(J=\vartheta _0(\textbf{j})\) for a unique \(\textbf{j}\in \mathcal{I}_N^{(0)}\).

Therefore, we can work with the exact solution of the PDE on the whole grid \(\Omega _h\),

with a notation that is clearer than the usual multi-index writing, and it is independent of the dimension N.

Then, the space discretization consists of looking for an approximation \(Y(t)=\left( Y_J(t)\right) _{J=0}^{M-1}\) where the derivative of each component \(Y_J(t)\) when \(J \in \textsf {Inner}\) must satisfy

where \(\left( d _i\right) _J:=d_i (\textbf{x}_J)\). On the other hand, when \(J\in \textsf {Outer}\), as usual, the boundary conditions are included in the discretization (8) by writing

In order to define the finite differences in (5), we need to take into account that for all \(J\in \textsf {Inner}\) and each \(E_i\) given in (6), \(J\pm E_i \in \{0,1,\dots ,M-1\}\), and when \(J\pm E_i \in \textsf {Inner}\), \(w_{\dots ,j_i\pm 1,\dots } (t)=Y_{J\pm E_i}.\)

In general, we can use several types of finite differences \(\Delta _J^{(i)}\), but, from now on, we will focus on second-order central differences on the uniform mesh \(\Omega _h\), which, for each \(J \in \textsf {Inner}\), can be written as

Writing with this notation the well-known property of these finite differences given in [11, p. 43], we have the following result which will be applied later in this paper.

Proposition 1

[11] If the solution u of the PDE (1)-(2)-(3) satisfies \(u\in C^4((0,t^*)\times {\Omega })\) with all partial derivatives of order less or equal to 4 continuous on \([0,t^*]\times \overline{\Omega }\), then, for each \(i=1,\dots ,N\) and for all \(J=\vartheta _0(j_1,\dots ,j_N)\) such that \(1\le j_i \le M_i\),

where U(t) is the vector given in (7).

In the sequel, to simplify the notation, we will denote \(\Delta _J^{(i)}=\Delta _J^{(i)}Y\) as long as it does not cause confusion. If the boundary conditions (2) are time-independent, i.e., \(Y_J=\beta (\textbf{x}_J),\,J\in \textsf {Outer}\), we can define \(Y_J'=0,\,J\in \textsf {Outer}\) and (8) can be written formally as the semi-discretized initial value problem (IVP)

where \(Y_0=(u_0(\textbf{x}_J))_{J=0}^{M-1}\) and \(\mathcal{F}(t,Y)=(\mathcal{F}_J(t,Y))_{J=0}^{M-1}\), whose components, for \(J\in \textsf {Inner}\) are

and \(\mathcal{F}_J(t,Y) = 0, \,\, \hbox {if } J\in \textsf {Outer}.\) Although it is not a proper semi-discretized IVP, since we need to involve the invariant (9), from the theoretical point of view, this is only a simpler way to write the usual semi-discretized IVP with these differences (see, for instance, [10]) when PDEs of high dimension are solved. This is the option used in [13] to integrate satisfactorily a model with spatial dimension six.

When time-dependent boundary conditions (2) are imposed, the system (8)-(9) can be seen as a differential-algebraic system (DAE) of index 1 of the form

where \(\mathcal{M}=(\mathcal{M}_{JK})_{J,K=0}^{M-1}\) is the matrix that has rank \(\hat{M}=\prod _{k=1}^N M_k\) with entries

with \(\mathcal{F}_J(t,Y)\) given in (11) when \(J\in \textsf {Inner}\) and, \(\mathcal{F}_J(t,Y)=Y_J-\beta (t,\textbf{x}_J)\) when \(J\in \textsf {Outer}\).

Lemma 1

For the DAE (12), the derivative function \(\mathcal{F}\) can be written as

where each matrix \(\mathcal{A}_i\), for \(1\le i\le N\), stores the corresponding terms to the discretization of the second derivative \(\partial ^2_{x_ix_i}u\). In detail, for \(i=1,\dots ,N,\) the entries of each matrix \(\mathcal{A}_i=((\mathcal{A}_i)_{JK})_{J,K=0}^{M-1}\in \mathbb {R}^{M\times M}\) are null except for the values

and the components of the non-homogeneous term \(G(t,Y)=(G_J(t,Y))_{J=0}^{M-1}\) are

Proof

The expressions for the finite differences give the components of the matrices \(\mathcal{A}_i\) for all \(1\le i\le N\). However, when \(J\in \textsf {Outer}\), we need to add the corresponding term \(Y_J\). It is obtained taking into account that

We are interested in applying AMF-W-methods that come from using the approximate matrix factorization (AMF) technique to reduce the computational cost of the linear systems involved in Rosenbrock methods. Then, let us first apply the latter on (12). From [9, p. 104 and 408], on (12) the step from \((t_n,Y_n)\) to \(t_{n+1}=t_n+\tau \) with time step size \(\tau \) with the \(s-\)stage Rosenbrock method \((\alpha ,\Gamma ,b)\), \(\alpha =(\alpha _{ij})_{i<j}, \Gamma =(\gamma _{ij})_{i\le j},\,\gamma _{ii}=\theta , c_r=\sum _{j=1}^{r-1}\alpha _{rj}\), results,

with advancing formula \(Y_{n+1}= Y_n+\sum _{r=1}^{s} b_rK_r.\) Then, applying Lemma 1, we need to solve the s linear systems

with

Taking into account that \(\mathcal{M}- \theta \tau (\mathcal{A}+I_M-\mathcal{M})= (I_M-\theta \tau \mathcal{A}) - (1+\theta \tau )(I_M-\mathcal{M}),\)

As a consequence, when \(J\in \textsf {Outer}\),

If we apply, as usual, the projection \((Y_n)_J=\beta (t_n,\textbf{x}_J)\), we obtain that for every \(i=1,\dots ,s\) and \(J\in \textsf {Outer}\),

Application of this to the one-stage Rosenbrock methods gives

In the case of the one-stage method of order 2, with coefficients \(c_1=0,\,b_1=1,\,\theta =1/2\) (see [10, p. 153]), this implies

resulting in a method of order 1 at the boundary values, which is a drawback.

In the case of the 2-stage Rosenbrock methods, the first stage always satisfies (18), so the second one results

Expanding \(\beta (t_n+\alpha _{21}\tau ,\textbf{x}_J)=\beta (t_n,\textbf{x}_J)+\alpha _{21}\tau \dot{\beta }(t_n,\textbf{x}_J)+\mathcal{O}(\tau ^2),\,\tau \rightarrow 0\),

If this is applied to the 2-stage third order Rosenbrock that satisfies the order conditions [10, p. 154],

it turns out that

and therefore suffers a severe order reduction at the boundary values.

As a consequence, we will try another approach, consisting in reducing the DAE system (12) to some ODE that can be solved numerically by applying AMF-W-methods with boundary corrections.\(\square \)

In order to explain clearly both options, we initially consider PDEs (1), (2) and (3) with a source term, that is, with

4.1 Explicit correction

In this option, we consider an extended IVP that involves the boundary points

where \(\mathcal{A}\) is the same matrix as in (17) (with \(\partial G/\partial Y=0\)), but the term G(t) is selected in such way that (9) is an invariant of this ODE system. Since \((\mathcal{A}Y)_J=0,\,\forall J\in \textsf {Outer}\), a natural choice for G(t) is then

where \(\dot{v}(t)\) means the derivative of a function v(t) regarding t.

Therefore, an IVP of type (10) is obtained where the derivative function can be split as

where the entries of each matrix \(\mathcal{A}_i\) are given in (14).

For the time integration of the resulting IVP, \(Y'=\mathcal{F}(t,Y),\,Y(0)=Y_0,\) we apply the \(s-\)stage AMF-W-method [6, p. A2915] (see also [5, 7, 10]), with time step size \(\tau >0\),

where, by denoting, \(A=(a_{rj})_{r>j}\), \(L=(\ell _{rj})_{r>j}\), \(\mathbbm {1}=(1,\dots ,1)^T\in \mathbb {R}^s\), \(\rho =(\rho _r)_{r=1}^s\), it holds that

If (24) is applied directly, the internal stages for each \(J\in \textsf {Outer}\), \(i=1,\dots ,N,\, r=1,\dots ,s,\) result

since \((\mathcal{A}_i K_r^{(i)})_J=0\). Denoting the stage vector as \(\textbf{K}=(K_1,\dots ,K_s)^T\) and

the equation (26) can be written as

Then, because of (25),

where \(\Phi =\Psi \rho =(\phi _r)_{r=1}^s\). More precisely, if we do not apply any boundary correction,

As a consequence, the advancing formula at the boundaries results

where \(\widetilde{\Psi }^T=(\widetilde{\psi }_{1},\dots ,\widetilde{\psi }_{s})=b^T\Psi ,\) \(\widetilde{\phi }=b^T \Phi .\)

If the AMF-W-method has order of consistency \(p\ge 2\), its coefficients satisfy \(b^T \Psi c^k=1/k, k=0,1,\dots ,p\) and \(b^T\Psi \rho =0\), so \((Y_{n+1})_J=\beta (t_{n+1},\textbf{x}_J)+\mathcal{O}(\tau ^{p+1})\ne \beta (t_{n+1},\textbf{x}_J),\) \(J\in \textsf {Outer}\), so it does not fulfil exactly the boundary condition at \(t=t_{n+1}\). The more natural way to meet these boundary conditions is to apply the projection

However, even with this projection, a severe order reduction appears in the \(l_{\infty }-\)norm. Roughly speaking, the main inconvenient is that the internal stages are not consistent with the boundary values.

Then, we are going to propose a technique to correct the values at the boundaries of the intermediate internal stages \(K_r^{(i)},\,i=1,\dots ,N-1,\) i.e., the values \((K_r^{(i)})_J,\,J\in \textsf {Outer}\), which is inspired by the one proposed for the Yanenko’s method in [10, p. 365].

In order to do that, we need to take into account that the linear systems of type \((I-\theta \tau \mathcal{A}_i)K=D\), involved in the splitting method, can be grouped into groups of \(M_i\) equations with \(M_i\) unknowns, in the following way.

For each multi-index \((\dots ,j_{i-1},j_{i+1},\dots )\) with \(N-1\) components, \(1\le j_k\le M_k,\,\forall k\ne i\), if \(J_k=\vartheta _0(\dots ,j_{i-1},k,j_{i+1},\dots )\),

where the values \(K_{J_0}\) and \(K_{J_{M_i+1}}\) on the right-hand side correspond to values of the stages at the boundary points. More precisely, \(J_0\) and \(J_{M_i+1}\) belong to

We will introduce the boundary correction through the linear systems (30). Following the idea of treating the boundary values as the inner ones, we know from Proposition 1 that under certain conditions on the regularity of the solution of the PDE, for each \(i=1,\dots ,N\) and U(t) as in (7),

Then, disregarding the \(\mathcal{O}(h_i^2)-\)term, each linear system \((I-\theta \tau \mathcal{A}_i)K_r^{(i)}=K_r^{(i-1)}\), \(i\ge 2\), in (24) can be extended to the boundary taking

If we write these relations starting from \(i=N\), bearing in mind that \(K_r=K_r^{(N)}\) and going backwards,

and so on. Disregarding the \(\mathcal{O}(\tau ^2)-\)terms and approximating the partial derivatives with the same differences applied on the inner points, we propose to correct the values of the intermediate stages at the boundary \(\textsf {Outer}_i,\) as

It is clear that with this correction, the values for all the iterations at the boundary depend only on \(\beta (t,\textbf{x})\). In Algorithm 1, a brief outline of the whole algorithm can be seen.

Explicit boundary correction of the AMF-W-method (24)

Full size image

Remark 1

When the reaction-diffusion PDE (1) is considered, where the non-homogeneous non-stiff term c depends on u, we must slightly adapt this algorithm. In this case, we must consider the extended IVP

with the matrices \(\mathcal{A}_i\) given in (14) and the vectors G and C with entries

Then, a \((N+1)-\)term splitting is applied \( \mathcal{F}(t,Y)= \sum _{i=1}^{N+1} \mathcal{F}_i(t,Y)\), with the same terms \(\mathcal{F}_i(t,Y),1\le i\le N,\) given in (23) and \(\mathcal{F}_{N+1}(t,Y)=C(t,Y)\), so the formulation of the \(s-\)stage AMF-W-method remains the same as in (24) except for the addition of the linear system

where \(\mathcal{A}_{N+1}= (\partial C/\partial Y)(t_n,Y_n).\) Since this matrix \(\mathcal{A}_{N+1}\) is diagonal, the latter linear system is solved directly and it does not require any boundary correction. Moreover, \((\mathcal{A}_{N+1}Y)_J=0,\forall J\in \textsf {Outer}, \forall Y\in \mathbb {R}^{M}\), so the explicit correction is applied in the same way as in Algorithm 1 only on the linear systems of type \((I-\theta \tau \mathcal{A}_i)K=D\) for \(i=1,\dots ,N\).

Remark 2

When \(N=1\) and no splitting is applied to (32), the formulation of (24) gives, in fact, the \(s-\)stage Rosenbrock method

with \(\widetilde{\mathcal{A}}=(\partial \mathcal{F}/\partial Y)(t_n,Y_n)\). In this case, the proposed explicit correction is related to the one introduced in [1, p. 6]. To understand this relation, it must be observed that, in that paper, the authors use the equivalent expression of Rosenbrock methods given in (16) (with \(\mathcal{M}\) being the identity matrix). If we denote \(\hat{\textbf{K}}=(\hat{K}_1,\dots ,\hat{K}_s)^T\) the internal stages there, they are related with the stages in (33) by

Then, if we apply our explicit correction, because of (27), we obtain that

On the other hand, to apply the correction proposed in [1] to the IVP (32), for \(1\le j\le 4\), the boundary values for the internal stages \(\hat{\textbf{K}}\) are chosen as

Therefore, expanding the Taylor’s series of \(\dot{B}_J(t_n \mathbbm {1}+\tau c)\) at \(t_n\) in (34), we get that \( \hat{\textbf{K}}_J=\left( K_n^{* [j]} \right) _J +\mathcal{O}(\tau ^{j+1}),\, \tau \rightarrow 0,\) at least up to \(j=4\). We must notice that one advantage of our approach is that only requires the computation of the time derivatives of order one and two of the boundary conditions, whereas the correction (34) must compute derivatives of higher order when \(j\ge 3\).

4.2 Implicit correction

Another approach is to extend the finite differences to the boundaries as long as they can be computed. On the PDE problem (1) with a source (19), it means to extend the second-order central differences to the boundary points, wherever it is possible, i.e., we define

and build an IVP of type \(Y'=\mathcal{A}Y+G(t)\) where (9) is an invariant of it. Because of (35), the diffusion matrices \(\{\mathcal{A}_i\}_{i=1}^N\) are no longer the same matrices as in (14). In this case, for each \(i=1,\dots ,N\), all the entries of \(\mathcal{A}_i\) are null but

On the other hand, clearly, \((G(t))_J=c(t,\textbf{x}_J)\) when \(J\in \textsf {Inner}\), but we need to select \((G(t))_J\) at the borders so that (9) is an invariant of the differential system, so we propose to take

where

and \(\Delta _J^{(i)}(\beta (t,\textbf{x}_J))\) is the same as in (35) when \(Y_J= \beta (t,\textbf{x}_J)\). From Proposition 1, if u is differentiable enough, for each \(J\in \textsf {Outer}\), \(\Delta _J^{(i)}(\beta (t,\textbf{x}_J))= \partial ^2_{x_ix_i} \beta (t,\textbf{x}_J)+ \mathcal{O}(h_i^2),\) when \(1\le j_i\le M_i\). Therefore, in some cases, the second-order derivatives of \(\beta (t,\textbf{x})\) can be used for correction instead of the differences.

As in the previous section, we integrate this IVP with the \(s-\)stage AMF-W-method (24). We must observe that now the linear systems of type \((I-\theta \tau \mathcal{A}_i)K=D\) involved can be decoupled in systems of dimension \(M_i\) but with many more groups than in the previous case.

For each multi-index \((\dots ,j_{i-1},j_{i+1},\dots )\) with \(N-1\) components and \(0\le j_k\le M_k+1,\,\forall k\ne i\), the following linear system similar to (30) is solved, denoting \(J_k=\vartheta _0(\dots ,j_{i-1},k,j_{i+1},\dots )\),

Moreover, as in the explicit correction, after applying the method at each step, the projection (29) is carried out.

It must be observed that, although the stages at the boundary points are computed with the same type of linear systems than the stages at the inner ones, this procedure is another correction of the values of the iterations \((K_r^{(i)})_{J\pm E_i}\) when \(J\pm E_i\in \textsf {Outer}\) and the linear systems (39) are solved for \(J\in \textsf {Inner}\).

For instance, when \(N=2\), \(\textsf {Outer}=\textsf {Outer}_1 \cup \textsf {Outer}_2 \cup \textsf {Corner}\) where \(\textsf {Corner}=\{ \vartheta _0(0,0),\vartheta _0(M_1+1,0),\vartheta _0(0,M_2+1),\vartheta _0(M_1+1,M_2+1)\}\), so in the case of the AMF-W-method of \(s=1\) stage we get, for \(J\in \textsf {Corner}\),

On the other hand,

The values of \((K_1^{(2)})_J\) for \(J\in \{ \vartheta _0(0,1),\vartheta _0(0,2),\dots ,\vartheta _0(0,M_2)\} \subset \textsf {Outer}_1\) are computed in the corresponding linear system (39) for \(i=2\) of dimension \(M_2\), where the right-hand side depends only on the function \(\beta (t,\textbf{x})\) and its derivatives. The same applies to the other \(M_2\) points \(J\in \{ \vartheta _0(M_1+1,1),\vartheta _0(M_1+1,2),\dots ,\vartheta _0(M_1+1,M_2)\} \subset \textsf {Outer}_1\). For the indexes in \(\textsf {Outer}_2\) we have that

whereas for the values of \((K_1^{(1)})_J\) for the two groups \(J\in \{ \vartheta _0(1,0),\vartheta _0(2,0),\dots \) \(\vartheta _0(M_1,0)\} \subset \textsf {Outer}_2\) and \(J\in \{ \vartheta _0(1,M_2+1),\vartheta _0(2,M_2+1),\dots ,\vartheta _0(M_1,M_2+1)\} \subset \textsf {Outer}_2\) are computed in the two corresponding linear systems (39) with \(i=1\), that in this case are of dimension \(M_1\), and the right-hand sides depend only on \(\beta (t,\textbf{x})\). Then, all the values of \((K_1^{(i)})_J\) when \(J\in \textsf {Outer}\) depend only on the function \(\beta (t,\textbf{x})\) given by the Dirichlet boundary conditions. The same holds for \(s\ge 2\) stages. The complete process is summarized in Algorithm 2.

Implicit boundary correction of the AMF-W-method (24)

Full size image

Remark 3

As in the Remark 1, the implicit correction must be adapted when we are dealing with reaction-diffusion PDEs (1). To extend the differences to the boundary points keeping (9) as an invariant, an IVP of type (32) is built where the entries of the matrices \(\mathcal{A}_i,\,i=1,\dots ,N,\) are given in (36), those of C(t, Y) are now \((C(t,Y))_J= c(t,x_J,Y_J)\), for all \(J=0,1,\dots ,M-1\), and

Then, as in Remark 1, a \((N+1)-\)term splitting is carried out, adding a diagonal linear system similar to the one given in the aforementioned remark.

In order to compare the two proposed boundary corrections with the methods without corrections, we will test two AMF-W-methods given in [5]. The first one is the one-stage method AMF-W1 without boundary corrections, with coefficients \(b_1=1,c_1=0, \rho _1=1,\theta =1/2\), and the second one is the two-stage method AMF-W2 with coefficients \(b_1=5/4, b_2=3/4, c_1=0, c_2=a_{21}=2/3, \ell _{21}=-4/3, \rho _1=1, \rho _2=-1/3\) and \(\theta =(3+\sqrt{3})/6\). It is known that these methods reach orders 2 and 3 respectively when time-independent boundary conditions are imposed. We will denote AMF-WiE and AMF-WiI the extension of the method AMF-Wi, i=1,2, with explicit boundary correction (Algorithm 1) and implicit correction (Algorithm 2) respectively.

5.1 Multi-dimensional linear problem

In order to check the behavior of the boundary corrections on multi-dimensional PDEs, we apply them on (1) with constant coefficients \(d_i(\textbf{x})=1,\forall i=1,\dots ,N,\) and non-homogeneous source term \(c(t,\textbf{x})\) selected in such a way that the PDE has the exact solution

taking the appropriate initial and Dirichlet time-dependent boundary conditions. Here, uniform spatial meshes (\(h_i=h,\,\forall i=1,\dots ,N\)) are considered.

In the following tables, the error of these methods at the interior points in both Euclidean (\(l_2\)) and uniform (\(l_{\infty }\)) norms, as well as their estimated PDE-order of convergence (when the time step size \(\tau =h_i=h,\,\forall i\)), is displayed. In Table 1, we can see clearly that the order of the standard AMF-W1 is reduced to order 1 in the uniform norm, whereas both boundary corrections in AMF-W1E and AMF-W1I reach order 2. The results for the 3-dimensional case, given in Table 2, are similar, being the order in \(l_{\infty }\)-norm for the standard AMF-W1 even less than 1. On the contrary, both corrected methods AMF-W1E and AMF-W1I recover the order 2. In Table 3, the results for dimension \(N=4\) are given only when the length of the spatial grid is greater or equal to 1/64 because of the limitation of capacity of our standard computers, but they suggest that the order 2 is recovered with both boundary corrections for any spatial dimension.

In Tables 4, 5 and 6, the results obtained for the two-stage methods are displayed. Clearly, without boundary correction, the AMF-W2 method suffers from a severe order reduction, from three to one in the 2-dimensional case and even less than 1 when the dimension N is greater than two. On the other hand, both boundary corrections recover the order three in \(l_2-\)norm as well as in \(l_{\infty }-\)norm, where the explicit one gives slightly smaller errors than the implicit ones.

5.2 Diffusion–reaction problem with traveling-wave solution

Let us consider now the 2D-PDE from [10]

that belongs to the class of PDEs (1) with a reaction term \(c(u)=\varepsilon ^{-1} u^2(1-u)\). The initial Dirichlet time-dependent boundary conditions are taken from the exact solution

This function represents a wave traveling diagonally across the domain, where the smaller the parameter \(\varepsilon >0\), the steeper the slope of the wave. As in the previous example, second-order central differences are applied to approximate the spatial derivatives on a uniform grid with \(h_1=h_2=h\).

Unlike the previous PDE, if we apply second-order differences, we will have a spatial error that will not allow us to see if the AMF-W2 method reaches order 3 with the boundary corrections. Since we know the exact solution (42), we can build a spatial discretization of order 4, having in mind that

by using,

Then, for the explicit boundary correction, taking into account Remark 1, we apply the AMF-W-methods on an IVP (32) with \(N=2\), with the matrices \(\mathcal{A}_1\) and \(\mathcal{A}_2\) given in (14), \((d_i)_J=\varepsilon ,\,i=1,2,\) and

In the case of the implicit boundary correction, proceeding as in Remark 3, we have an IVP similar to (32) where the matrices \(\mathcal{A}_1\) and \(\mathcal{A}_2\) are given now by (36) and \((C(Y))_J= \varepsilon ^{-1} Y_J^2(1-Y_J)\) for all \(J=0,\dots ,M-1\),

In both cases, as indicated in Remarks 1 and 3, to apply the AMF-W-methods (24), a three-term splitting must be carried out, taking \(Y'=\sum _{i=1}^3 \mathcal{F}_i(t,Y)\), with \(\mathcal{F}_1=\mathcal{A}_1 Y+G(t),\, \mathcal{F}_2=\mathcal{A}_2 Y,\, \mathcal{F}_3=C(Y).\)

In Tables 7 and 8, the performance on this PDE of the boundary corrections for the methods of one and two stages, respectively, is shown. As in the previous example, we present the results obtained in both norms when \(\tau =h\). Clearly, the boundary corrections recover the orders 2 and 3 respectively, in \(l_{\infty }-\)norm, for both methods, as in the previous example. In addition, significant improvement in the accuracy for all the step sizes is observed with regard to the uncorrected AMF-W-method.

Our next goal is to extend the proposed boundary corrections to the numerical solution of semi-linear parabolic PDEs with mixed derivatives on \([0,t^*]\times \Omega \),

with the same initial and time-dependent Dirichlet boundary conditions (3)-(2), where \(\partial _{x_j} u\equiv \frac{\partial u}{\partial x_j},\; \partial ^2_{x_i x_j} u\equiv \frac{\partial ^2 u}{\partial x_i \partial x_j},\) the reaction term c satisfies the same requirements as in (1), the coefficients \(a_i(\textbf{x})\) correspond to advection terms on each space variable and the diffusion terms are those corresponding to the coefficient matrix \(\mathbb {D}=(D_{ij}(\textbf{x}))_{i,j=1}^N\), where \(D_{ii}=d_i,\,\forall i,\) and \(D_{ij}=m_{ij}/2,\,\forall i\ne j\), which will be assumed to be symmetric and positive semi-definite for every \((t,\textbf{x})\in [0,t^*]\times \overline{\Omega }\).

Proceeding in a similar way as with the PDE (1), the uniform spatial grid (4) is considered and second-order central finite differences are applied on each interior point of the mesh obtaining the approximations \(Y_J(t) \approx u(t,\textbf{x}_J),\,J\in \textsf {Inner},\) that must satisfy the equations, for all \(J\in \textsf {Inner},\)

as well as the boundary conditions (9). Here, \(\Delta _J^{(i)}\) are the same as in (8), but

In this case, the part of the Jacobian of the derivative function corresponding to the mixed derivatives terms does not have a nice tridiagonal structure, so the application of methods as the AMF-W-method (24) for the time integration is no longer appropriate. In this case, we are going to apply another type of AMF-W-methods, the AMFR-W-methods that were introduced in [6], to mitigate the reduction in the convergence order suffered by AMF-W-methods when mixed derivatives are included in the PDE problem. It should be mentioned that AMFR-W-methods also suffer from order reduction on the boundary values in case of time-dependent boundary conditions; hence, we will apply the two boundary corrections presented in the previous sections to avoid such a reduction.

In order to simplify the presentation of the boundary corrections, we will assume (19) again. In the case of having reaction terms, one more term should be added to the splitting of the derivative function, as in Remarks 1 and 3.

6.1 Explicit correction for AMFR-W-methods

With mixed derivates and the explicit boundary correction, the spatial discretization with second-order central differences leads to an ODE system (20) with Jacobian \(\mathcal{A}=\sum _{i=0}^N \mathcal{A}_i\), where the components of the matrices \(\mathcal{A}_i\), for \(i=1,\dots ,N,\) are given by (14) and the non-homogeneous term G(t) by (21). The additional matrix \(\mathcal{A}_0\) stores the terms corresponding to the discretization of the mixed derivatives and advection terms and, for all \(Y\in \mathbb {R}^M\),

In this case, the derivative function \(\mathcal{F}\) is now split as \(\mathcal{F}(t,Y)=\mathcal{F}_0(t,Y)+ \sum _{i=1}^N \mathcal{F}_i(Y)\) with

To solve this IVP with splitting, we apply an \(s-\)stage AMFR-W-method, in which each time step from \(t_n\) to \(t_{n+1}=t_n+\tau \) is computed by ([6, p. A2916])

where \(\nu \ge N \kappa _N\,\theta \) for certain values \(\kappa _N\) that can be seen in [6, Table 2] and are selected for stability reasons.

In a similar way as in Sect. 4.1, we propose that the internal stages at the boundary points be corrected following the same procedure as for AMF-W-methods, which is summarized in Algorithm 1, on the two sets of AMF-W-systems for the computation of \(K_r^{(i)}\) and \(\hat{K}_r^{(i)}\). We need to take into account that now the parameter \(\theta \) in Algorithm 1 is \(\nu \), the term \(\theta \rho _r \tau ^2\dot{G}(t_n)\) is not added in the first allocation of \(K_r\), and that for the stages \(\hat{K}_r^{(i)}\) of the refinement, the initial definition of \(K_r\) must be changed by the expression of \(\hat{K}_r^{(0)}\) given in (47).

6.2 Implicit correction for AMFR-W-methods

As in Sect. 4.2, expanding the differences to the boundaries as in (35), we obtain an ODE system of type (20) with \(\mathcal{A}=\sum _{i=0}^N \mathcal{A}_i\), where, for each \(i=1,\dots ,N\), the entries of \(\mathcal{A}_i\) are given in (36). Due to the presence of mixed derivatives and advection in the PDE, the term \(\mathcal{A}_0Y\) appears, corresponding to the discretization of them, whose entries are

where the second-order differences \(\Delta _J^{(ik)}\) and \(\nabla _J^{(i)}\) are expanded as in (35). Moreover, the time-dependent term G(t) has the same form as the one given in (37) but now

Then, to correct the values of the stages at the boundaries, we propose a similar procedure to the one given in Algorithm 2, but adapted to the AMFR-W-method, as we proposed for the explicit correction. As a consequence, all the values of \((K_r^{(i)})_J, (\hat{K}_r^{(i)})_J,\) when \(J\in \textsf {Outer}\) depend only on the function \(\beta (t,\textbf{x})\), as we saw in Sect. 4.2.

We present the numerical results obtained by the methods in the class of AMFR-W-methods that are the refinement of the AMF-W-methods seen in Sect. 5. In other words, we denote AMFR-W1 the uncorrected one-stage method AMFR-W1 that is the refinement of AMF-W1 in Sect. 5, and so on. We have tested them on uniform spatial meshes \(h_i=h,\,i=1,\dots ,N\) on the PDE (44) with constant coefficients and mixed derivatives terms where

and with a source \(c(t,\textbf{x})\) and initial and boundary conditions selected so that the PDE has the exact solution (40).

In Tables 9, 10 and 11, the results obtained for the one-stage methods are displayed, where the columns store the errors and the estimated PDE-orders in the \(l_2\) and \(l_{\infty }\) norms, as for AMF-methods. It can be seen that without correction, similar results as in the case without mixed derivatives are obtained, that is, its order in the \(l_{\infty }-\)norm decays to one when \(N=2\) and even less when \(N\ge 3\). On the contrary, both proposed boundary corrections recover order two.

Furthermore, we have also tested the boundary corrections on the AMFR-W2 method of two stages and ODE-order three. Numerical experiments not presented in this paper indicate that when this method integrates the PDE problem (44)–(50) with solution (40) and time-independent boundary conditions (\(\kappa =0\)), it only reaches PDE-order 2.75 in \(l_2-\)norm and 2.5 in \(l_{\infty }-\)norm. The results obtained on problems with time-dependent boundary values without and with corrections can be seen in Tables 12, 13 and 14. There, it can be observed that the PDE-order in \(l_2-\)norm drops to 2, whereas it is less than one in \(l_{\infty }-\)norm. Clearly, when the two boundary corrections are applied, the order is recovered in both norms and for all spatial dimensions.

In this paper, we have introduced two different boundary corrections to mitigate the order reduction suffered by AMF-W-type methods when applied on PDEs with time-dependent boundary conditions, one explicit and the other implicit. Although both are novel, the former can be understood as an extension to the one applied to Yanenko’s method in [10]. Numerically, we have seen that on several linear problems, the order is recovered, both in \(l_2\) and \(l_{\infty }\) norms, obtaining a substantial improvement in accuracy compared to the uncorrected methods. The work in progress is to explain the mitigation of the order reduction with the approaches here considered. We also intend to extend these corrections to other types of time-dependent boundary conditions.