Mesh selection strategies of the code TOM for Boundary Value Problems

This paper presents new hybrid mesh selection strategies for boundary value problems implemented in the code TOM. Originally the code was proposed for the numerical solution of stiff or singularly perturbed problems. The code has been now improved with the introduction of three classes of mesh selection strategies, that can be used for different categories of problems. Numerical experiments show that the mesh selection and, in the nonlinear case, the strategy for solving the nonlinear equations are determinant for the good behaviour of a general purpose code. The possibility to choose the mesh selection should be considered for all general purposes codes to make them suitable for wider classes of problems.


Introduction
The development of a general purpose code for the solution of two point Boundary Value Problems (BVPs) described by g(y(a), y(b)) = 0 ( 1 ) with f (x, y) ∈ IR m and g continuous functions, requires to pay a special attention to: the error estimation, the mesh selection strategy and the method for solving the nonlinear equations. Codes based on the same numerical scheme show completely different behaviour depending on how these aspects are exploited. An interesting work related to the choice of the error estimation used in the mesh selection strategy has been made in [1] where the authors show that properly changing the mesh selection, changes the behavior of the code and its efficiency. In [1] three strategies for estimating the global error are compared: the one based on Richardson Extrapolation (RE), the one based on the use of a higher order formula (HO), the one based on the use of deferred correction (DC). Additionally, the authors consider also a global error bound based on the estimation of the conditioning constant and the defect. The mesh selection and the method for solving the nonlinear problems remain essentially the same, but the final results differ. The conclusion is that the HO and DC strategies give the best results compared to RE. The authors suggest that the best is to have both the strategies available and to leave to the user the choice. The last release of the code BVP_SOLVER includes four different error estimation routines. A comparison of the efficiency of the code changing the mesh selection has been made for the code TOM in [2] where an hybrid mesh selection strategy based on conditioning has been analyzed and in [3,4] where the mesh selection introduced in [2] has been adapted to the Fortran codes TWPBVP and TWPBVPL. The last release of these codes allow to choose two different mesh selections: the hybrid one based on conditioning suggested for stiff or singularly perturbed problems and the one based on the estimation of the global error using a deferred correction strategy. In [5] the Fortran codes have been translated in Matlab and the mesh selection based on conditioning has been added to the code ACDC, that uses an automatic continuation strategy. We observe that the error estimation and the strategy for solving the nonlinear equations have not been changed in all the codes. The mesh selection strategy needs a special attention when the problem is considered stiff or singularly perturbed. For these classes of problems many studies have been carried out, often associated to the definition of new numerical methods. In this case the solution presents regions of rapid variation that can be located at the boundary (boundary layers) or in any points internal to (a, b) (turning points). In particular we recall the Shishkin mesh and the Bakhvalov mesh strategies, that have special properties of uniform convergence [6,7]. These strategies are effective when the location and the width of the layer are known, even if recently, adaptive meshes for weakly coupled system of singularly perturbed delay differential equations have been proposed in [8]. Robust hybrid methods of higher order for singularly perturbed convection-diffusion problems have been proposed in [9]. An efficient automatic grid control strategy has been proposed in [10] and implemented in the code bvpsuite, a Matlab code designed for the numerical solution of singular boundary value problems [11], which is the new version of the code sbvp [12]. The authors compute the grid density function used for the mesh variation with a feedback control law that adjusts the grid. The importance of generating a smooth mesh is highlighted in [13] were the authors develop a semi-analytic technique called W-grids for generating a suitable smooth mesh for singularly perturbed boundary value problems. The advantages of these meshes over the piecewise uniform Shishkin meshes are also demonstrated with computational experiments.
Recently a comparison of some of the codes available in Matlab for the solution of optimal control problems has been made in [14]. General purpose codes for BVPs are used for the solution of many problems in applications see, for example, [15][16][17] for optimal control problems and path planning, [18][19][20] for numerical simulation, and [21,22] for analysis of bifurcation. General purpose codes are really useful for many real world problems. The users just want to find the solution to their model using functions available in a problem solving environment, like Matlab, Python and R. The R package bvpSolve [23,24], allows the use of the Fortran codes TWPBVP, TWPBVPL and ACDC, and of the collocation codes colsys, colnew and colmod [25,26]. The Matlab environment includes the codes bvp4c, bvp5c and the codes bvptwp [5] and bvpsuite [11]. A collection of Fortran and Fortran 90 codes is available in the website "Testset for BVP solvers", together with many numerical examples [27].
The aim of this paper is to describe the mesh selection strategies that have been included in the code TOM, with the main purpose to make the code comparable in time with the other available codes when the problem is not stiff. Moreover, the computation of the conditioning constants has been modified and theoretically justified, making the code more robust. The paper is structured as follows. In Sect. 2 we review the code TOM, the error estimation criteria and the conditioning parameters used in the mesh selection. In Sect. 3 we present the new mesh selection strategies and in Sect. 4 the method used for solving nonlinear problems. In Sect. 5 we show some numerical experiments with a discussion, comparing the results with some of the available Matlab codes.

The code TOM
We first describe the code TOM for the solution of the linear problem where J (x), v(x) are respectively IR m×m and IR m continuous functions, w ∈ IR m , C a , C b ∈ IR m×m are constant matrices related to the boundary conditions. For linear problems the solution can be easily defined using the fundamental solution matrix Y (x), the Green's function G(x, t) and the nonsingular matrix Q = C a Y (a)+C b Y (b) as: and the continuous conditioning parameters are defined by where κ 1 is related to the perturbation of the boundary conditions and κ 2 is related to the perturbation of the differential equations. The following parameter takes into account both the perturbations see [25,28,29] for a description of the conditioning of linear problems. We also observe that it is possible to compute the conditioning using a different norm and in particular with the scaled 1-norm the parameter corresponding to κ 1 is Moreover we can have information about the stiffness by computing a stiffness ratio related to the conditioning parameters computed in two different norms: The approximation of the conditioning parameters in both norms is one of the main characteristic of the hybrid mesh selection algorithm used in the code TOM and in the code bvptwp. The procedure used in bvptwp has been described in detail in [30] for one step methods. The numerical methods used in the code TOM are however linear multistep methods, so the algorithm needs to be adapted for the computation of κ and κ 2 .
The linear multistep methods used as boundary value methods implemented in the code TOM are all symmetric formulae and in particular are in the class of Top Order Methods (TOM), Extended Trapezoidal Rules (ETR), Extended Trapezoidal Rules of second type (ETR2) and BS methods (BS) [31,32]. The last class turns out to be particularly interesting because generates a collocation method using the B-spline basis. For all the class of methods the discretization of the linear problem (2) with a method of order p generates a linear system of the following form: where, given π = {a = x 0 , x 1 , . . . , x N = b} the discretization of with J := diag(J (x 0 ), . . . , J (x N )) and where, we denote by h i = x i − x i−1 for any x i ∈ π , H := diag(h 1 , · · · , h N ) and A p and B p are banded matrices of size N × (N + 1) whose entries are all the coefficients of the linear multistep methods involved. The number of elements different from zero in each rows are k + 1, where k is the number of steps of the methods. The linear multistep method is used with k 1 < k initial and k 2 = k − k 1 final conditions and the matrices have a band structure with k 1 lower diagonals, except for the first k 1 rows and the last k 2 rows, where the number of elements different from zero is k + 1.
The BS method has a spline continuous extension computed using the algorithm presented in [32] so that the numerical solution satisfies the differential equation at the mesh points. The same continuous extension defines a quasi-interpolation method that, in this release of the code, is used as a continuous extension for the other methods [33]. The coefficients of the BS-methods and of the corresponding quasi-interpolation are computed using the Matlab interface of the QIBSH++ library [34].
Given s p the BS Hermite spline quasi-interpolant associated to ( and it is estimated using only its evaluation at the mesh points: We estimate the local truncation error using a method of higher order q = p + 2 as and we estimate the absolute global error as e E = |y ( p) − y (q) | where y (q) is computed by using an approximated deferred correction method solving the linear system: Depending on the user choice, the code computes either the defect, the local truncation error, or the global error. We consider as approximation of the absolute global error one of the following: We have suitably scaled the defect and the truncation error in order to make the related error comparable with the global error. The code tries to compute a solution with the mixed componentwise relative approximation of the error that satisfies the following inequality where Z is R, T or E depending on the approximation used to estimate the error, atol, a vector of m entries, and rtol are the requested absolute and relative tolerances and the operators are defined componentwise. The estimation of the conditioning parameters need some care. The value of κ 1 can be calculated just computing the first m columns of the inverse of M p , since it depends only on the initial conditions. The values of κ and κ 2 are instead related to the differential problems and they involve all the other block columns of the inverse. For one-step methods it has been proved that the approximation of κ is given by M −1 p 1 [25]. In [29] the authors have defined an efficient algorithm for computing κ, κ 1 and κ 2 .
In order to have a good approximation of the Green function for boundary value methods we scale the linear system (17) by the inverse of the following matrix: where the matrix B p = (B 0 p , B 1 p ) with B 0 p denoting the first column of the matrix and B 1 p all the other columns. The inverse of this function exists for quasi uniform meshes since the methods are symmetric and A k 1 ,k 2 -stable, that is the stability region is equal to the negative half plane (see [31]). After the scaling, the resulting linear problem is: where I e = diag(0 ⊗ I m , I N m ). (2) is discretized by a symmetric BVMs, then

Theorem 1 If the boundary value problem
The proof follows by considering that, denoting by i, j , 0 ≤ i, j ≤ N the blocks of size m of the matrix , the numerical solution is given by: and considering that S p ((0, , the vector B 0 has only k 1 + 1 elements different from zero, so only k 1 + 1 columns of with K a constant that does not depend on N and h = max h i , therefore, y = :,0 w + :,1: This relation shows that the matrix has the same role of the fundamental matrix and of the Green function for the discrete solution and its 1-norm is an O(h) approximation of κ.
Using the results in Theorem 1 we can approximate the value of κ, κ 1 and κ 2 with an algorithm similar to the one used in [29]. If we have not computed the first block column of the inverse we propose, however, a different initialization of the algorithm, saving computational time. This initialization is based on the approximation of the conditioning parameters presented in the next section.

Mesh selection strategies
In this section we introduce the mesh selection strategies that are possible to use in the code TOM. The numerical experiments will show that are all useful, depending on the problem under consideration. The variants will be denoted by taking a combination of three or four letters, the first two letters will be HS (High stiff), MS (Middle Stiff), NS (Non Stiff) and give information to the class of problems for which the mesh is suggested. If we compute the conditioning parameters by computing the first block column of the inverse we add the letter (C). If we estimate them using the computed solution we add the letter (S). The three new parameters computed using the approximation y of the solution in the current mesh π are called κ y , γ y , σ y . We have The parameter κ y gives indication about the maximum norm of the approximated solution. The parameter γ y gives information about the integral of the absolute value of the functions associated to each components of y. If the values of κ y and γ y with different meshes are stabilized, then we have a good preliminary approximation of the solution and a good mesh. The value of σ y gives indication about the stiffness of the problem: it is the ratio between the maximum absolute value of each component of the approximated solution and the corresponding approximated scaled integral. This value is very high if the solution has rapid variation. In the mesh selection κ y , γ y , σ y are used in substitution of the approximations of κ 1 , γ 1 and σ obtained computing the first block column of the matrix , the inverse of S −1 p M p .
When the mesh selection is denoted with only three letters we are working with nonstiff problems and the conditioning parameters are not used in the monitor function. The last letter is E if we use an estimation of the global error (E E ), R if we use a scaled estimation of the residual (E R ) and T if we use a scaled estimation of the local truncation error (E T ).
The mesh selection denoted by HSCE is the hybrid mesh selection introduced in [2] where for high stiff problems, we use the order two method to select the mesh and, when a good approximation is computed, the selected order p is used. The mesh selection is based on the computation of the first block column of the matrix in combination with the error estimation. In particular we choose as monitor function If the conditioning parameters are stabilized, that is the relative error between the approximation computed with two different meshes has a relative error less than 10 −2 , also the error estimation is used in the mesh selection. Moreover the border and internal values of the biggest values of the monitor function are multiplied by a factor, in order to add more points in the regions of rapid variation, we use 100 for the values at the two extremes of the interval and 125 for the one in the center of the interval. This monitor function is extremely useful for high stiff problems, but when the problem is only middly stiff the procedure used for changing the order does not work well and for relatively simple problems the code uses the order two, thus requiring much more points. The first new mesh selection added is the one called MSCE for middle stiff problems, that does not use the order two method to selecting an initial appropriate mesh, but uses the same order in all the steps. The multiplication by a factor of the maximum value of the monitor function is useful only for high stiff or middly stiff problems, for general problems the monitor function gives soon information about the error, so the third considered new mesh selection is the one called NSCE, that uses the monitor function as it is. The HSCR, MSCR, NSCR and HSCT, MSCT, NSCT mesh selection strategies use respectively the residual and the truncation error in the mesh selection.
Since the computation of the first block column of the matrix requires the solution of m linear systems, we have introduced a set of mesh selection that uses κ y , γ y and σ y instead of κ 1 , γ 1 and σ . Moreover the new monitor function φ s (x) is a scaled linear combination of the following monitor functions φ 1 (x i ) = y i 1 + y i−1 1 and φ 2 (x i ) = | y i ∞ − y i−1 ∞ |, where y is an approximation of the derivative of the approximated solution at the mesh points of the same order p . It is interesting to observe that, if we have computed a good approximation of the solution, we are going to add more points when its derivative is large, and these are the regions of rapid variation. We have implemented the same variants of the mesh selection based on conditioning, just considering these new monitor functions. The new conditioning parameters based on the computed solution are also used to decide if the parameters are stabilized and we can use the error in the mesh selection. The resulting mesh selection strategies are denoted respectively as HSSE, MSSE and NSSE if we use the global error. A preliminary version of MSSE has been used in [14] for the solution of optimal control problems. The HSSR, MSSR, NSSR, HSST, MSST, NSST mesh selection strategies use respectively the residual and the truncation error in the mesh selection.

Algorithm 1 Mesh selection
r = 0 s (r ) approximated initial solution, π r initial mesh, TMF = string with the chosen monitor function, Z = (R or T or E). while the maximum estimated error max(E Z ) is higher than the input tolerances do compute the values of the new approximate solution in π r by solving the linear problem (21) using s (r ) compute an approximation of the componentwise mixed relative approximation of the error E Z compute an estimate of the conditioning parameters κ, κ 1 , κ 2 , γ 1 , σ (optional for TMF= *SS* or TMF= NS*) compute κ y , γ y , σ The mixed relative error (15) is approximated using either the global error E E (last letter E), the scaled residual E R (last letter R) or the scaled local truncation error E T (last letter T).
The algorithm for the hybrid mesh selection starts using as monitor function φ c or φ s . The monitor function is firstly used to move the mesh points and after (if needed) to add or remove points. If the estimation of error is considered in the mesh selection, the monitor function used to move the mesh points is a linear combination of the two, the one used to add and/or remove points is only based on the error. The algorithm for linear problems is briefly described in Algorithm 1. The empirical parameters used to decide how to add and remove points are omitted, they depend on the values of the conditioning parameters.

Nonlinear iteration
For linear problems the most important step is the mesh selection, once we have a good information on the error we need just to move mesh points, or to add and remove them, in order to have the error less than the input tolerances. In the nonlinear case we have another big issue: the convergence of the method used for the solution of the nonlinear equations. By the author knowledge usually general purposes codes implement the Newton iteration or some of its variants that enlarge the region of convergence. For difficult problems, like singularly perturbed one, a continuation strategy is considered. Some attempt to make the continuation strategy automatic has been made in the codes colmod and ACDC [26]. In general, the mesh is adapted only if the non linear method is able to compute an approximated solution. The same holds for the quasi-linearization strategy implemented in colsys, colnew [25,35,36], the mesh is the same for all the quasi-linearization steps. The code TOM implements the quasi-linearization algorithm described in [32], solving at each step a linear differential boundary value problem, and moving to the next linear problem when some properties are satisfied. This is the most important decision, because solving the linear problem with high accuracy is not required, but if the linear problem is not sufficiently accurate the convergence is not assured in the discrete case. Given y (n) (x) an approximation of the continuous linear problem at step n, the linear continuous problem that we solve at the next step is given by: where The discrete problem obtained after the discretization with a BVM is the following: where f = ( f (x 0 , y n 0 )), · · · , f (x N , y n N )). In this way, if we use a fixed mesh for each iteration, the algorithm is equivalent to the Newton method applied to the discretized equations. If we change the mesh, the vector y n is computed by evaluating the spline quasi-interpolant of the previous computed solution, that is y n i = s (n) (x i ). A difference with respect to the implementation described in [32] is given by the criteria to decide when to start the next step. Given s (n+1) the spline quasi-interpolant associated to the current approximation the computed approximation of the residual in the mesh points is approximated as: Following [37,38] we check that the relative residual choosing η n+1 < min(0.01, (s (n+1) ) − f (x, s (n+1) ) L 2 ). This assures the quadratic local convergence of Newton method for linear boundary conditions. To make the procedure more robust we add more checking before to accept the solution: • the relative error E Z must be less than 1 using as tolerances for the linear problem: rtol lin := min(0.05, max(rtol, f l (x, y (n+1) (x)) − f (x,ȳ (n+1) ) L 2 )), atol lin := min(10 −3 , max(atol, f l (x, y (n+1) (x)) − f (x,ȳ (n+1) ) L 2 )); • the conditioning parameters κ y , γ y , σ y (and if available κ 1 , γ 1 , σ ) are stabilized; The nonlinear iteration is described by Algorithm 2.
Algorithm 2 nonlinear iteration n = 0, r = 0 s (n) approximate initial solution π r initial mesh while the error of the nonlinear problem is higher than the input tolerance do compute the values of the new approximate solution in π r by solving the linear problem (21) using s (n) compute an approximation of the error (E E or E R or E T ) of the linear problem and the conditioning parameters if the linear problem satisfies the stopping criteria then compute the spline extension s (n+1) , n = n + 1, π 0 = π r , r = 0 else change the mesh using Algorithm 1, r = r + 1 end if end while It is interesting to observe that the mesh selection and the nonlinear iteration are strictly connected. An important role in checking if the numerical solution has been approximated in the correct way is given by the evaluation of the conditioning parameters, in fact if they are stabilized we can suppose that we have reached a good approximation of the continuous solution of the current linear problem.

Numerical examples
In this section we present some numerical examples to show that the appropriate mesh selection is important to have a code that works in difficult situations and that the use of conditioning and the mesh selection embedded in the nonlinear iteration can help the convergence of the nonlinear method. We also observe that if the problem is not stiff, the use of the conditioning parameters is not necessary and the method for solving the nonlinear problem usually converges without any issue.
In the following we compare the code TOM with the Matlab codes bvp4c, bvp5c, bvptwp. Even if in the code TOM it is possible to choose the order and the method, we use for all the examples the orders 4 and 6 and the BS method. In summary the codes used in the experiments are: For all the compared codes we choose as maximum number of mesh points N=20000. All the experiments have been performed using an Intel Core I9 10 core 3.6 GHz with Matlab R2021b for MacOS. The work precision diagrams are computed choosing the tolerances atol = rtol ranging from 10 −3 to 10 −8 and by estimating the mixed error as by computing a more precise solution y * with the code twpbvp_m with a doubled number of fixed points and a relative tolerance 10 times smaller. We choose for comparison some linear and nonlinear examples available in the website testset for BVP solvers [27], that contains a list of singularly perturbed test problems in Fortran and in Matlab. The Matlab version has been used to check the behaviour of the code bvptwp in [5].
The first example is the problem bvpT1, a linear second order boundary value problem, written as a system of first order equations. This problem has a boundary layer in a, but if we choose λ = 10 −2 the solution does not have any rapid variations. In this case it is not useful the use of the conditioning in the mesh selection and the code is more efficient if we use the mesh selection for nonstiff problems.  shows that the NSE, NSR, NST, NSSE, NSST and NSSR mesh selections give similar results, so they can be used indifferently, HSCE gives the highest computational time. Moreover we observe that the behaviour of the compared methods is mainly dictated by the order of the numerical method used and not by the mesh selection. In fact the methods with a smaller execution time are the ones with variable orders 4, 6 and 8 (twpbvp_m and twpbvp_l).
The second test problem is the problem bvpT6, a linear problem with a turning point. We choose λ = 10 −2 and we observe that the behaviour of the codes is similar to the one of the first example, the mesh selection HSCE being the most expensive (see Fig. 2). This is confirmed by the values of the stiffness parameter reported in Table 1, that shows that the problem is nonstiff.

Conclusion
In this paper we have presented the new mesh selection strategies implemented in the code TOM, together with a proper approximation of the conditioning constants of the discretized problem when the underlying method is not one step. The results show that for linear problems the hybrid mesh selections are efficient only for the stiff case, this justifies the introduction of standard mesh selection for the solution of nonstiff problems. Moreover, for nonlinear problems, the conditioning parameters, or their approximations, allow to implement a method for their solution that is efficient for some problems for which the standard Newton type iteration methods do not converge using a fixed mesh. The proposed quasi-linearization strategy for the solution of non linear problems needs some more theoretical studies in order to have a robust technique for choosing the acceptance of the solution of an intermediate linear system. Moreover it would be interesting to see how the generation of a final smooth grid using the algorithm presented in [13] would impact on the performance of the code.