On the a posteriori estimates for inverse operators of linear parabolic equations with applications to the numerical enclosure of solutions for nonlinear problems

We consider the guaranteed a posteriori estimates for the inverse parabolic operators with homogeneous initial-boundary conditions. Our estimation technique uses a full-discrete numerical scheme, which is based on the Galerkin method with an interpolation in time by using the fundamental solution for semidiscretization in space. In our technique, the constructive a priori error estimates for a full discretization of solutions for the heat equation play an essential role. Combining these estimates with an argument for the discretized inverse operator and a contraction property of the Newton-type formulation, we derive an a posteriori estimate of the norm for the infinite-dimensional operator. In numerical examples, we show that the proposed method should be more efficient than the existing method. Moreover, as an application, we give some prototype results for numerical verification of solutions of nonlinear parabolic problems, which confirm the actual usefulness of our technique.


Introduction
Setting L t := ∂ ∂t − ν + b · ∇ + c, for f ∈ L 2 J ; L 2 (Ω) , consider the following linear parabolic partial differential equations (PDEs) with homogeneous initial and boundary conditions: u(x, 0) = 0, in , where Ω ⊂ R d , (d ∈ {1, 2, 3}) is a bounded polygonal or polyhedral domain, J := (0, T ) ⊂ R, (T < ∞) is a bounded interval, ν is a positive constant, b ∈ L ∞ J ; L ∞ (Ω) d , and c ∈ L ∞ J ; L ∞ (Ω) . As is well known, for any f ∈ L 2 J ; L 2 (Ω) , there exists a unique weak solution u ∈ L 2 J ; H 1 0 (Ω) to the problem (1). Denoting the solution operator of (1) by L −1 t , it is a bounded linear operator from L 2 J ; L 2 (Ω) to L 2 J ; H 1 0 (Ω) . The main aim of this paper is to obtain the concrete value C L 2 L 2 ,L 2 H 1 0 > 0 satisfying the following estimates: The constant C L 2 L 2 ,L 2 H 1 0 plays an important role in the verification of solutions for the initial-boundary-value problems for the nonlinear parabolic PDEs, and we usually need to estimate it as small as possible. The concrete value C L 2 L 2 ,L 2 H 1 0 > 0 satisfying (2) can be calculated by the Gronwall inequality or other theoretical considerations (e.g., [16]), which we call the "a priori estimates." However, in general, C L 2 L 2 ,L 2 H 1 0 obtained by such a priori estimates is exponentially dependent on the length of the time interval J unless the corresponding elliptic part of the operator L t is coercive [4,5]. Thus a priori estimates often lead to an overestimate for the norm of L −1 t , which yields worse results for some purposes.
In order to overcome this difficulty, we proposed a method to calculate C L 2 L 2 ,L 2 H 1 0 by numerical computation with guaranteed accuracy in [10], which we called "a posteriori estimates." The method is based on combining the a priori error estimates for a semidiscretization with the a priori estimates for the ordinary differential equations (ODEs) in time. It has proven to be more efficient than the existing a priori method; some numerical examples show that this a posteriori method can remove the exponential dependency on the time interval J . However, it has a very large computational cost, because the semidiscretization of (1) causes stiff ODEs that require a very small step size. Also, it is not clear what time-space ratio to use in the discretization process.
In this paper, we propose a new a posteriori method with a fully discretized Newtontype operator, which uses the Galerkin approximation in the space direction and the Lagrange-type interpolation in the time direction. In the case of the simple heat equations, some fundamental properties (e.g., the stability and a priori error estimates) for this full-discretization scheme have already been obtained in [11]. In the desired estimation of the inverse operator norm L −1 , the matrix norm estimates corresponding to the discretized inverse operator and the constructive error analysis for the simple heat equations are important and essential. By constructive analysis, we can also guess an appropriate time-space ratio prior to the actual computation. Moreover, by using numerical examples, we will show that the proposed method succeeds in obtaining a posteriori estimates with less computational cost than the previous method in [10]. This means that the present method is very robust compared with the previous one. The contents of this paper are as follows: In Sect. 2, we introduce some function spaces, operators, and other notation. In Sect. 3, we introduce the results of stability and a priori error estimates for the full-discretization scheme for the simple heat equations, which were obtained in [11]. In Sect. 4, we consider the approximate quasi-Newton operator that corresponds to the full-discretization scheme for problem (1). In Sect. 5, we derive the new a posteriori estimates of (2) by combining the results in Sect. 3 with the property of the approximate quasi-Newton operator defined in the previous section. In Sect. 6, we compare the computed values for C L 2 L 2 ,L 2 H 1 0 by three methods, namely, the a priori method, the a posteriori estimates in [10], and the new a posteriori method obtained in Sect. 5. In this section we also show some prototype results of the numerical enclosure of solutions for nonlinear parabolic problems as an application of our method.

Notation
In this section, we introduce some function spaces, operators, and other notation. Let L 2 (Ω) and H 1 (Ω) be the usual Lebesgue and Sobolev spaces on Ω, respectively, and define the natural inner product of u, v in We will sometimes refer to the following Sobolev inequality on H 1 0 (Ω). Namely, for a suitable constant p ≥ 1, which is dependent on the dimension of Ω, there exists a constant C s, p > 0 such that When p = 2, (3) is called the Poincaré inequality. Let We denote the function space L 2 J ; L 2 (Ω) as L 2 L 2 , for short. Let L 2 J ; H 1 0 (Ω) be a subspace of L 2 L 2 defined by Then, Then, Moreover, we define the partial differential operator t : Then, the inverse of t exists (e.g., [3]), and we denote it by −1 From the compactness of the embedding Let S h (Ω) be a finite-dimensional subspace of H 1 0 (Ω) dependent on the discretization parameter h. For example, S h (Ω) is considered to be a finite element space with mesh size h. Let n be the number of degrees of freedom of S h (Ω), and let {φ i } n i=1 ⊂ H 1 0 (Ω) be the basis functions of S h (Ω). Moreover, we denote a vector of the basis functions of S h (Ω) by φ := (φ 1 , . . . , φ n ) T . We also assume the inverse estimates on S h (Ω) like as follows: We need the following assumptions as the a priori error estimates for P 1 h .

Assumption 2.2 There exists a positive constant C
For example, if Ω is a bounded open interval in R, and S h (Ω) is the P1 finite element space, then Assumption 2.2 is realized as C Ω (h) = h π , where h is the mesh size (see e.g., [1,7]).
Let V 1 k (J ) be a finite-dimensional subspace of V 1 (J ) dependent on the discretization parameter k. For example, V 1 k (J ) is considered to be a finite element space with mesh size (time step size) k. Let m be the number of degrees of freedom for V 1 k (J ), We need the following assumption as the a priori error estimate for Π k .

Assumption 2.3 There exists a positive constant C J
For example, if V 1 k (J ) is the P1 finite element space, then Assumption 2.3 is realized by C J (k) = k π (see e.g., [15,Theorem 2.4]). Let V 1 J ; S h (Ω) and V 1 k J ; S h (Ω) be the semidiscretization and the fulldiscretization subspaces of V , respectively. We now define the semidiscretization operator P h : V → V 1 J ; S h (Ω) by the following weak form for any u ∈ V Then the full-discretization operator P h,k : is defined as the composition of P h and Π k , that is, by P h,k := Π k P h .

Constructive a priori error estimates
In this section, we introduce some results for the stability of, and a priori error estimates for, the full-discretization operator P h,k . Since the results of this section are given in [11], we omit the proofs.
Moreover, if V 1 k (J ) is the P1 finite element space then we have the following estimates: Since the full-discretization scheme proposed in [6,9] has no V 1 L 2 stability, we can say that the present full-discretized approximation has better properties, in an analytical and practical sense.
Finally, we introduce the constructive a priori error estimates for P h,k . where

Discretized quasi-Newton scheme
In this section, we consider a full-discretized approximation scheme for solutions of (1) by using a quasi-Newton operator. Since the full-discretization scheme in this paper uses interpolation in time, its computational method is somewhat complicated. However, it enables us to get an efficient and accurate estimation of the inverse operator norm in (2), as well as the verified computation of solutions to nonlinear problems. We first describe an easy, but an important operation of matrix-vector multiplication.
We call this transformation a "row-major matrix-vector transformation".

Definition 4.2
Let M be an m 1 -by-m 2 matrix. Then, we define the block diagonal matrix (I n ⊗ M) as follows: Here, I n is the n-by-n identity matrix, and the operator ⊗ denotes the Kronecker product.
From these definitions, we have the following lemma.

Lemma 4.3 For an arbitrary n-by-m matrix M and m-dimensional vector x, the following equality holds
: Proof The elements of M x are calculated by On the other hand, the elements of I n ⊗ x T vec (M) are calculated by Therefore, the corresponding components coincide with each other.
Next, we consider the quasi-Newton operator of (1) and its full-discretization. Let A be an integral operator defined by A : Then, the differential operator of the left-hand side of (1a) can be represented as L t = t (I − A), where I denotes the identity operator on D( t ). We define the quasi-Newton operator as the inverse of I − A, i.e., (I − A) −1 : L 2 H 1 0 → L 2 H 1 0 . We now define the symmetric and positive definite matrices L φ and D φ ∈ R n×n by φ be the Cholesky factors of L φ and D φ , respectively, i.e., the following equalities hold φ are lower triangular matrices, and L T /2 φ and D T /2 φ are those matrices transposed. Let L ψ ∈ R m×m be the symmetric and positive-definite matrix whose elements are defined by L ψ,i, j := ψ j , ψ i L 2 (J ) . We define Z φ ∈ L ∞ (J ) n×n as the matrix function on J whose elements are defined by For any i ∈ {1, . . . , m}, we define the matricesG Moreover, we define G φ,ψ ∈ R nm×nm as G φ,ψ := I nm −G φ,ψ . We obtain the Theorem 4.4 as a full-discretization scheme of the quasi-Newton operator.
For each v h ∈ S h (Ω), and almost everywhere t ∈ J , from the definition of P h and the operator A, we have From the arbitrariness of v h ∈ S h (Ω), the variational equation (20) is equivalent to the following system of first-order linear ODEs with homogeneous initial conditions Since (21) is an initial-value problem for an ODE system with constant coefficients, by using its fundamental matrix, w can be presented as where we have used (17) to make the deformation from (22) to (23). And, from (18), we have Thus, from (23), we obtain the following relation between U and w: Note that, for any nodal points t i , we have Since we assume that V 1 k (J ) is the finite element space constituted by the Lagrange elements, ψ j (t i ) = δ j,i is satisfied, where δ j,i denotes the Kronecker delta. Therefore, we get From the arbitrariness of x and i, the variational equation (25) is equivalent to the following simultaneous linear equations: Namely, we have Therefore, from the arbitrariness of f h,k , the nonsingularity of I nm −G φ,ψ follows. The converse of this proposition is easily obtained by reversing the discussion.
When we apply the proposed a posteriori estimates, it is necessary to confirm that G φ,ψ is nonsingular, which will be able to verify by validated computations such as [14]. Therefore, in what follows, we always assume the nonsingularity of G φ,ψ . Moreover, we define the linear operator by the solution of (19). We call this operator a "fully discretized quasi-Newton operator".

A posteriori estimates
In this section, we derive a new a posteriori estimate to obtain C L 2 L 2 ,L 2 H 1 0 , which satisfies (2) by using the fully discretized quasi-Newton operator.
First, we describe a method to calculate the norm of the elements in the fulldiscretization space. Let K φ,ψ be a matrix in R nm×nm defined by From the symmetric positive definiteness of D φ and L ψ , it is readily seen that K φ,ψ is also symmetric positive definite. Therefore, K φ,ψ is Cholesky decomposable such Similarly, we define the matrix L φ,ψ in R nm×nm as L φ,ψ := L φ ⊗ L ψ .

Lemma 5.1 For an element u h,k
, taking U ∈ R n×m such that u h,k = φ T U ψ, then the following equalities hold where | · | denotes the Euclidean norm of a vector.
Proof Since the proofs of (27) and (28) are almost the same, we will prove only (27). From (17), we have which proves equation (27).
, where · 2 denotes the matrix two-norm. The following theorem for M φ,ψ holds.

Theorem 5.2 It holds that
Proof This completes the proof.
Let C 0 and C 1 be the nonnegative constants defined by respectively. Moreover, we define the constant κ φ,ψ as follows:

Theorem 5.3 Assume that
Then under the same assumptions as in Theorem 3.2, we have the following constructive a posteriori estimates Proof For any f ∈ L 2 J ; L 2 (Ω) , we set u := L −1 t f ∈ D( t ). Then we make the following decomposition of (1) into two parts, e.g., the finite-and infinite-dimensional parts, using the projection P h,k . Namely, in the space L 2 J ; H 1 0 (Ω) , using the following equivalency we have the decomposition: We set u ⊥ := u − P h,k u for short. From (34a), using the definition of the operator A, we have Therefore, from (29) and (11), we have the following estimates, From the definition of C 0 , we have By calculating the L 2 L 2 norm of (34b) using (14), we have which yields Thus (35) is estimated as Setting nonnegative constants R 1,1 , R 1,2 , and b 1 as follows: is rewritten as On the other hand, by considering the L 2 H 1 0 norm of (34b), from (13) we have From (36), we obtain We set nonnegative constants R 2,1 , R 2,2 , and b 2 as follows: where we note that the positivity of R 2,2 follows by the condition (31). Thus (39) can be rewritten as From (38) and (40), we have the following simultaneous inequalities, By assumption (31), we obtain Therefore, the simultaneous inequalities can be solved as follows: ⎛ Finally, from (41), we have which proves the desired estimates.

Numerical examples
In this section, we show several rigorous numerical results for C L 2 L 2 ,L 2 H 1 0 satisfying (2) for test problems by three kinds of methods, namely, a priori estimates (the Gronwall inequality), a posteriori estimates proposed in [10], and the new method in Theorem 5.3. Moreover, we also show several rigorous error bounds of the numerical solutions for the nonlinear parabolic equations as an application of the estimates of (2).
We considered the norm estimates for an inverse operator of the following L t : that is, b = 0 and c = −2u k h in (2). Here, u k h is assumed to be an approximate solution of the following nonlinear parabolic problem: Therefore, (42) becomes a linearized operator of (43) at u k h . We only considered one-space-dimensional case (d = 1) with Ω = (0, 1). Furthermore, the function f was chosen so that the problem (43) had the following exact solutions: Note that Example 1.1 and Example 2.1 are studied in [10]. In each example, the function u k h was computed as an approximation of the corresponding u by using a piecewise-cubic Hermite interpolation in the space direction with a piecewise-linear interpolation in the time direction. Therefore, u k h belongs to V 1 J ; H 1 0 (Ω) ∩ H 2 (Ω) . We used the finite-dimensional spaces S h (Ω) and V 1 k (J ), spanned by piecewise linear functions with uniform mesh size h and k, respectively, so that they satisfied k = h 2 . Then, it was seen that the constants in previous sections could be taken as C Ω (h) = h/π , C inv (h) = √ 12/ h, C J (k) = k/π = h 2 /π , and C s,2 = 1/π , respectively. Moreover, we had T (Example 1.1 and 1.2) 2 (Example 2.1 and 2.2).

A posteriori estimates of the inverse parabolic operator
We now present the results computed for C L 2 L 2 ,L 2 H 1 0 by using three kinds of method. Here, we used the following a priori estimate, which comes from the Gronwall inequality In this subsection, we will refer to the a posteriori estimates studied in [10] and our present estimates (32) as "a posteriori estimate I" and "a posteriori estimate II," respectively. To compute a posteriori estimate I, we used the same parameters as in [10], i.e., (n, m) = (5, 700 · T 2 ) for Example 1.  Figures 1, 2 show the values of C L 2 L 2 ,L 2 H 1 0 for Example 1.1-1.2, plotted out on loglinear coordinates. For T > 1, the values of the proposed estimates are smaller than the other estimates. The two kinds of a posteriori estimates require the validated upper bound for the matrix two-norm of the corresponding unsymmetric dense matrices (e.g., M φ,ψ ), and most of the computational costs is due to this task. In Example 1.1, for T = 2, a posteriori estimate I a matrix of size 14000, but in a posteriori estimate II, we can attain our purpose with a matrix of size 896 for h = 1/8, and 7680 for h = 1/16. This fact shows, in the case of a posteriori estimate II, that it is not necessary to take special account of the stiff property of the ODEs coming from the semidiscretization. for Example 2.1-2.2 (log-linear coordinates). For T > 1/2, the values of the proposed estimates with h = 1/16 are smaller than the other estimates. In Example 2.1, for T = 2, a posteriori estimate I requires a matrix of size 8000, but a posteriori estimate II requires only one of size 896 for h = 1/8 and size 7680 for h = 1/16. It is notable that the results of the proposed estimates show no exponential dependency for T . On the other hand, due to the stiffness of the corresponding ODEs, we were not successful in computing the inverse operator of a posteriori estimate I, except for the case where T was very small.

Verification results for solutions of nonlinear parabolic equations
Applying the estimates (2), we implemented a numerical verification method to prove the existence of solutions for the nonlinear parabolic problems. As a prototype application, we considered the nonlinear parabolic initial-boundary-value problems of the form (43). In a similar way as in [8] for the elliptic case, we defined the fixed-point equation for a compact operator, which is equivalent to (43) with the Newton-type residual form, and derived a verification condition by applying the Schauder fixedpoint theorem.
First, we considered the following residual equation for (43): where Note that if the approximate solution u k h is close to the exact solution of (43), then w ≈ 0, ε ≈ 0, and g(w) ≈ 0. Thus (44) can be rewritten as the following fixed-point equation of the compact map F: Next, for any positive constants α and β, we define the candidate set W α,β as W α,β := w ∈ V ; w L 2 H 1 0 ≤ α, w V 1 L 2 ≤ β .
From the Schauder fixed-point theorem, noting that the continuity of the map F in the space L 2 H 1 0 , if the set W α,β satisfies then a fixed point of (45) exists in the set W α,β , where W α,β stands for the closure of the set W α,β in L 2 H 1 0 . Now, for any w ∈ W α,β , we obtain the following estimates: Taking notice of the Sobolev embedding constants in one domensional case with Ω = (0, 1) and J = (0, T ) (cf. [15, p. 8]), we have Here, in order to get the last inequality, we have used the fact that in case of the function with 0 only at one endpoint, i.e., t = 0, we can also apply the usual embedding theorem by considering a symmetric extension of the function on (0, T ) to (0, 2T ). Therefore, from w ∈ W α,β , we have On the other hand, from [10, Lemma 2], the V 1 L 2 norm of F(w) is estimated by: Noting that the Poincaré constant can be taken as 1 π with Ω = (0, 1), we obtain, by (47) Therefore, from (48), the following inequality holds: By these inequalities, we have the following sufficient condition for (46): Thus, we obtain the verification condition for the existence of the solutions of (44).
Since it holds that w = u−u k h , by solving the above simultaneous algebraic inequalities in α and β, we have error bounds of the form u − u k h L 2 H 1 0 ≤ α. We now present the verification results for the solutions of (44), namely, α, β, and the residual norm, ε L 2 L 2 . In Fig. 5, we chose the function f so that (43) has the exact solution u(x, t) = 0.5t sin(π x), with ν = 0.1 and ν = 1.0, which correspond to Example 1.1 and Example 1.2, respectively, in the previous subsection. We show more results in Fig. 6, in which the function f is chosen so that (43) has the exact