A new formulation using the Schur complement for the numerical existence proof of solutions to elliptic problems: without direct estimation for an inverse of the linearized operator

Infinite-dimensional Newton methods can be effectively used to derive numerical proofs of the existence of solutions to partial differential equations (PDEs). In computer-assisted proofs of PDEs, the original problem is transformed into the infinite-dimensional Newton-type fixed point equation w=-L-1F(u^)+L-1G(w)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w = - {\mathcal {L}}^{-1} {\mathcal {F}}(\hat{u}) + {\mathcal {L}}^{-1} {\mathcal {G}}(w)$$\end{document}, where L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}$$\end{document} is a linearized operator, F(u^)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {F}}(\hat{u})$$\end{document} is a residual, and G(w)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {G}}(w)$$\end{document} is a nonlinear term. Therefore, the estimations of ‖L-1F(u^)‖\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert {\mathcal {L}}^{-1} {\mathcal {F}}(\hat{u}) \Vert $$\end{document} and ‖L-1G(w)‖\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert {\mathcal {L}}^{-1}{\mathcal {G}}(w) \Vert $$\end{document} play major roles in the verification procedures . In this paper, using a similar concept to block Gaussian elimination and its corresponding ‘Schur complement’ for matrix problems, we represent the inverse operator L-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}^{-1}$$\end{document} as an infinite-dimensional operator matrix that can be decomposed into two parts: finite-dimensional and infinite-dimensional. This operator matrix yields a new effective realization of the infinite-dimensional Newton method, which enables a more efficient verification procedure compared with existing Nakao’s methods for the solution of elliptic PDEs. We present some numerical examples that confirm the usefulness of the proposed method. Related results obtained from the representation of the operator matrix as L-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}^{-1}$$\end{document} are presented in the “Appendix”.


Introduction
In this paper, we study a new approach to proving the existence of solutions to elliptic problems. The proposed approach offers an improvement over existing Nakao's methods that use a finite-dimensional projection, that is, FN and IN methods in [11] and [9,Part I]. Particularly, an important aspect of our approach is the formulation using the Schur complement without direct estimates of an inverse of the linearized operator. Our approach inherits the advantages of both Nakao's FN and IN methods using a finite-dimensional projection, which also indicates that the disadvantages of both methods are resolved.
In the FN method, with an appropriate setting of the finite-dimensional subspace V h ⊂ H 1 0 (Ω), we first consider the Ritz projection R h : H 1 0 (Ω) → V h defined as for u ∈ H 1 0 (Ω), where I denotes the identity operator on H 1 0 (Ω). As a property of R h , we assume that there exists a positive constant C(h) which can be numerically estimated satisfying with the property that C(h) → 0 as the parameter h → 0. Using this projection, the problem is decomposed into two parts: finite-dimensional and infinite-dimensional.
Let ψ 1 ,. . . ,ψ N be a basis of V h , and let V ⊥ : be an orthogonal complement of V h . For a given approximate solutionû ∈ V h , setting w := u −û, w h := R h w, and w ⊥ := (I − R h )w, the FN method uses the following fixed point formulation: where A : is equal to minus the weak Laplace operator. In [4], the candidate set of solutions was set to and the fixed point theorem was applied to (3) to verify a solution in the set W = W h + W ⊥ . To confirm the verification condition in the fixed point theorem, the verified computation of solutions for a linear system of equations and the constructive a priori error estimates (2) for the Ritz projection play an essential role. This was the first approach to numerical verification, but is based on a contraction assumption which is more restrictive compared with the FN method. The FN method applies Newton's method to the finite-dimensional part of (3) as follows: Let f [û] be the Fréchet derivative atû of the nonlinear term f (u), and let L : be a linear operator defined as Furthermore, the finite-dimensional operator T : V h → V h is defined as and T is assumed to be invertible. Then, we can rewrite the finite-dimensional part of (3) as This FN method was developed in [5], and it has been confirmed that Newton's method is more effective for the verification of w h than the method in [4]. Additionally, because the computation of w h by the iterative use of (8) is sufficient to simply solve the linear system of equations, the cost is not large compared with the matrix norm estimation of T −1 . However, when f (u) includes the first-order derivative ∇u, because of the property of interval arithmetic, it sometimes causes unexpected inefficiencies in the estimation of the right-hand side of (8); that is, the inefficiency appears as a result of the effect of the corresponding estimation of ∇w ⊥ and leads to the failure of verification because of an explosive expansion of intervals in the iteration (e.g., see [9,10]). This implies that the effect of Newton's method on w ⊥ has not been achieved. To overcome this difficulty, in [9,10], the FN-Norm method was introduced. This ensures a more effective verification procedure by setting the candidate set W h of the finite-dimensional part to By contrast, the IN method assumes that the linearized operator L is invertible and considers the following fixed point equation: where F(û) := Aû − f (û) (10) and As (9) is a Newton-type formulation in the infinite-dimensional sense, the linear term with respect to w, which is a shortcoming of the FN method, is no longer present on the right-hand side of the equation. Therefore, a Newton-Kantorovich-like theorem can be derived using the fixed point equation (9). However, as the inverse operator L −1 cannot be calculated directly, it is necessary to show the invertibility of L and estimate the operator norm L −1 . Thus, the evaluation of ||L −1 || is the major task in the verification procedures of the IN method, and has been studied by many researchers since 1991. For example, [3,8,12,13,[19][20][21][22] and [9, Part I] used the Ritz projection R h and the projection error bound C(h) by limiting Ω to bounded domains. Other methods exist that avoid the Ritz projection R h and its error bound C(h), and hence can be applied even for unbounded domains Ω, where the projection error bound C(h) is not accessible or does not even exist (e.g., [14][15][16] and [9, Part II]). IN methods that do not need a projection error bound and methods that use projection error bounds have strengths and weaknesses, and a comparison of their qualities on a general level is impossible. Generally, the IN method defines the candidate set as Therefore, the IN method overestimates the finite-dimensional parts compared with the FN method. We now describe a new approach that incorporates the advantages of both Nakao's FN and IN methods using the Ritz projection R h in [11] and [9,Part I]. In this paper, we essentially consider the IN method based on (9), but we propose a verification method without estimating the norm L −1 . Through a computational procedure that avoids the direct evaluation of the operator norm L −1 , highly accurate and efficient verification can be expected. We decompose (9) into finite-dimensional and infinitedimensional parts; that is, However, we cannot directly calculate L −1 for the same reason as for the existing IN method. To overcome this difficulty, we introduce an operator matrix, whose existence will be proved later, satisfying Here, we apply the fixed point theorem to the following equation: with candidate sets (4) and (5). Note that the right-hand side of this fixed point equation no longer has linear terms in w h nor w ⊥ . Therefore, the proposed method overcomes the disadvantages of the FN method. Additionally, note that we can directly compute the finite-dimensional part (e.g., H 11 R h A −1 F(û)) by solving the linear system of equations, which is an advantage of the FN method. Therefore, the proposed method also overcomes the disadvantages of the IN method. Thus, our approach inherits the advantages of both Nakao's FN and IN methods using a finite-dimensional projection in [11] and [9, Part I], which also indicates that the disadvantages of both methods are resolved; that is, the disadvantage of the Newton method not working for the infinitedimensional part of the FN method and the disadvantage of the poor computational efficiency of the finite-dimensional part in the IN method are removed. For the actual implementation of the verification procedure, it is necessary to obtain a more detailed form of the operator matrix H . Thus, we consider a specific construction of this matrix below. The remainder of this paper is organized as follows: In Sect. 2, we describe how to construct the operator matrix H . In Sect. 3, we present the results of numerical experiments using the operator matrix H . We also describe why the proposed method offers an improvement over previous techniques based on some useful results in the "Appendix".

Constitution of the inverse operator matrix H
In this section, we present a detailed description of the actual construction of the operator matrix H satisfying (14) and (15). The basic idea comes from the concept of block Gaussian elimination and its corresponding 'Schur complement' for matrix problems.
We consider a solution φ that satisfies the linear equation for a given g ∈ H −1 (Ω). We denote We multiply A −1 from the left on both sides of (17), and decompose the result into the finite and infinite-dimensional parts using the Ritz projection R h as follows: where I V ⊥ is the identity operator on V ⊥ . Furthermore, transforming the above equation using the operator matrix yields respectively. Additionally, we define the 2 × 2 block operator matrix D as Moreover, if the operators L and D are invertible, then we have and H is equal to D −1 from (15). We first present a sufficient condition for the invertibility of the operator D, which also provides a detailed expression for D −1 . (7) is assumed to be invertible. Let S : V ⊥ → V ⊥ be a linear operator defined as which is the so-called Schur complement corresponding to the operator T of the block operator matrix D. If S is bijective, then the operator D defined by (20) is also bijective and we have where I V h is the identity operator on V h .
Proof . We first apply block Gaussian elimination to the operator matrix D. We multiply M 1 from the left of D as follows: This procedure is so-called block forward elimination. Next, we multiply M 2 from the left of M 1 D as From the assumption that T and S are bijective, we have Note that the multiplication of M 3 M 2 corresponds to so-called back substitution.
Next, to show that M 3 M 2 M 1 is D −1 , we multiply M 3 M 2 M 1 from the right of D as follows: We obtain the inverse operator of D using Lemma 1. However, the operator matrix H satisfying (15) cannot be introduced without the invertibility of L. Therefore, the following theorem provides a sufficient condition for the invertibility of L.
Theorem 1 Under the same notation and assumptions as in Lemma 1, the linearized operator L defined by (6) is bijective, and the operator matrix H satisfying (15) is given by Proof We show that the linearized operator L is bijective. From the assumptions and Lemma 1, we have an inverse operator of D. Multiplying D −1 from the left of (18), we obtain To prove that L is injective, we show that φ = 0 is the only solution of the equation Lφ = 0. This can be readily seen from the fact that, using the bijectivity of operators D and A, the solution of (18) with g = 0 implies that φ = 0.
Finally, we prove that the linearized operator L is surjective. For this, it is sufficient to show that there exists a solution φ ∈ H 1 0 (Ω) satisfying Lφ = g for any g ∈ H −1 (Ω). Now, for g ∈ H −1 (Ω), define (φ h , φ ⊥ ) T as the left-hand side of (22) and set φ := φ h + φ ⊥ . Then, by the invertibility of the operator matrix H , the function φ satisfies (18), and therefore it also implies (17), which proves the desired surjectivity. Thus, the operator L is bijective.
Theorem 1 provides a specific expression of the operator matrix H satisfying (14). Moreover, Theorem 1 provides a new expression for the solution φ of the linear noncoercive elliptic PDE Lφ = g. Therefore, for example, the exact expression of the Ritz projection error for the solution of the linear noncoercive elliptic PDE Lφ = g is also derived from Theorem 1 (see Appendix C for details).
At the end of this section, we describe how to verify that T and S in the assumptions of Theorem 1 are invertible operators. Let G ∈ R N ×N be a real matrix defined as Then, the matrix G is invertible if and only if T is the invertible operator (e.g., [9, Lemma 2.1, p.44]). Therefore, by confirming the invertibility of the matrix G using numerical computation with guaranteed accuracy, we can easily check the invertibility of T .
To confirm the invertibility of the operator S, the following operator norm is used. For two Banach spaces X and Y , the set of bounded linear operators from X to Y is denoted by L(X , Y ) with the operator norm B L(X ,Y ) := sup{ Bu Y / u X : u ∈ X \ {0}} for B ∈ L(X , Y ). When X = Y , we simply use L(X ). Then, we can confirm the invertibility of S using the following well-known theorem. (21).

then S is bijective and its norm satisfies
Proof The Ritz projection R h is a continuous projection; therefore, V h and V ⊥ are Hilbert spaces with the inner product (·, ·) H 1 0 , respectively. The definition S : | V ⊥ and the assumption κ < 1 yield the bijectively of S as a result of the well-known theory of the Neumann series.
To compute the constant κ, we use the inequality (2) and the Aubin-Nitsche trick There are several ways to compute the constant κ. For example, in the case of f [û] ∈ L(L 2 (Ω)), we define the constants C f ,L 2 and τ L 2 as Then, we can estimate κ as Note that τ L 2 can be estimated as the matrix norm L where τ L 2 ,H 1 0 can be also estimated using verified computing (e.g., see [9, (3.63), p.93]). Then, we can estimate κ as The estimation (24) is the same as [9, (3.64), p.93]. It is also possible to derive an evaluation similar to [9, (3.46), p.89] from Lemma 2. Moreover, since C(h) has the property that C(h) → 0 as h → 0, it is expected to satisfy the sufficient condition κ < 1.

Verification procedure
In this subsection, we describe a verification procedure to realize computer-assisted proofs using the fixed point formulation (16) with the operator matrix H . The proposed verification method is a combination of the FN and IN methods developed above. Letû ∈ V h ⊂ H 1 0 (Ω) be a solution that satisfies the following equation: For example,û may be obtained by numerical computations with guaranteed accuracy for finite-dimensional nonlinear equations, such as the Krawczyk method. Note that (16) becomes zero. Next, we verify that the matrix G defined as and the operator S are invertible. As Theorem 1 implies that the linearized operator L defined as (6) is invertible, we can transform problem (1) into the fixed point equation (16) using the operator matrix H .
Schauder or Banach's fixed point theorem may be applied to the equation (16) with the candidate sets (4) and (5), as in [4,5]. Here, the fixed point theorem may be selected as follows. If the nonlinear term f (u) is similar to that in [15] (e.g., f (u) ∈ L 2 (Ω) ∀u ∈ H 1 0 (Ω)), then because L −1 G is a compact operator, Schauder's fixed point theorem can be used. If this is not the case, or if we need to prove the local uniqueness of the solution, it is preferable to use Banach's fixed point theorem. A survey of the FN method [6] makes it easy to select the appropriate method.

Example
In this subsection, we present an example in which our method is used to verify a solution of the elliptic boundary value problem with Ω = (0, 1) 2 . This is Emden's equation, which is a good test problem for comparing the proposed method with other approaches. Additionally, because the results for ρ in (12) given by many existing IN methods [8,15,16,18] are almost the same, it is not necessary to compare numerous existing IN approaches. Emden's equation has also been discussed in terms of the FN method [23].
All computations were implemented on a computer with 2.20 GHz Intel Xeon E7-4830 v2 CPU × 4, 2 TB RAM, and CentOS 7.2 using C++11 with GCC version 4.8.5. All rounding errors were strictly estimated using kv library [1]. This guarantees the mathematical correctness of all the results.
We constructed approximate solutions for (26) using a Legendre polynomial basis. Specifically, we defined the set {ψ 1 , ψ 2 , . . . ψ N } of Legendre polynomials as and defined the finite-dimensional subspace as a tensor product See [2] for information on how to compute C(h) in (2) corresponding to the Legendre polynomial basis. Then,û ∈ V N h satisfying the finite-dimensional nonlinear problem (25) can be written asû whereû i, j are real numbers. Note that, in a real computation,û i, j are obtained as intervals that include the exact solution of (25) using the Krawczyk method. For example, when N = 10, a solutionû satisfying (25) can be obtained (see Table 1). Here, 1.23 789 456 denotes the interval [1.23456, 1.23789]. As the solution has symmetry, note thatû i, j is zero when i and j are even. Under this setting, exact solutionsû ∈ V 10 h satisfying (25) were computed numerically. Their graphs are displayed in Fig. 1. The computational results of the constants C(h) and κ using (23) were 0.037268163011998514 and 0.31088290279527165, respectively. Therefore, because κ < 1, the operator S was bijective. Thus, taking the candidate set as the method proposed in this paper succeeded in the numerical verification of problem (26). The verified W h results for the finite-dimensional parts are presented in Tables  1 and 2, and the α results for the infinite-dimensional parts are given in Table 2. Furthermore, we can prove that the exact solution u * of (26) exists inû + W h + W ⊥ , and we can estimate FN-Int (e.g., [5,6,11]) was also applied using similar candidate sets (27) and (28) and performing a detailed evaluation of the finite-dimensional and infinite-dimensional parts. However, as FN-Int only applies Newton's method to the finite-dimensional  parts, the verification failed for N = 10, which implies that N = 10 was too small to achieve a successful verification.
As the error bound of the form u * −û H 1 0 ≤ ρ was obtained in the course of a successful verification by the IN method for N = 10, we compared the proposed method with the IN method [8,18] from this viewpoint in Table 2. It is clear from the table that the value of ρ given by the proposed method was smaller than that produced by the existing IN method [8,18]. Additionally, as described in Sect. 3.1, because we partially incorporated the detailed calculations of the IN method into the proposed method, enhancing the IN method could lead to improvements to the results of the proposed method. The reason why the proposed method produced smaller values of ρ than [8,18] is discussed in Appendix B.
From Table 1, which presents results for N = 10, the error interval of the finitedimensional part appears to be rather large. This is because N was too small. In fact, very sharp results were obtained in the case of N = 40 (see Tables 3 and 4). Additionally, we obtained C(h) = 0.01150109366291357 and κ = 0.029612643887446451 using (23). Because C(h) was smaller than the case of N = 10, the constant κ was also small . Here, ±1.

Conclusion
In this paper, we described a new formulation (16) for the numerical proof of the existence of solutions for elliptic problems. The proposed approach has advantages over both Nakao's FN and IN methods using the Ritz projection R h in [11] and [9, Part I]. In particular, we derived a specific formula for the operator matrix H , which is needed to compute (16), in Theorem 1. As a result, while using the infinite-dimensional Newton method, the error evaluation for each coefficient interval of the finite basis is also enclosed, as demonstrated by the results in Tables 1 and 3. This is considered an advantage of the FN method. Furthermore, the proposed method produced better results than the IN method [8,18], even for the norm estimation, as demonstrated in Tables 2 and 4. Moreover, it is possible to derive, from the operator matrix D in Theorem 1, the results in previous papers [8,12,21] for the norm evaluation L −1 L(H −1 ,H 1 0 ) of the inverse operator L −1 (see Appendix B). However, there has been no previous discussion (e.g., [9]) of the norm evaluation of the operator L −1 : H −1 (Ω) → H 1 0 (Ω) using the operator matrix D in Theorem 2.  [8,9,[12][13][14][15]19,21]). In the following, we show that it is also possible to obtain the upper bound of the inverse operator norm L −1 L(H −1 ,H 1 0 ) using Theorem 1.

Corollary 1 (of Theorem 1) Under the same assumptions as in Theorem 1, L is invertible, and it follows that
where · E denotes a matrix norm induced by the Euclidean vector norm | · | E .
Proof We define the norm of the direct product space Then, from Theorem 1, the solution φ of the linear equation (17) can be evaluated as Therefore, we have Moreover, using the structure of the operator matrix H satisfying (14), we have

Remark 1
The estimation of the inverse operator norm in Corollary 1 is closely related to the method used in previous studies [8,9,12,21]. However, it is noted that Theorem 1 can be applied directly without using Corollary 1. For example, for the first term on the right-hand side of (9), the existing IN method evaluates .

Appendix C An exact expression for the Ritz projection error in the solution of the noncoercive equation (17)
In this appendix, we derive an expression for the Ritz projection error φ − R h φ in the exact solution φ of the linear noncoercive elliptic PDE (17). Using this result, we can, for example, determine a constant C that satisfies the inequality φ − R h φ H 1 0 ≤ C Lφ L 2 , φ ∈ {φ ∈ H 1 0 |Δφ ∈ L 2 (Ω)}, even though the elliptic operator L is noncoercive (cf. [7]).

Corollary 3 (of Theorem 1)
Under the same assumptions as in Theorem 1, the Ritz projection error of the solution φ ∈ H 1 0 (Ω) of the linear noncoercive elliptic PDE (17) is Proof The proof is obtained immediately from φ ⊥ in the expression of the proof for Theorem 1.
Using Corollary 3 as g ∈ L 2 (Ω), it is easy to derive the constant C satisfying φ − R h φ H 1 0 ≤ C Lφ L 2 in [7]. Similarly, we can derive the following proposition from Theorem 2.

Corollary 4 (of Theorem 2)
Under the same assumptions as in Theorem 2, the Ritz projection error in the solution φ ∈ H 1 0 (Ω) of the linear noncoercive elliptic PDE (17) is represented as follows: