A Posteriori Error Bounds for Two Point Boundary Value Problems: A Green's Function Approach

We present a computer assisted method for generating existence proofs and a posteriori error bounds for solutions to two point boundary value problems (BVPs). All truncation errors are accounted for and, if combined with interval arithmetic to bound the rounding errors, the computer generated results are mathematically rigorous. The method is formulated for $n$-dimensional systems and does not require any special form for the vector field of the differential equation. It utilizes a numerically generated approximation to the BVP fundamental solution and Green's function and thus can be applied to stable BVPs whose initial value problem is unstable. The utility of the method is demonstrated on a pair of singularly perturbed model BVPs and by using it to rigorously show the existence of a periodic orbit in the Lorenz system.

1. Introduction.We propose a new computer assisted method for rigorously proving invertibility and bounding the norm of the inverse of operators of the form 1   (1.1) Operators of this form correspond to linear two point boundary value problems (BVPs) and arise naturally as the Fréchet derivative of operators defining nonlinear BVPs.Bounds on the norm of the inverse are required in a posteriori proofs of existence and local uniqueness via the Newton Kantorovich theorem [10].
Computer assisted methods have been successfully applied to many problems in differential equations and dynamical systems, including shadowing results for initial value problems [25,3,8], chaos in the Rössler [32] and Lorenz [6] systems, the existence of the Lorenz attractor [27], and flow (in)stability [2,30].Generalizations to elliptic and parabolic PDEs are discussed in, for example, [18,17,20,23].As it pertains to two point boundary value problems, computer assisted proofs are typically based on a fixed point theorem [31,22,19] or various forms of the Newton Kantorovich theorem, as used in [14,12,26].A method based on differential inequalities as presented in [7].
The methods in [14,7,22,26,19] focused on second order scalar equations; [14] worked with an equivalent integral formulation while [26,19] used a finite element based method.The results in [29,12] utilized the theory of initial value problem (IVP) fundamental solutions to treat n-dimensional two point BVPs, the former focusing on the periodic case.
Our method most closely resembles that of Ref. [12] and applies to n-dimensional systems without any special assumptions on the form of the vector field.The present method differs from [12] in the use of the BVP fundamental solution and Green's function as the primary tool, as opposed to the IVP fundamental solution.A general outline of the method is as follows.
(i) Generate approximations, Ψ j , to the BVP fundamental solution on a mesh.
(ii) Extend these to approximations Ψ(t) and Ψ −1 (t) on the whole interval.
(iv) Use G to define an approximate inverse, H, to the BVP operator F .(v) Using estimates involving H, prove existence of an exact inverse to F and derive a machine computable formula that bounds its norm.The use of the BVP fundamental solution makes the method presented here applicable in a range of situations for which IVP based methods are ill suited, namely when the BVP is stable but the corresponding IVP is not.This is often the case in singularly perturbed problems and so we demonstrate the effectiveness of the method on several such examples in section 9.
To obtain completely rigorous results, some means of bounding the floating point rounding error is needed, typically based on the theory of interval arithmetic.There are many references from which to learn more about validated numerics and interval arithmetic, for example [15,11,9,16,28].We present results so that if interval arithmetic is used for all arithmetic operations and interval valued extensions are available for the vector field of the differential equation and its derivatives then the resulting computer generated existence proofs and error bounds are mathematically rigorous (we perform tests using the MATLAB package INTLAB [24]).However, in light of the additional cost of interval arithmetic as compared to floating point, if complete rigor is not needed then the method can just as easily be used with traditional floating point operations to obtain approximate error bounds.
In section 2 we fix some notation and discuss a general strategy for a posteriori existence proofs via the Newton-Kantorovich theorem.For comparison, we outline the method of [12] in section 3. Background on the Green's function of a BVP is given in section 4. Development of our method begins in section 5. (2.1) y ′ (t) = f (t, y(t)), g(y(0), y(1)) = 0, where f : R n+1 → R is continuous, differentiable in y, D y f is continuous on R n+1 , and g : R 2n → R n is C 1 .In particular, we are guaranteed local existence and uniqueness of the IVP.

A Posteriori Existence
In order to apply functional analytic methods to this problem it is convenient to follow [12] and define the Banach spaces In this setting, the BVP Eq. (2.1) can be put in a form more amenable to analysis as follows.Define the map G : Solutions of the BVP are in one to one correspondence with zeros of G.As shown in [12], G is Fréchet C 1 with derivative at y 0 ∈ X 1 given by Dyf (s, y0(s))v(s)ds, (2.4) Having reformulated the problem in terms of finding zeros of a C 1 map between Banach spaces, a natural tool is the Newton-Kantorovich theorem.Invertibility of DG(y 0 ) along with a bound on the norm of its inverse are required to apply the theorem and these are typically the most difficult ingredients to obtain in practice.Once this has been done, and assuming f is C 2 , the remainder of the estimates are relatively straightforward and have been discussed in, for example, [12].For this reason we will now focus on the invertiblility of DG(y 0 ) i.e. on a posteriori methods for proving the existence and uniqueness to solutions of linear BVPs.
3. IVP Fundamental Solution Method.Motivated by the previous section, consider a linear BVP on [0, 1] as given by the bounded linear operator F : X 1 → X 2 defined by the formula Eq. (1.1), where A(t) is a continuous R n×n valued map.
In Theorem 1 of [12], the author shows how the solvability of F [v] = (r, w) is related to the matrix of fundamental solutions, defined as the solution to the initial value problem (IVP) We reproduce this result using our notation as the following theorem.Theorem 3.1.Let A(s) and r(t) be continuous matrix and vector valued resp.with r(0) = 0 and consider the two point BVP Let Y (t) be the corresponding matrix of fundamental solutions and define R = B 0 + B 1 Y (1).Then the BVP has a unique solution iff R is nonsingular and the solution is given by It was also shown in [12] how to apply this theorem to rigorously prove solvability of the BVP and derive guaranteed error bounds using numerically constructed approximations to Y (t) and Y −1 (t).Because of the use of Y , we will term this the IVP fundamental solution method or simply the IVP method.
A limitation of this method is that the stability properties of an IVP take center stage; the quantities |Y (t)| and |Y −1 (t)| feature critically in the error estimates derived in [12].Unfortunately there are situations in which a BVP of interest is stable whereas the corresponding IVP is not.In such situations the bounds obtained by the IVP fundamental solution method will be very pessimistic due to the large size of |Y (t)| and/or |Y −1 (t)|.
As discussed in [1], perhaps the simplest test problem that exhibits these features is the BVP The corresponding fundamental solution matrix and its inverse are and have exponentially increasing norms as a function of b.This is in spite of the fact that the value and derivative of the exact solution Even if one tries to avoid separately bounding the norms of Y and Y −1 and take advantage of the high degree of cancellation that can occur in Y (t)R −1 B 0 Y −1 (s) and Y (t)(R −1 B 0 − I)Y −1 (s), the above example is sufficient to show that this cancellation does not in general occur over the entire (t, s) domain; the (1, 1) component of the former quantity is given by which diverges at t = 0, s = b as b → ∞.This problem is in addition to the numerical difficulty of generating Y (t) and Y −1 (s) for unstable problems to high enough accuracy that one could even take advantage of cancellation.In summary, while the IVP fundamental solution method suffices for problems whose IVP is stable, if one is interested in BVPs whose corresponding IVP is unstable then a new method is needed.Our contribution is designed to address such situations by utilizing the BVP fundamental solution and Green's function in place of the IVP fundamental solution.Before we give the details of the method we first recall the general theory of Green's functions for linear BVPs.The treatment is based on [1] but in slightly greater generality.
4. Green's Function for Linear BVPs.Let Y be the matrix of fundamental solutions of the BVP and R = B 0 + B 1 Y (1) be invertible, so the BVP Eq. (3.2) has a unique solution, as guaranteed by 3.1.Define Φ(t) = Y (t)R −1 .As outlined in, for example, [1] one can verify that Φ satisfies Φ is called the matrix of fundamental solutions to the BVP.We note that any solution Φ of this BVP with Φ(0) invertible (and hence Φ(t) is invertible for all t) is related to the IVP fundamental solution by an invertible R in this way and hence proves the existence and uniqueness of solutions.The following theorem goes further by providing an analytic formula for the solution.
Theorem 4.1.If Φ exists then the unique solution to where the Green's Function is defined as Therefore Φ provides us with an inverse of the operator F from Eq. (1.1) given by (4.5) This result is essentially found in [1].The only complication beyond the presentation there is that r(t) must be taken to be an arbitrary continuous function that vanishes at t = 0, and not only one of the form for some continuous q, as would be the case when transforming from the differential to the integral form of the BVP.It is essentially then an application of integration by parts that leads to Eq. ( 4.3).One can directly verify the formula by composing Eq. (4.5) with Eq. (1.1) and vice versa.

5.
A More General Framework.We are now in position to begin discussing an alternative method for producing a posteriori existence proofs and error bounds based on the BVP fundamental solution.We refer to this method as the BVP fundamental solution method, or BVP method for short, to contrast it the IVP method.
The input to any a posteriori method consists of a certain (numerically generated) approximate solution.For us this will consist of function values on a mesh 0 = t 0 < ... < t N = 1 from which we will have to manufacture an approximate solution on the entire interval [0, 1] before we begin producing error bounds.While one could work in spaces of continuous functions by interpolating the values at the nodes, as done in the IVP method, we will find things more convenient (and the resulting error bounds more transparent) when this assumption is removed.In this section we outline the slightly generalized framework in which we will formulate the equations.
To this end, fix a mesh, t i , as above and let Y 1 be the space of piecewise continuous functions.More precisely, Y 1 is the direct sum On each factor we choose the norm to be a weighted ∞ norm where the weight, W , is some invertible matrix.In practice, we will take W to be diagonal and choose the diagonal elements to balance the magnitude of the errors in the different components of the solution.One could consider the generalization where the weight is allowed to vary between subintervals but we do not explore this here.Y 1 is a Banach space and the map that sends any y ∈ Y 1 to its value at any t ∈ {0, 1} ∪ N j=1 (t j−1 , t j ) or at t ± j (+ denots limit from above, − limit from below) is a well defined bounded linear map.Elements of Y 1 can be identified with elements of L ∞ and integration or differentiation (of piecewise C 1 ) elements of Y 1 will be defined using this identification.
Mirroring the continuous case, Eq. (2.2), we define a second space This is a Banach space under the norm (v, b) W = max{ v W , |b| W }. Finally, we point out that the X i 's are closed subspaces of the Y i 's.We will drop the subscript W from the above norms whenever the meaning is clear (i.e. when we are only working with the Y i 's).
Having specified the functional analytic framework, we can extend the definition of a linear BVP to this setting by defining F : Y 1 → Y 2 by the same formula as in Eq. (1.1).This extension is still a bounded linear map.
The formula Eq. ( 1.1) makes it clear that if r(t) is continuous on [0, 1] and though in applications it is probably wise to continue using the weighted norms in the Newton Kantorovich estimates (i.e. on X 1 and X 2 ), as we find a substantial improvement in the bounds from use of the weight.Eq. (5.4) shows that working in the extended framework will suit our purpose, namely proving invertibility and bounding the norm of the inverse of the operator F : X 1 → X 2 .
6. Proving Invertibility.We now suppose that we are given a piecewise C 1 map Ψ : [0, 1] → R n×n consisting of invertible matrices which we think of as an approximate fundamental BVP solution (but without any a priori knowledge that an exact Φ actually exists) and was obtained by non-rigorous numerical means.In practice, we will have some collection of nodes 0 = t 0 < ... < t N = 1 and approximations Ψ j at the nodes and will construct Ψ(t) as a piecewise polynomial on the intervals [t j−1 , t j ].However, for the present purposes, the only relevant fact is that Ψ is piecewise C 1 .
Our strategy for proving the existence of F −1 : Y 2 → Y 1 is to mimic the formula Eq. (4.5) for F −1 : X 2 → X 1 in terms of the fundamental solution but using the approximate fundamental solution instead.More precisely, if we replace Φ with Ψ in the formula for F −1 then we obtain a bounded linear map H : Y 2 → Y 1 (since Ψ is piecewise C 1 on the given mesh).
We will use H to find conditions under which we can prove the invertiblility of F .To do this we use the following well known result concerning perturbations of the identity, see for example [5].
Lemma 6.1.Let A : X → X be a bounded linear map on a Banach space.If I − A < 1 then A has a bounded inverse that satisfies (6.1) In fact, we will need a slight generalization of this result, given below.Lemma 6.2.Let A : X → Y , B : Y → X be bounded linear maps where X, Y are Banach spaces.If I − AB ≤ α for α < 1 where I is the identity on Y then A is surjective.If A is also injective then A, B are invertible, (6.2) Proof.By the above lemma, AB has a bounded inverse with In particular A is surjective.It is injective by assumption so by the open mapping theorem it has a bounded inverse.This implies B has a bounded inverse as well.
We have (6.4) 6.1.Applying the Fredholm Alternative.In finite dimensions and when X and Y have the same dimension then the hypothesis of injectivity of A in Lemma 6.2 is not needed but in infinite dimensions it is required in general.However, when we have the appropriate element of compactness, the situation in infinite dimensions mirrors the finite dimensional case.More specifically, we recall one consequence of the Fredholm alternative, see for example [4].Theorem 6.3.Let X be a Banach space and K : X → X be a compact operator.Then I − K is injective iff it is surjective.
If a Fredholm alternative-like condition holds for F then surjectivity will imply injectivity, and hence I − F H < 1 will imply the invertiblity of F .In addition, this will provide us with a bound on F −1 per Lemma 6.2.These are the two ingredients that are needed in order to utilize the Newton-Kantorovich theorem.With slight modifications of the above we could just as well work with I − HF but find this less appealing for reasons we point out later.
In order to more easily apply the Fredholm alternative, we reformulate the BVP operator F in terms of a map from a single Banach space to itself.Define F : ] and so F is surjective or injective if and only if F is.
The second component of K has finite rank and the first is an integral operator with bounded kernel over a bounded domain and whose image consists of continuous functions, hence the compactness of K follows from an application of the Arzela-Ascoli theorem.This implies that the Fredholm alternative applies to F and hence proving that I − F H ≤ α < 1 is sufficient to prove that F has a bounded inverse with (6.9) We will now discuss how to bound I − F H and H .
6.2.Bounding I −F H .In Appendix A we outline the derivation of a formula for I − F H. For each t ∈ {t i } N i=0 let j t = max{j : t j < t}.The two components of I − F H are given by G(s, z)A(z)r(z)dzds (6.10) where superscripts + and − denote the limits from above and below respectively.As we show in Appendix A, the second component vanishes exactly when the boundary conditions are satisfied and so bounding the first component is usually the more crucial estimate.For comparison, we give the analogous expression for I − HF in Appendix B and discussion why we prefer to formulate the BVP method in terms of I − F H.
The most straightforward way to use Eq.(6.10) to bound the norm of I − F H is to use the submiltiplicative property of induced matrix norms and bound each of the quantities in the Green's function individually.This leads to an estimate whose operation count scales linearly in the number of subintervals, N .If |Φ| and |Φ −1 | are not too large over [0, 1] then this is a good strategy.However, this is often not the case, even in situations where the BVP is stable.We illustrate this phenomenon with the same test problem introduced in Eq. (3.4).

Behavior of the Test Problem.
If we again consider the simple test problem, Eq. (3.4), for motivation, we find that the fundamental solution and its inverse are (6.12) While Φ is bounded, Φ −1 diverges as b → ∞.
The general notion of BVP stability only requires the norm of the Green's function to be small, not the fundamental solution and its inverse independently, see for example [1] for details.Therefore, in general we must abandon any estimate based on bounding Ψ and Ψ −1 separately and use one based on bounding the entire Green's function.The test problem Eq. (3.4) does have a bounded Green's function (6.13) We emphasize that G = O(1) does not arise from a large cancellation of the computed results.Rather, at each endpoint the boundary conditions annihilate the growing modes of Φ −1 .This suggests that we should not run into the same degree of numerical cancellation issues using the BVP method as opposed to the IVP method.
7. A Sharper Bound.In this section we show how to compute a sharper bound on I − F H by bounding the Green's function as a whole.Bounding the second component of I − F H is straightforward so we focus on the first.
To this point the only assumption we have made about Ψ is that it is piecewise C 1 .In practice, a non-rigorous numerical algorithm will provide us with approximate solution values at the mesh points, Ψ j , j = 0, ..., N , and we will construct Ψ| [tj−1,tj ] to be the polynomial of degree m that satisfies the ODE (7.1) up to order m − 1. Formulas for the coefficients and the ODE residual are given in Appendix C. Note that by our convention, the rightmost value Ψ N is ignored but the opposite choice could have been made.If the boundary conditions are very sensitive to perturbations then it may be wise to use a Taylor series centered at 1 over the rightmost interval, but we don't employ this technique here.
We emphasize that the resulting Ψ(t) will be piecewise smooth but not continuous.This is the reason we generalized the notion of solution in section 5 to allow for discontinuities.Though we could have modified Ψ to ensure that it interpolated continuously between the given data points, we feel that the error bounds one obtains from the method we have employed are clearer.One can see in Eq. (6.10) the separation into error due to Taylor series truncation on each interval [t j−1 , t j ], which is primarily controlled by the chosen polynomial degree m but also by the size of the mesh (especially if the exact solution is not analytic on the subinterval), and error due to boundary conditions and jumps, which are heavily influenced by the quality of the initial approximation.
We emphasize that in the absence of jumps (i.e. if we replace Ψ j with Ψ(t − j )) this method reduces to a piecewise Taylor model IVP method, see for example [21].While well suited for initial value problems, such methods perform extremely poorly for BVPs with an unstable IVP such as the examples in section 9. Therefore, having accurate Ψ j−1 's as inputs is crucial to the success of the method; the Taylor series simply ensure sufficiently small error over the interior of each subinterval of the mesh.
Using Appendix C we can write AΨ(t) − Ψ ′ (t) = R j (t)Ψ j−1 (t − t j−1 ) m on each [t j−1 , t j ] and break the integrals in Eq. (6.10) into sums of integrals over the subin-tervals of the mesh to obtain 7.1.Bounding the Inverse Fundamental Solution.The bound Eq. (7.2) is not yet in a machine computable form as it involves the exact inverse of Ψ, something we will not have access to in general.What we can do is compute approximate inverses Ψj of the Ψ j 's and define an approximate inverse Ψ(t) to be the polynomial on each [t j−1 , t j ] obtained by ensuring that I − Ψ Ψ vanishes up to some prescribed order in h j (apart from the additional floating point error in Ψj ).See Appendix D for details.
Terms of the form |D 1 Ψ −1 D 2 | for some D i can then be bounded using Lemma 6.1, where we require |I − Ψ(t) Ψ(t)| < 1 for all t so that the indicated inverses are guaranteed to exist.When being automated on a computer, if this bound cannot be proven by the program then a failure message should be returned.

A Machine Computable
Bound.With this we obtain a form of the bound that is machine computable where the subscript j on • j indicates that the supremum is taken only over [t j−1 , t j ].In bounding R j on [t j−1 , t j ] we assume that we have an interval matrix valued extension for A(t) and its derivatives so that we can rigorously bound the error term in its Taylor series expansion, or have some alternative means of computing an upper bound analytically.To bound the polynomial terms we use the simple method of summing the bounds on the norms of each monomial.
The leading order term in N of the complexity, which comes from the loops involving the Green's function like terms such as where k is the degree of Ψ and n is the dimension of the ODE system.Of course, as discussed above we can drop the algorithm from N 2 to N by using the submultiplicative property to break up the approximate Green's function terms, but this result in a huge loss in accuracy in cases where the fundamental solution is ill conditioned.Intelligently choosing when such an approximation is acceptable is one potential avenue for improving the performance of this method.
A factor of k speedup can be achieved if we use Appendix D to write Ψ(t) = Ψj−1 (I + E j (t)) and bound With the use of an appropriately chosen weight W in the norm, which can be adaptively chosen by algorithm as discussed below, the overestimation due to this step is often minor compared to the speedup achieved and so we will employ this when testing the method.Finally, to use Eq. ( 6.9) to bound the norm of F −1 we must have a machine computable bound on H , which we now present.
The leading order complexity in N of this formula is O(mkn 3 N 2 ) where m is the degree of Ψ.Similar to I − F H , a factor of mk speedup can be achieved by the estimate Eq. (7.6) together with the analogous estimate for Ψ using Eq.(C.3).Again, with the proper choice of weight this is a relatively insignificant overestimation and so we will use it when testing the method.With these approximations, the overall cost of the method to leading order in N is O(n 3 N 2 ).In particular, the degree of Ψ can usually be set as high as is necessary for the Taylor series errors over the the subintervals to be negligible, making the boundary conditions and jumps the main sources of error.We also note that it is trivial to parallelize the computation of these bounds.Now that we have machine computable bounds on I −F H and H the method (coupled with the Newton Kantorovich theorem) can be applied to rigorously prove existence and derive error bounds to nonlinear two point BVPs.In the linear case most of the complications in the statement of the Newton Kantorovich theorem vanish and the derivation of error bounds is much simpler.We discuss how this in done in the following section.
First we summarize these results with a theorem.Theorem 7.1.Let Y 1 and Y 2 be defined by Eq. (5.1) and Eq.(5.3) respectively and consider the two point BVP F where where G is the approximate Green's function as computed by Eq. (7.4) and Eq.(7.5), then both F −1 : Y 2 → Y 1 and F −1 : X 2 → X 1 exist (X i are defined in Eq. (2.2)), are bounded, and A bound on H can be computed using Eq.(7.7).
8. Validated Solution of Inhomogeneous Linear BVPs.In this section, we show how to apply the BVP method to produce a posteriori error bounds for solutions of an inhomogenous linear BVP.Suppose we have proven I − F H ≤ α < 1 so that F is invertible.Let v be the exact solution of F [v] = (r, w) and suppose we are given an approximate solution ṽ.Then by Eq. (6.9), Note that we don't require ṽ to be continuous, even though v is whenever the inhomogeneity r is.We only need it to be piecewise continuous on the mesh 0 = t 0 < ... < t N = 1.We suppose that A(t) and r(t) are C m and so on each subinterval [t j−1 , t j ] we have Taylor series with remainders (8.2) With this we bound the components of the residual in Eq. (8.1) as where as above, • j denotes the sup norm over [t j−1 , t j ].
Similar to the approximate fundamental solution, we assume we are given values ṽj at the nodes.We construct ṽ(t) on [t j−1 , t j ] to be the polynomial which causes the degree 1 through m terms in Eq. (8.4) to vanish.We don't have the freedom to choose the zero order term; it is the analog of the jump terms from Eq. (6.10).
If we analytically evaluate the integrals in Eq. (8.4) then the bound on the first component of the residual becomes where The second component of the residual, Eq. (8.6), is much simpler.One just needs to remember to evaluate ṽ(1) using the polynomial Eq. (8.7) and not the mesh value ṽN .The bound is now in a form that can easily programmed for automated rigorous evaluation using interval arithmetic or estimation in traditional floating point.We also note that the cost is just O(N ) as opposed to the O(N 2 ) cost of bounding F −1 .9. Numerical Tests.In Eq. (7.4) and Eq.(7.7) we gave bounds on I−F H and H that involve finitely representable objects and finite operations with which we can perform validated computations, namely the piecewise polynomials Ψ(t) and Ψ(t) and the interval valued extension of A(t) and its derivatives.In practice, we implemented this using the MATLAB interval analysis package INTLAB [24] by creating simple routines that perform algebraic operations with polynomials whose coefficients are interval matrices.With such functionality, the code that produces the desired bound on I − F H is a simple translation of Eq. (7.4) and Eq.(7.5)In this section, we assess the applicability and accuracy of the BVP method using several singularly perturbed test problems.As the extension to nonlinear problems of a similar method has been discussed elsewhere [12], we focus on the linear case.We take the test problems from [13].
If we let err i j be the ith component of the jth jump error from Eq. (8.4) then weight we used in these tests has diagonal elements given by (9.1) This choice prevents overestimation of the error in certain components when they have widely different scales and makes a large difference in the quality of the bounds, up to several orders of magnitude; it was often the difference between the existence test I − F H < 1 passing or failing.Since the jump errors can be computed in O(N ) operations prior to the main O(N 2 ) computation, this is an effective and practical way to improve the estimates.We also note that to obtain usable bounds, even the exponentially small parts of the approximate fundamental solution need to be resolved with the same relative accuracy as the larger portions, otherwise the loss of accuracy will prevent the growing modes from being annihilated when forming the approximate Green's function.In the tests, we set the absolute tolerance of the non-rigorous solver near the underflow limit in order to achieve this.
We performed all tests using a uniform mesh.There is room for improvement by intelligently selecting the mesh but we did not employ such methods here.The mesh size was chosen by adding points until the error approximately stabilized.We used k = m = 15 with the aim of making the Taylor series and approximate inverse errors negligible, allowing us to focus on how the error in the initial approximation supplied to the BVP method translates into the bound produced by the method.9.1.Example 1: Turning Point.Our first test problem is a singularly perturbed Airy equation which is a model of the quantum mechanical wave function near a turning point in the potential, here at t = 1/2.The exact solution is a linear combination of Airy functions and is shown in figure 1, where ǫ = 10 −6 .The solution exhibits dense oscillations for t < 1/2 and a boundary layer near t = 1.We ran the BVP method for the cases ǫ = 10 −4 , 10 −5 , 10 −6 .The test results are shown in table 1.The fourth column gives the error in solution values on the mesh that are provided to the BVP method, the fifth column is the bound on ||I − F H|| computed by the method, and the final column is the computed bound on the absolute value of the error in v.
The a posteriori BVP method does overestimate the error by several orders of magnitude, but the error bounds are still practical.The IVP fundamental solution method would not be usable in this situation, as the fundamental solution matrix has a norm |Y (1)| ≈ 10 12 even for ǫ = 10 −4 .We also note that the use of an adaptive weighted norm is crucial for the success of the method and greatly improves the error bounds.Without the weight, the existence test I − F H < 1 fails even at ǫ = 10 −5 .
The BVP method fails on this problem for ǫ = 10 −7 , as Ψ becomes so ill conditioned that its decaying components cannot be resolved to sufficient relative accuracy due to underflow.This point will reappear in the subsequent examples as well; the exponentially decaying modes of Φ (which tend to be the least 'interesting' features) cause the most trouble for the BVP method.9.2.Example 2: Viscous Shock.As a second example we consider a type of steady advection diffusion equation  which has the exact solution where erf is the error function.The solution has a transition layer at t = 1/2, shown in figure 2 (left pane).Again we use a uniform mesh in the tests for illustration.It is important to recall that the traditional choice of mesh, dense near the transition and sparse elsewhere, is not appropriate here as even the decaying modes of the BVP fundamental solution need to be resolved to the same relative error as the growing modes .This is illustrated by the results in table 1.We obtain good error bounds for ǫ = 10 −3 but the method fails for ǫ = 10 −4 due to overflow in computing the approximate inverses Ψj .
We also tested the cusp problem from [13] with similar results to this example, and so we don't discuss it further.Both are examples are very stiff problems with which the BVP method has difficulty for small ǫ due to under/overflow.9.3.Example 3: Potential Well.As a final test, we consider a quantum mechanical potential well problem (9.5) where we take ω = 1/4.The solution is oscillatory on [0, 1/4] ∪ [3/4, 1] and exponentially decaying as one moves towards the center of [1/4, 3/4] as seen in figure 2 (right pane) for ǫ = 10 −6 .As with the other test problems, the method fails once the decaying modes start to underflow, which occurs for ǫ = 10 −7 .
For this example we don't have an analytical formula for the exact solution, so the input error reported for in table 1 for this problem is obtained by comparing approximate solution values after refining the mesh.In the first two examples, the supremum norm of the exact solution was relatively consistent between the different parameter values, but that is not true for the current example, where v ≈ 300 for ǫ = 10 −5 and is single digits ǫ = 10 −6 .This explains why the absolute error is so similar for the cases tested despite the difference in ǫ.

Conclusion.
We have presented an a posteriori method for proving existence to and deriving rigorous error bounds for two point boundary value problems on a finite interval, summarized in Theorem 7.1.The method applies to general ndimensional systems without any assumption on the form of the vector field.It uses a (non-rigorous) approximation to the Green's function of the BVP to generate the bounds and can be evaluated using interval arithmetic for mathematically rigorous results or in traditional floating point for faster generation of approximate bounds.Because the method is based on the Green's function of a BVP rather than the fundamental solution to an IVP the method is applicable to cases where the BVP is stable but the corresponding IVP is unstable.We also utilize an adaptively chosen weight matrix in the norms of the estimates.This improves the quality of the final L ∞ error bounds by several orders of magnitude in many of the the tests.
In section 9 we tested the BVP method on a several singularly perturbed model problems that exhibit such problematic features as dense oscillation, boundary layers, and sharp transition regions.The method is successful in proving existence and producing reasonable and usable error bounds, even when the relevant perturbation parameter becomes moderately small.The large improvement we obtain in the domain of applicability and quality of the bounds does come at an increased cost as compared to other methods, with complexity scaling with the square of the number of subintervals in the mesh, N .However, in less numerically delicate situations where other methods are applicable the complexity can be reduced to O(N ) by weakening the bounds.
In the cases tested, the main limiting factor on the success of the method is the extreme ill conditioning of the approximate BVP fundamental solution Ψ as the parameter that controls the singular perturbation is made smaller.Although Ψ itself has moderate sized norm in the tests, it is close to being singular.The BVP method requires computation of the fundamental solution and its inverse to sufficient relative accuracy; poorly resolved decaying modes will spoil the method and the method always fails once the underflow limit is reached.This commonly occurs for very stiff problems and so, while the BVP method works for moderately stiff problems as we have shown, it fails in the very stiff limit.This issue can be addressed in a brute force manner by developing a floating point arithmetic package that allows for substantially larger and smaller exponent values but such an approach is likely to be quite inefficient.The use of asymptotic formulas in the exponentially decaying regions is perhaps a more promising avenue Proofs via Newton-Kantorovich.Consider a general non-linear two point BVP on [0, 1] where we use the ∞-norm on R n , denoted by | • |, the sup-norm y ≡ sup [0,1] |y(t)| on the function spaces, and on X 2 we use the norm (y, b) = max{ y , |b|}.

Table 1
Test results for the BVP method on the example problems.