Constructive a priori error estimates for a full discrete approximation of periodic solutions for the heat equation

TL;DR

This paper develops constructive a priori error estimates for a fully discrete numerical method solving the heat equation with periodic boundary conditions, enhancing the understanding of solution accuracy.

Contribution

It introduces new a priori error estimates specifically tailored for full discrete approximations of periodic heat equation solutions.

Findings

Provides rigorous error bounds for the numerical scheme.

Demonstrates the effectiveness of the estimates through analysis.

Improves reliability of numerical simulations for periodic heat problems.

Abstract

We consider the constructive a priori error estimates for a full discrete numerical solution of the heat equation with time-periodic condition.

Equations

$${\Vert u\Vert }_{L^2{H}_0^1}^2\equiv {\Vert u\Vert }_{L^2((0,T);H_0^1(\Omega ))}^2:={\int }_0^T{\Vert u(t)\Vert }_{H_0^1(\Omega )}^{2}dt,$$

$$\frac{\partial u}{\partial t}-\nu \Delta u=f(x,t)$$

$$u(x,t)=0$$

$$u(x,0)=u(x,T)$$

$$\frac{\partial \varphi }{\partial s}-\nu \Delta \varphi =0$$

$$\varphi (x,s)=0$$

$$\varphi (x,0)=v(x)$$

$$\frac{\partial \psi }{\partial s}-\nu \Delta \psi =f(x,s)$$

$$\psi (x,s)=0$$

$$\psi (x,0)=0$$

$$\psi (t)={\int }_0^{t}E(t-s)f(s)ds.$$

$$u(t)$$

$${\Vert E(t)v\Vert }_X$$

$$u(0)=E(T)u(0)+\psi (T).$$

$$u(0)=(I-E(T){)}^{-1}\psi (T).$$

$${\Vert u(0)\Vert }_{L^2(\Omega )}\le (1-e^{-\nu {\lambda }_{1}T}{)}^{-1}\frac{C_p}{\sqrt{\nu }}{\Vert f\Vert }_{L^2{L}^2},$$

$${\Vert u(0)\Vert }_{L^2(\Omega )}\le (1-e^{-\nu {\lambda }_{1}T}{)}^{-1}\sqrt{T}{\Vert f\Vert }_{L^2{L}^2}.$$

$${\Vert \nabla \psi (T)\Vert }_{L^2(\Omega )}\le \frac{1}{\sqrt{\nu }}{\Vert f\Vert }_{L^2{L}^2},$$

$${\Vert \nabla u(0)\Vert }_{L^2(\Omega )}\le (1-e^{-\nu {\lambda }_{1}T}{)}^{-1}\frac{1}{\sqrt{\nu }}{\Vert f\Vert }_{L^2{L}^2}.$$

$${\Vert u_t\Vert }_{L^2{L}^2}$$

$${\Vert u(T)\Vert }_{L^2(\Omega )}^2+\nu {\Vert \nabla u\Vert }_{L^2{L}^2}^2$$

$${\Vert u_t\Vert }_{L^2(\Omega )}^2+\nu \frac{d}{dt}{\Vert \nabla u\Vert }_{L^2(\Omega {)}^d}^2\le {\Vert f\Vert }_{L^2(\Omega )}^2.$$

$$\frac{d}{dt}{\Vert u(t)\Vert }_{L^2(\Omega )}^2+\nu {\Vert \nabla u(t)\Vert }_{L^2(\Omega {)}^d}^2$$

$${\left(\nabla (u-P_h^{1}u),\nabla v_h\right)}_{L^2(\Omega )}=0\forall v_h\in S_h,$$

$${\Vert u-P_h^{1}u\Vert }_{H_0^1(\Omega )}$$

$${\Vert u-P_h^{1}u\Vert }_{L^2(\Omega )}$$

$${\left(\frac{\partial }{\partial t}(u-P_{h}u),v_h\right)}_{L^2(\Omega )}+\nu {\left(\nabla (u-P_{h}u),\nabla v_h\right)}_{L^2(\Omega {)}^d}=0\forall v_h\in S_h,t\in J,$$

$$(P_{h}u)(0)=(P_{h}u)(T).$$

$$\frac{\partial {\varphi }_h}{\partial s}-\nu {\Delta }_h{\varphi }_h=0$$

$${\varphi }_h(x,0)=v_h(x)$$

$$(({\varphi }_h{)}_t,{\eta }_h{)}_{L^2(\Omega )}+\nu (\nabla {\varphi }_h,\nabla {\eta }_h{)}_{L^2(\Omega )}=0\forall {\eta }_h\in S_h,t>0.$$

$$\frac{\partial {\psi }_h}{\partial s}-\nu {\Delta }_h{\psi }_h=P_h^{0}f$$

$${\psi }_h(x,0)=0$$

$$P_{h}u(t)=E_h(t)u_h(0)+{\int }_0^t{E}_h(t-s)P_h^{0}fds.$$

$${\Vert (I-E_h(T){)}^{-1}\Vert }_{\mathcal{L}(H_0^1)}\le (1-{\kappa }_1{)}^{-1}.$$

$${\Vert \nabla u_h(0)\Vert }_{L^2(\Omega )}\le (1-{\kappa }_1{)}^{-1}\frac{1}{\sqrt{\nu }}{\Vert f\Vert }_{L^2{L}^2}.$$

Full text

1 Introduction

Many works have been done concerning the error estimates for the approximate solutions of linear parabolic initial boundary value problems. Particularly, in [4], [2], they treated the time-periodic problems of the heat equation. On the other hand, recently, there are many results on the numerical enclosing the closed orbits corresponding to the periodic solutions by mainly using spectral techniques, [12],[3] etc., as part of the study in dynamical systems. In their works, the spectral properties for the simple operator restricted to the rectangular domains are effectively used. In the present paper, we consider the finite element approach instead the spectral method. Such a technique seems to be more complicated and the error estimates are not so easy compared with spectral method. But, there is no limit to the shape of the domain at all. The method we describe here basically extends the results of the previous paper [7] to the time-periodic problem of a heat equation.

In the followings, we use the time-dependent Sobolev spaces with associated norms of the form $L^{p}((0,t);X)$. For example, $u\in L^{2}((0,T);H^{1}_{0}(\Omega))$, then

[TABLE]

also use the notation such that $\left\|u\right\|_{L^{2}L^{2}}\equiv\left\|u\right\|_{L^{2}((0,T);L^{2}(\Omega))}$ for short and so on. For other notations and properties of function spaces, see e.g. [1], [11].

2 Problem and basic properties

In this section, we introduce the time-periodic problem and give the basic properties of the solution.

We consider the following heat equation with time-periodic condition:

[TABLE]

where $\nu$ is a positive constant, $J:=(0,T)\subset\mathbb{R}~{}(T<\infty)$ and $\Omega\subset\mathbb{R}^{d}$ $(d=1,2,3)$ a convex polygonal or polyhedral domains. Also we define $\displaystyle\Delta_{t}\equiv{\frac{\partial{}}{\partial{t}}}-\nu\Delta$ and assume that ${{f\in L^{2}((0,T);L^{2}(\Omega))}}\equiv L^{2}(\Omega\times J)$. On the existence and uniqueness of solution for (1), see e.g. [1], [11].

Now, for any $v\in L^{2}(\Omega)$ and $t>0$, we define the evolution operator $E(t):L^{2}(\Omega)\rightarrow L^{2}((0,t);H^{1}_{0}(\Omega))$ as a solution $\phi\in L^{2}((0,t);H^{1}_{0}(\Omega))$ of the following equation. Namely, $E(t)v\equiv\phi$ satisfies

[TABLE]

Next, consider the solution $\psi\in L^{2}((0,t);H^{1}_{0}(\Omega))$ satisfying the following parabolic problem with homogeneous initial condition

[TABLE]

Then note that by using the notation in semigroup theory, e.g., [8], we can rewrite (3) as follows:

[TABLE]

Taking notice that, using a solution $\phi$ of (2) for an appropriately chosen initial function $v=u(0)$ and $\psi$ in (3), the solution $u$ of (1) can be represented as $u(t)\equiv u(\cdot,t)=\phi(t)+\psi(t)$. Namely, we have

[TABLE]

Now, by the well known arguments using spectral theory in [1] or semigroup approches in [8], for the minimal eigenvalue ${\lambda}_{1}$ of $-\Delta$ on $\Omega$, it holds that for the spaces $X=L^{2}(\Omega)$ or $X=H_{0}^{1}(\Omega)$

[TABLE]

where $\left\|u\right\|_{H_{0}^{1}(\Omega)}\equiv\left\|\nabla u\right\|_{L^{2}(\Omega)}$. Then, from the periodic condition, we have by (4)

[TABLE]

Hence, from the contraction property of $E(T)$ due to the estimates (5), the invertibility of the operator $I-E(T)$ follows and the initial value $u(0)$ is determined by

[TABLE]

Furthermore, by the fact that $\psi$ is a solution of (3), it is readily seen that, by (5) and (7) (cf. in the proof of Lemma 4.2 of [7]):

[TABLE]

where $C_{p}$ is a Poincaré constant on $\Omega$. Also, if we use the fact that $\left\|\psi(T)\right\|_{L^{2}(\Omega)}\leq\sqrt{T}\left\|\psi_{t}\right\|_{L^{2}L^{2}}$ and the estimates $\left\|\psi_{t}\right\|_{L^{2}L^{2}}\leq\left\|f\right\|_{L^{2}L^{2}}$(Lemma 4.2 in [7]), we have another estimates as follows:

[TABLE]

By the similar arguments, from (5), (7) and the following estimates (cf. in the proof of Lemma 4.1 of [7])

[TABLE]

we have the bound for $\nabla u(0)$ as

[TABLE]

The following lemma can be similarly obtained.

Lemma 2.1**.**

For the solution $u$ of (1), it holds that

[TABLE]

Proof. As in the proof of Lemma 4.1 in [7] we have

[TABLE]

Integrating this on $J$, by taking notice of the periodic condition, yields (11).

Similarly, from the proof of Lemma 4.2 in [7] we get

[TABLE]

which proves (12) by combining with the estimates (9). $\Box$

3 Semidiscrete approximation

In the present section, we define the semidiscrete approximation by the finite element method and derive the constructive error estimates. These results play important and essential roles in the error estimates for a full-discretization of the problem (1).

Let $S_{h}\equiv S_{h}(\Omega)\subset H^{1}_{0}(\Omega)$ be a finite dimensional subspace in spatial direction with $\dim S_{h}=n$ and let $V_{k}^{1}\equiv V_{k}^{1}(J)\subset V^{1}(J)\equiv H^{1}(J)\cap\{u~{}|~{}u(0)=u(T)\}$ be a piecewise linear Lagrange type finite element space in time direction with $\dim S^{k}=m$. Also define V:=H^{1}\bigl{(}J;L^{2}(\Omega)\bigr{)}\cap L^{2}\bigl{(}J;H_{0}^{1}(\Omega)\bigr{)}\cap\{u~{}|~{}u(0)=u(T)\;\text{in}\;H_{0}^{1}(\Omega)\}.

Now, let $P_{h}^{1}:H_{0}^{1}(\Omega)\to S_{h}$ be an $H_{0}^{1}$-projection satisfying

[TABLE]

with the following assumptions on the approximation property:

[TABLE]

Here C$Y(\Omega):={\left\{{u\in L^{2}(\Omega)}~{};~{}{\Delta u\in L^{2}(\Omega)}\right\}}$.

Now, we define the semidiscrete projection P_{h}:V\to H^{1}\bigl{(}J;S_{h}(\Omega)\bigr{)}\equiv V^{1}\bigl{(}J;S_{h}(\Omega)\bigr{)}\cap\{v_{h}(0)=v_{h}(T)\} by the following weak form:

[TABLE]

Note that $P_{h}u$ implies the semidiscrete approximation of a solution $u$ for (1) with given function $f\in L^{2}(J;L^{2}(\Omega))$. Therefore, we denote $(P_{h}u)(t)$ by $u_{h}(t)$, i.e., $u_{h}\equiv P_{h}u$ in the below.

Next we consider the constructive error estimates for $P_{h}u$ defined by (18).

For any $v_{h}\in S_{h}$ and $t>0$, we define the semidiscrete evolutional operator $E_{h}(t):S_{h}\rightarrow S_{h}$ by the solution \phi_{h}\in H^{1}\bigl{(}(0,t);S_{h}(\Omega)\bigr{)} of the following equation. Namely, $E_{h}(t)v_{h}\equiv\phi_{h}$ corresponds to a semidiscretization of the solution $E(t)v\equiv\phi$ defined by (2).

[TABLE]

Here, $\Delta_{h}$ means the discretization of a weak Laplacian on $S_{h}$ and (19a) is equivalent to the following variational form:

[TABLE]

Similarly, as an semidiscretization for (3), we consider a solution \psi_{h}\in H^{1}\bigl{(}(0,t);S_{h}(\Omega)\bigr{)} of the following equation

[TABLE]

where $P_{h}^{0}f$ means the $L^{2}$-projection of $f$ to $S_{h}$. Also by using the similar symbol and arguments as in the previous section we get the following expression:

[TABLE]

Here, note that we can numerically compute the norm $\kappa_{1}:=\left\|E_{h}(T)\right\|_{\mathcal{L}(H^{1}_{0})}$ by matrix norm computations to confirm it is actually less than one, namely, contraction map on $S_{h}$. On the actual estimation of $\kappa_{1}$, see Remark 4.1 in the next section. And we can also compute the following inverse operator norm for $(I-E_{h}(T))^{-1}$

[TABLE]

Thus, from the definition and discrete analog to the previous section, we have $u_{h}(0)=(I-E_{h}(T))^{-1}\psi_{h}(T)$ and obtain the following estimates:

[TABLE]

Now, in order to get the error estimates for the semidisctrete approximation defined by (18) or equivalently by (22) for the problem (1), first we consider the constructive error estimates for the semidiscretization of the nonhomogeneous parabolic initial boundary value problem with initial condition $\xi_{0}\in H^{1}_{0}(\Omega)$ of the form :

[TABLE]

Let $\xi_{h}\in S_{h}$ be a semidiscrete approximation of (25) given by the following weak form:

[TABLE]

Here, $\zeta_{h}\in S_{h}$ is an appropriate approximation of $\xi_{0}$. Then we have the following estimates for solutions of (25) and (26).

Lemma 3.1**.**

[TABLE]

Proof. @These results are obtained by the similar arguments to that in the proofs for Lemma 4.1-4.4 in [7] with some additional considerations.

First, by the same argument to derive (13), we have

[TABLE]

which implies (27). Next, by the similar manner of getting (14) in the proof of Lemma 2.1, we have

[TABLE]

Thus integrating both sides in $t$ yields the estimates (28).

We now take $v_{h}:=(\xi_{h})_{t}$ for $t>0$ in (26a) and integrate it in $t$, we have

[TABLE]

which proves the assertion (29). Finally, the estimates (30) can be easily derived by the argument analogous to proving (28). $\Box$

Also, setting $\xi_{\perp}:=\xi-\xi_{h}$, we obtain the following two kinds of error estimates, which are obtained similar arguments in the proof of Theorem 4.6 in [7].

Theorem 3.2**.**

The following estimates for $\xi_{\perp}:=\xi-\xi_{h}$ hold:

[TABLE]

also $L^{2}$-estimates at $t=T$,

[TABLE]

Proof. Applying the same arguments in the proof of Theorem 4.6 in [7], we have

[TABLE]

Integrating this on $J$, from (31) and (32), we get

[TABLE]

which yields the desired conclusions (33) and (34). $\Box$

4 Full-discrete approximation and error estimates

In this section, we define the full-discrete approximation of solutions for the problem (1) by using an interpolation procedure in time direction for the spatial discretized solution. We also show a computational scheme for this full discretization by the effective use of the fundamental matrix for an ODE system corresponding to semidiscretized problem. The constructive and optimal order $H^{1}$ and $L^{2}$ error estimates are established, which are main results in the present paper.

4.1 A full discretizaion scheme

Now, defining the interpolation operator $\Pi^{k}:$ $~{}V^{1}(J)$ $\to$ $V_{k}^{1}$ in time direction by

[TABLE]

we define the full discrete projection P_{h}^{k}:V\to V_{k}^{1}\bigl{(}J;S_{h}(\Omega)\bigr{)}\equiv S_{h}\otimes V_{k}^{1} as

[TABLE]

which corresponds to the full discretization of (1).

In order to present the actual computation procedure of the above full discretization scheme, we first consider a representation of the semidiscretization defined in (18). Let $\{\phi_{i}\}_{i=1}^{n}$ be a basis of $S_{h}$ and define the $n\times n$ matrices $L_{\phi}$, $D_{\phi}$ by

[TABLE]

respectively. Since they are symmetric and positive definite, we get the Cholesky decomposition as $L_{\phi}=L_{\phi}^{1/2}L_{\phi}^{T/2}$ and $D_{\phi}=D_{\phi}^{1/2}D_{\phi}^{T/2}$, respectively. Also note that there exists a vector valued function $\vec{u}_{h}\in V^{1}(J)^{n}$ satisfying

[TABLE]

where $\Phi(x)\equiv(\phi_{1},\cdots,\phi_{n})^{T}$.

Thus by using $\vec{u}_{h}$ Cthe semidiscretization (18) is equivalently presented as ODEs:

[TABLE]

where $\tilde{f}=(\tilde{f}_{i})\in\mathbb{R}^{n}$ with $\tilde{f}_{i}=\left(f,\phi_{i}\right)_{L^{2}(\Omega)}.$ For simplicity we denote as ${\vec{b}}(t)\equiv L_{\phi}^{-1}\tilde{f}(t)$. Then note that using the fundamental matrix $\Theta(t)=\exp(-\nu L_{\phi}^{-1}D_{\phi}t)$ of the equation (37a), we can represent (37) as

[TABLE]

Therefore, assuming that the invertibility of $\left(I-\Theta(T)\right)$ Cfrom (38a) Cwe have

[TABLE]

which yields the following expression of the solution of (38) F

[TABLE]

Hence, we obtain

[TABLE]

Thus the full discrete approximation $P_{h}^{k}u\equiv\Pi^{k}P_{h}u$ for the solution $u$ of (1) can be numerically computed by using this procedure.

Remark 4.1:

For any $v_{h}\in S_{h}$, using the definition (36), by some simple consideration on the $H^{1}_{0}$ norm for the element $E_{h}(T)v_{h}\in S_{h}$, we have readily seen that

[TABLE]

where $||\cdot||_{2}$ means the matrix 2-norm. This immediately yields the estimate of $\kappa_{1}$ in (23).

4.2 $H^{1}$ error estimates

In this subsecton, we present an error estimate in the $L^{2}H^{1}_{0}$ sense on $\Omega\times J$ for the full discretization (35). Denoting again the semidiscrete projection $P_{h}u$ defined in (18) as $P_{h}u\equiv u_{h}$, the semidiscrete approximation $u_{h}$ for (1) is written by

[TABLE]

In order to obtain the desired estimates, we use the following decomposition

[TABLE]

The second term of the above is estimated by using the standard interpolation estimates, e.g., [10], we have from (29) and (24)

[TABLE]

Furthermore, using an inverse estimation constant $C_{inv}(h)$, which makes possible to bound the $H^{1}$ norm by the $L^{2}$ norm in $S_{h}$, we get

[TABLE]

Note that using the definition of the operator $E_{h}(t)$, we have by (40)

[TABLE]

Therefore, using $\psi(t)$ defined by (3), we have

[TABLE]

which implies

[TABLE]

Note that, for any $t\in J$, setting

[TABLE]

[TABLE]

then $\xi$ and $\xi_{h}$ are solutions corresponding to (25) and (26), respectively.

Hence, setting $\xi_{\perp}:=\xi-\xi_{h}$, the right-hand side of (45) coincides with $\xi_{\perp}(T)$. Therefore, we have

[TABLE]

By the argument in the section 2, we have the following estimates

[TABLE]

Next, applying the error estimates (34) in Theorem 3.2 with taking notice of $\xi_{0}=\zeta_{h}$, by using (24) we have

[TABLE]

Therefore, from (46)-(48), we obtain

[TABLE]

where

[TABLE]

On the other hand, we have by (33) in Theorem 3.2

[TABLE]

Thus, from the estimates (10), (24) and (49), we obtain the following estimation for the semidiscrete solution:

[TABLE]

where we set as

[TABLE]

Combining (43) and (51) with (41), we have the following desired $H^{1}$ error estimates.

Theorem 4.1**.**

Let $P_{h}^{k}u$ be a full-discrete approximation defined by (35) for the periodic solution $u$ of the heat equation (1). Then, it holds that

[TABLE]

Here, the constant $K_{2}$ is defined in (51).

4.3 $L^{2}$ error estimates

In this subsection, we consider the error estimates in the $L^{2}L^{2}$ sense for the full-discrete approximation $P_{h}^{k}u$, which enable us higher order estimates than the $L^{2}H^{1}$ error bound in Theorem 4.1. As in the previous subsection, we use a semidiscrete approximation $u_{h}$ with decomposition (41). Note that, if we take as $k\approx h^{2}$, then by applying the $L^{2}$ estimates (42), we immediately obtain $O(h^{2})$ estimates for the latter term in (41). Hence, it suffices to derive the $O(h^{2})$ estimates for the former part.

Theorem 4.2**.**

It holds that

[TABLE]

where $K_{2}$ is the same constant defined in the estimates (51).

Proof. @For any function $g\in L^{2}(Q)$, where $Q\equiv\Omega\times J$, let $v$ be a solution of (1) with the right-hand side $g(T-t)\equiv g(\cdot,T-t)$. Here, $t$ is a variable such that $t\in J$. Then $v$ satisfies the following weak form:

[TABLE]

Particularly, taking $w=u-u_{h}$ in (54) and transform the variable as $t\rightarrow T-s$, we have

[TABLE]

Integrating both sides of the above in $s$ on $(0,T)$ yields that

[TABLE]

Taking notice of the periodic condition, observe that

[TABLE]

Therefore, by the definition of $u_{h}$ and (55) we have for any $v_{h}(s)\in S_{h}$

[TABLE]

Moreover, by the similar derivation process of (42) using (29) in the previous subsection and Lemma 2.1, we obtain

[TABLE]

Furthermore, for any $t\in(0,T)$, taking $v_{h}(t):=P_{h}^{1}v(t)$ to apply the approximation properties (16) and (17), by considering the estimates in Lemma 2.1, we have

[TABLE]

and

[TABLE]

Therefore, combining (57)-(59) with (51), we have the estimates

[TABLE]

which proves the theorem by (42). $\Box$

5 Numerical examples

In this section, we show several numerical examples which confirm us the optimal rate of convergence. We used the interval arithmetic toolbox INTLAB 11 [9] with MATLAB R2012a on an Intel Xeon W2155 (3.30 GHz) with CentOS 7.4.

Here, we only consider $d=1$, $\Omega=(0,1)$ and $J=(0,1)$, then the lower bound of eigenvalue of $-\Delta$ on $\Omega$ can be taken as $\lambda_{1}=\pi^{2}$. Furthermore, we set $f$ to be the problem (1) have the exact solution $u(x,t)=\sin(2\pi x)\sin(2\pi t+\beta)$. Here, $\beta$ is a given constant. Since the exact solutions are known, the upper bounds of the exact errors for approximate solutions can be validated in the a posteriori sense.

We used the finite dimensional subspaces $S_{h}$ and $V_{k}^{1}$ spanned by piecewise linear basis functions with uniform mesh size $h$ and $k$, respectively. Therefore, the constants can be taken as $C_{\Omega}(h)=h/\pi$, $C_{J}(k)=k/\pi$, ${{C_{inv}(h)}}=\sqrt{12}/h$, and $C_{p}=1/\pi$, respectively. We set $k=h^{2}$ then Theorem 4.1 and 4.2 are $O(h)$ and $O(h^{2})$ error estimates . In Figure 1-2, the a priori error estimates and the exact errors of this example are shown. These Figures show the estimates presented in Theorem 4.1-4.2 give the optimal order estimates.