Balancing truncation and round-off errors in FEM One-dimensional analysis

In finite element methods, the accuracy of the solution cannot increase indefinitely since the round-off error related to limited computer precision increases when the number of degrees of freedom (DoFs) is large enough. Because a priori information of the highest attainable accuracy is of great interest, we construct an innovative method to obtain the highest attainable accuracy given the order of the elements. In this method, the truncation error is extrapolated when it converges at the asymptotic rate, and the bound of the round-off error follows from a generically valid error estimate, obtained and validated through extensive numerical experiments. The highest attainable accuracy is obtained by minimizing the sum of these two types of errors. We validate this method using a one-dimensional Helmholtz equation in space. It shows that the highest attainable accuracy can be accurately predicted, and the CPU time required is much smaller compared with that using successive grid refinement.


Introduction
Many problems in engineering sciences and industry are modelled mathematically by initial-boundary value problems comprising systems of coupled, nonlinear partial and/or ordinary differential equations. These problems often consider complex geometries, with initial and/or boundary conditions that depend on measured data [1]. In some applications, not only the solution, but also its derivatives are of interest [1,2]. For many problems of practical interest, analytical or semi-analytical solutions are not available, and hence one has to resort to numerical solution methods, such as the finite difference, finite volume, and finite element methods. The latter will be adopted throughout this paper and applied to one-dimensional boundary value problems.
The accuracy of the numerically obtained solution is influenced by many sources of errors [3]: firstly, modelling errors in the set-up of the models, such as the simplification of realistic domains and governing equations and the approximation of initial and boundary conditions; next, truncation errors due to the discretization of the computational domain and the use of basis functions for the function spaces defined on it; then, round-off errors due to the adoption of finite-precision computer arithmetics, rather than exact arithmetics; finally, iteration errors resulting from the artificially controlled tolerance of iterative solvers.
One tacitly assumes that most errors are well-balanced and/or negligibly small. In this paper, the focus is on the truncation error (E T ) and the round-off (E R ), by considering idealized problems, which do not introduce modelling errors, and using a direct solver, which avoids the introduction of iterative errors. In particular, the round-off error is often ignored based on the argument that it will be 'sufficiently small' if just IEEE-754 double-precision floating-point arithmetics [4] are adopted. Therefore, to improve the accuracy, i.e. to decrease E T , one often reduces the mesh width (h-refinement), increases the approximation order (p-refinement), or applies both strategies simultaneously (hp-refinement) [5,6].
The common characteristic of these methods is to increase the number of degrees of freedom (DoFs). However, E R increases with the number of DoFs, and dominates the total error if more and more DoFs are employed [7,8]. While typically an impractically large number of DoFs is required for E R to dominate the total error if low(est)order approximations are used, the number can be very small if high-order approximations are adopted, which are nowadays becoming more and more popular. This shift to higher order approximations makes the results more prone to be dominated by round-off errors. Despite this alarming observation, to the authors' best knowledge, only very few publications address the impact of accumulated round-off errors on the overall accuracy of the final solution or take them into account explicitly in the error-estimation procedure. The general rule of thumb is still to perform as many h-refinements as possible considering the available computer hardware.
The aim of this paper is to systematically analyse the influence of the round-off error on the total error when using h-refinements for different orders of p. Not only the solution but also its first and second derivatives are investigated for one-dimensional, second order model problems, assuming the second derivative exists in the weak sense [9]. Both the standard finite element method (FEM) and the mixed FEM [10] are analysed for multiple p's. Furthermore, the following factors are considered: types of boundary conditions and methods of implementing them, choices and configurations of the linear system solver, orders of magnitude of the variables and coefficients. Based on the statistics of the evolution of the round-off error, we propose an algorithm to predict the best accuracy E min that occurs when the sum of E T and E R is the smallest, and the corresponding number of DoFs (N opt ).
The paper is organized as follows. The model problem, finite element formulation and numerical implementation are described in Section 2. The approach to predicting E min is illustrated in Section 3. The statistics on the evolution of the round-off error are given in Section 4. The algorithm for realizing the approach is put forward in Section 5, followed by its validation by a Helmholtz problem in Section 6. The conclusions are drawn in Section 7.

Model problem
Consider the following one-dimensional second-order differential equation: with u denoting the unknown variable, which can either be real or complex, f (x) ∈ L 2 (I) a prescribed right-hand side, and d(x) and r(x) continuous coefficient functions. By choosing d(x) = 1 and r(x) = 0, Eq. (1) reduces to the Poisson equation; for d(x) > 0 and not constant, the diffusion equation is found when r(x) = 0, and the Helmholtz equation [11] is found when r(x) ̸ = 0. The boundary conditions are u(x) = g(x) on Γ D and d(x)u x = h(x) on Γ N , where Γ D and Γ N are the boundaries where Dirichlet and Neumann boundary conditions are imposed, respectively.

Finite element formulation
For convenience, we introduce two inner products [12]: where f 1 (x) and f 2 (x) are continuous functions defined on the unit interval I, Γ denotes the boundary of I, and x 0 denotes the values of x on Γ .
2 Imposing the Dirichlet boundary conditions in the weak sense [13], the weak form reads: Weak form 2 Find u ∈ H 1 (I) such that: where ρ is a positive value that serves as the penalty parameter. (4) In both forms, η denotes the test function, n is equal to 1 at x = 1, and −1 at x = 0; the terms on the right-hand sides consist of information of Neumann boundary conditions which vanishes if no Neumann boundary conditions are prescribed. We approximate u by a linear combination of a finite number of basis functions: Here, p is the element degree, m is the number of DoFs, which equals p × t + 1, with t denoting the total number of grid cells; u i 's are the values of u i (x j ) = δ ij , with x j denoting the support point. This type of element will be referred to as P p . Taking η equal to ϕ (p) k , k = 1, 2, . . . , m, the resulting linear system of equations is read as where A is the stiffness matrix, F the right-hand side and U = [u 1 , . . . , u m ] ⊤ .

The mixed FEM
As a first step, we introduce the auxiliary variable allowing Eq. (1) to be rewritten as The weak form of Eq. (1) using the mixed FEM, derived in Appendix A.2, is given by: Weak form 3 Find v ∈ H 1 N (I) and u ∈ L 2 (I) such that: with H 1 In this form, w and q denote the test functions of v and u, respectively, and n has the same interpretation as before. We approximate v and u by: where m is the number of DoFs for v (p) h , which is equal to p × t + 1, v i 's are the values of v (p) h at the DoFs, and ϕ (p) i 's are of the same type of basis functions used in Eq. (5); u sj 's are the values of u (p−1) h at the DoFs, ψ (p−1) j 's are discontinuous Lagrange basis functions, which implies that two independent u sj 's have been assigned at the cell interfaces. This pair of elements will be referred to as P p /P disc p−1 . Replacing w and q by ϕ (p) k , k = 1, 2, . . . , p×t +1, and ψ (p−1) e , e = 1, 2, . . . , p×t, respectively, the resulting coupled linear system of equations that has to be solved is read as: where the mass matrix M, discrete gradient operator B, and its transpose, the discrete divergence operator B ⊤ , comprise the left-hand side; G and H are the components of the right-hand side; V = [v 1 , . . . , v m ] ⊤ and U = [u 11 , . . . , u 1p , . . . , u t1 , . . . , u tp ] ⊤ , respectively. For the sake of readability, we will drop the superscript (p) or (p − 1) whenever the approximation order is clear from the context.

Solution technique
All results are computed in IEEE-754 double precision [4] using the deal.II finite element library [14]. Unless stated otherwise, the computational mesh is obtained by globally refining a single element that covers the interval I, and the Dirichlet boundary conditions are imposed strongly. The former means that, when the solution is real valued, the number of DoFs equals 2 R × p + 1 using the standard FEM and 2 × 2 R × p + 1 using the mixed FEM, at the Rth refinement; when the solution is complex valued, the above numbers double since deal.II does not provide native support for complex-valued problems and, hence, all components need to be split into their real and imaginary parts.
To compute the occurring integrals, sufficiently accurate Gaussian quadrature formulas are used. To solve the systems of equations, the UMFPACK solver [15], which implements the multi-frontal LU factorization approach, is used unless stated otherwise. This solver results in relatively fast computations of the problems considered in this paper, and prevents the iteration errors of iterative solvers. The derivatives of the numerical solution, which are u h,x and u h,xx in the standard FEM and only v h,x in the mixed FEM, are computed in the classical finite element manner, e.g.
an approximation to u x using standard FEM. Note that, each differentiation decreases the element degree by one.

Error estimation
For the numerical results var h , where var can be u, u x and u xx of the standard FEM, and u, v and v x of the mixed FEM, the error measured in the L 2 norm is used. It is defined as when the exact solution var exc is available, and [16] otherwise, where var h/2 is the numerical solution computed on a mesh once refined with grid size h/2.

Convergence of the solution
When the number of DoFs is relatively large, but the round-off error does not exceed the truncation error, the error converges at a fixed rate, known as asymptotic convergence rate, of which the value is one order higher than the approximation order [6]. In practice, the convergence rate in the numerical experiments can be calculated from either using Eq. (11b).

Approach to predicting the highest attainable accuracy
A conceptual sketch of E h against the number of DoFs (N h ) in a log-log plot can be found in Fig. 1, also see [17]. When N h is relatively small (N h < N c ), E h does not decrease at the aforementioned asymptotic order of convergence, and only when N h is large enough (N c ⩽ N h < N opt ) this asymptotic order of convergence is attained. The transition from the first phase, denoted by black circles, to the second phase, denoted by green circles, is usually fast, cf. Fig. 1. E h in both phases is controlled by the truncation error E T ; in the second phase E h can be represented by Here α T is the offset, and β T is the slope of the line approximating E h , equalling the asymptotic order of convergence, see Appendix B for the proof. Note that, α T can be inverted by using at the beginning of the second phase, where E c equals the corresponding error E h .

Table 1
Description of the evolution of E h .
, the round-off error E R starts to dominate and E h increases, illustrated by orange circles. At this phase, the slope of the line approximating E h , denoted by β R , tends to be fixed [8,18]. The parameter β R and the associated offset, denoted by α R , are investigated in detail in Section 4. As will be shown there, α R and β R are fixed constants, which allow us to estimate E h as In summary, the evolution of E h is described in Table 1, and depicted in Fig. 1.
Since the evolution of E h (= E T + E R ) is known after entering the second phase, by solving we can predict the optimal number of DoFs and hence, the highest attainable accuracy is given by

Results
In this section, we assess the general values of α R and β R . We start with a benchmark Poisson equation, for which the influences of solution strategies and boundary conditions are investigated, and then consider more general parameters for the Poisson equation, as well as the diffusion and Helmholtz equations.

Benchmark Poisson equation
We consider the Poisson equation with . The boundary conditions are imposed as follows: The error E h of the resulting solution u, and its first and second derivatives, using the standard FEM and the mixed FEM, with p ranging from 1 to 5 are in Fig. 2 and Fig. 3, respectively. In these figures, α R and β R are denoted. It is found that, for all dependent variables and their derivatives, the values of α R and β R of different element degrees are the same. The statistics of the former can be found in Fig. 4(a). The values of β R only depend on the FEM method, and are 1 using the mixed FEM and 2 using the standard FEM. Notably, α R is of order 10 −16 , which is as expected when using double precision, and tends to increase slightly with increasing order of derivative. Furthermore, α R of the mixed FEM is smaller than that of the standard FEM.
For larger p, E T decreases faster such that smaller E min can be obtained, see Fig. 4(b). In general, smaller E min can be obtained using the mixed FEM compared to using the standard FEM.
In Sections 4.1.1-4.1.2, the sensitivity of the above results will be investigated, using P 2 elements for the standard FEM and P 4 /P disc 3 elements for the mixed FEM.

Solution strategy
In this section, we investigate the influence of the solution strategy on the accuracy of the numerical solution. In particular, we compare the outcome when applying the direct solver UMFPACK with that of using the iterative Conjugate Gradient (CG) method [19], which can be applied when the left-hand side, e.g. A in Eq. (6), is symmetric and positive definite. The tolerance of the CG solver is a small number, denoted by tol prm , in the standard FEM, and the product of a small number tol prm and the L 2 norm of the corresponding right-hand side in the mixed FEM. When the L 2 norm of the residual, e.g. ∥F − Au∥ 2 in Eq. (6), is smaller than the tolerance, the iteration is stopped.
The standard FEM. For tol prm = 10 −10 and 10 −4 , the absolute errors of u, u x and u xx using the CG solver are shown in When tol prm is adequately small, i.e. tol prm = 10 −10 , the round-off error for the solution and the first derivative using the CG solver is the same with that using the UMFPACK solver; the round-off error for the second derivative using the CG solver increases faster than that using the UMFPACK solver. When tol prm is too large, i.e. tol prm = 10 −4 , the error contribution due to the iterative solver dominates both truncation and round-off errors for an intermediate number of DoFs. The mixed FEM. Since the resulting matrix Eq. (10) is indefinite, a widely used alternative is to decouple the fully coupled monolithic approach using Schur's complement and solve both equations in segregated manner, i.e. Eq. (18a) is solved in the first place to obtain U, which is then substituted into Eq. (18b) to obtain V .
Eq. (18a) involves the term M −1 G on the right-hand side, which is computed by solving the auxiliary linear system MY = G by using either the UMFPACK or the CG solver. The same options are available for solving Eq. (18b).
The difficulty in solving Eq. (18a) lies in not assembling the Schur complement matrix explicitly since it comprises M −1 . The CG solver only makes use of matrix-vector products of the form (B ⊤ M −1 B)W , which can be computed by the following three-step algorithm: X = BW , MY = X and Z = B ⊤ Y . As before, the linear system MY = X can be solved by the UMFPACK or the CG solver.
We first investigate the influence of tol prm of the CG solver on the accuracy of the solutions when the left-hand side is B ⊤ M −1 B. In this case, the UMFPACK solver is used to solve the matrix equations when the left-hand side is M. For tol prm being 10 −16 and 10 −10 , the results are shown in Fig. 6, in comparison with that obtained from solving the monolithic Eq. (10) directly using the UMFPACK solver. It shows that, for the problem at hand, the monolithic solution approach yields by far the most accurate solution and derivative values. The round-off error for v x increases fastest using the Schur complement approach even though tol prm is sufficiently small, i.e. tol prm = 10 −16 , which makes the highest attainable accuracy much lower. When tol prm is less strict, i.e. tol prm = 10 −10 , the iteration error dominates the total error instead of the round-off error.
Next, we investigate the influence of tol prm of the CG solver when the left-hand side is M. In this case, the CG solver with tol prm being 10 −16 is used to solve the matrix equation with the left-hand side being B ⊤ M −1 B. For tol prm being 10 −16 and 10 −10 , the results are shown in Fig. 7, in comparison with that obtained from solving the monolithic Eq. (10) directly using the UMFPACK solver. It also shows that, when the tolerance is less strict, i.e. tol prm = 10 −10 , the iteration error dominates the total error before the round-off error.
In summary, in comparison with the UMFPACK solver, the CG solver gives the same accuracy for u and u x but less accuracy for u xx using the standard FEM, and gives smaller accuracy for all the three variables using the mixed FEM when   tol prm is strict enough; the CG solver introduces iteration errors for both the standard FEM and the mixed FEM when tol prm is less strict. Thereby, we continue with the UMFPACK solver in our algorithm to obtain higher accuracy.

Boundary condition
In this section, two aspects of the influence of the boundary conditions on the round-off error are investigated: first the method of implementing the Dirichlet boundary conditions, and secondly types of boundary conditions. For the first aspect, using Weak form 2 for ρ = 50 and 10 6 , the errors are depicted in Fig. 8, in comparison with that using Weak form 1. As can be seen, both the weak imposition and the strong imposition of the Dirichlet boundary condition yield the same trend line for the round-off error for the solution and its derivatives, and the magnitude of the penalty parameter in the weak imposition makes no difference. In addition, small penalty parameters might lead to larger truncation errors for u, but the difference diminishes when the penalty parameter is large enough.
To construct the problem for the second aspect, the Dirichlet boundary condition at the left boundary (x = 0) is kept while the Dirichlet boundary condition at the right boundary (x = 1) is replaced by the Neumann boundary condition u x (1) = −e −1/4 , leading to the same solution and derivative profiles.
The standard FEM. Using the standard FEM, the offsets α R for the two types of boundary conditions are depicted in Fig. 9(a). For the Dirichlet/Neumann boundary condition, the offsets α R for u and u x are slightly larger than that for the Dirichlet/Dirichlet boundary condition by a factor of 3.5 and 2, respectively. The offsets α R for u xx are identical for the two types of boundary conditions.
The mixed FEM. Using the mixed FEM, the offsets α R for the two types of boundary conditions are depicted in Fig. 9(b). As can be seen, the type of boundary conditions plays a more important role for α R for the solution than α R for other variables.
In summary, α R are relatively independent of the variations in the type of boundary conditions and the method Dirichlet boundary conditions are implemented, which is an important prerequisite for our a posteriori refinement strategy to be applicable for a wide range of problems. However, since the asymptotic behaviour of the truncation error is essential in our algorithm, we continue with the strong imposition of the Dirichlet boundary condition. In Sections 4.2-4.3, where the influences of u(x), d(x) and r(x) are investigated, we only consider the Dirichlet boundary conditions, and still use P 2 elements for the standard FEM and P 4 /P disc 3 elements for the mixed FEM.

General Poisson equation
In this section, we will again consider the Poisson equation, but now focus on the influence of the order of magnitude of the solution and right-hand side on α R . To cover a wide range of scenarios, we choose the cases shown in Table 2.
Each case contains a coefficient c i , i = 1, 2, . . . , 5, which is varied over several orders of magnitude so that the L 2 norm of the solution, denoted by ∥u∥ 2 , and the L 2 norm of the right-hand side, denoted by ∥f ∥ 2 , extend over a wide range of magnitudes. Fig. 10 gives an overview of the distribution of ∥u∥ 2 and ∥f ∥ 2 for Cases 1-4. For these Poisson problems, the error basically evolves according to that shown in Fig. 1. To summarize, β R is 2 using the standard FEM and 1 using the mixed FEM. α R against ∥u∥ 2 is shown in Fig. 11(a) for the standard FEM; α R against ∥u∥ 2 or ∥v∥ 2 , depending on the dependent variable of interest, is shown in Fig. 11(b) for the mixed FEM.
For all the variables using both the standard FEM and the mixed FEM, the distribution of α R can be approximated by a straight line. This implies that α R has a power-law relation with the horizontal coordinate. The power term, represented by the slope of the line, is 1 for all the relations; the constant term, represented by the intercept of the line, is denoted in the figure. The latter reads 2e-17, 5e-17 and 5e-16 for u, u x and u xx , respectively, using the standard FEM, and 1e-18, 1e-16 and 5e-16 for u, v and v x , respectively, using the mixed FEM. Note that, the intercept is the value of α R when the horizontal coordinate is 1. Therefore, α R has a linear relation with the L 2 norm of the dependent variable. We express α R as the product of a constant and the magnitude of ∥u∥ 2 or ∥v∥ 2 that is shown in Table 3.  Table 2.  Table 3 α R in terms of the product of a constant and the L 2 norm ∥u∥ 2 or ∥v∥ 2 , dependent on the specific dependent variable, for one-dimensional Poisson equations.

General diffusion and Helmholtz equation
The coefficient d(x) is taken from Table 4 for the diffusion equations, and d(x) = 1 is taken for the Helmholtz equations, with r(x) taken from Table 5. The analytical solution u is the same as the benchmark Poisson equation, of which ∥u∥ 2 = 0.92 and ∥v∥ 2 = 0.5.
Again the errors show the same behaviour as that shown in Fig. 1. The coefficient β R is 2 using the standard FEM and 1 using the mixed FEM.  The parameter α R against ∥d∥ 2 is shown in Fig. 12 for the diffusion equations. Using the standard FEM, α R is independent of ∥d∥ 2 for u, u x and u xx , with the upper bound reading 2e-17, 5e-17 and 1e-15, respectively. Using the mixed FEM, α R is independent of ∥d∥ 2 for u, but is linearly dependent on ∥d∥ 2 for both v and v x . The upper bound for u reads 5e-17; the intercept for v and v x reads 2e-16 and 5e-16, respectively. The parameter α R against ∥r∥ 2 is shown in Fig. 13 for the Helmholtz equations. Visibly, α R is independent of ∥r∥ 2 using both the standard FEM and the mixed FEM. Using the standard FEM, α R reads 2e-17, 5e-17 and 2e-16 for u, u x and u xx , respectively; using the mixed FEM, α R reads 1e-19, 1e-16 and 2e-16 for u, v and v x , respectively.
In Figs. 12-13, we do not take the influence of ∥u∥ 2 or∥v∥ 2 into account. Since ∥u∥ 2 is of order 1, we omit its influence on α R for the standard FEM and on α R of u and v for the mixed FEM. Since ∥v∥ 2 is half of 1, we multiply the intercept of α R of v x in Fig. 12(b) and the upper bound of α R of v x in Fig. 13(b) by 2, for validating the constant term of α R in Table 3.
Comparing the resulting α R of the diffusion equations with the constant term of α R in Table 3, for the standard FEM, the constant term of α R of u and u x does not change, while that of u xx increases to 1e-15; for the mixed FEM, the constant term of α R of u increases to 5e-17, that of v becomes 2e-16 ×∥d∥ 2 , and that of v x increases to 1e-15. Therefore, we amend α R in terms of the product of a constant and the L 2 norm ∥u∥ 2 or ∥v∥ 2 , dependent on the specific dependent variable, for one-dimensional second order differential equations.   Table 3 to be that shown in Table 6. Note that, for the constant term of α R of u of the mixed FEM, considering most of its values of the diffusion equations are smaller than 2e-17, which is the value we obtain for the standard FEM, we choose 2e-17 as its value. Furthermore, since the constant term of α R of the Helmholtz equations is smaller than that in Table 6, its effect can be ignored. Summarizing Sections 4.2-4.3, the order of magnitude of the coefficient r(x) is basically irrelevant; the order of magnitude of the dependent variable and the coefficient d(x) can be mitigated when their values are known. Specifically, the value of the order of magnitude of the dependent variable can be obtained from a built-in function in our algorithm, and that of d(x) can be obtained by carrying out simple integration.

A posteriori algorithm for finding the optimal number of degrees of freedom
Based on the validation experiments from the previous section, we introduce a novel a posteriori algorithm for determining E min and its associated N opt for the solution and its first and second derivatives without performing brute-force mesh refinement. We call the algorithm DoFinder.
In DoFinder, we define the following coefficients and use them in the steps given below.
-a minimal number of h-refinements before carrying out 'NORMALIZATION' and 'PREDICTION', denoted by R min , with the following default values: We choose this parameter mainly because the error might increase, or decrease faster than the asymptotic order of convergence for coarse refinements, especially for lower-order elements. -the allowed maximum N h : 10 8 , denoted by N max .
-a stopping criterion c s for seeking the L 2 norm of the dependent variable, of which the value is 0.001 by default. We choose this parameter because the analytical solution does not exist for most practical problems. -a relaxation coefficient c r for seeking the asymptotic order of convergence, with the following default values: c r = ⎧ ⎨ ⎩ 0.9 for p < 4, 0.7 for 4 ⩽ p < 10, 0.5 otherwise. (20) -the offset α R , see Table 6 for the default values.
The procedure of DoFinder consists of four steps, which are explained below: Step-1. 'INPUT '. In this step, the custom input shown in Table 7 has to be provided.
Step-2. 'NORMALIZATION'. The function of this step is to find the L 2 norm of the dependent variable, in which elements of degree p min are used. The specific procedure can be found in Algorithm 1. Step-3. 'PREDICTION'. This step finds E min for each var and p of interest, as illustrated in Fig. 1. The procedure for carrying out this step can be found in Algorithm 2. Step-4. 'OUTPUT '. In this step, we output E min , N opt , etc., obtained from Step-3.

Validation
In what follows, we validate the strategy discussed in Section 3 by using the following Helmholtz problem: with homogeneous Dirichlet and Neumann boundary conditions imposed as follows: u(0) = 0 and u x (1) = 0. Both the standard FEM and the mixed FEM are investigated, with the element degree p taken in {1, 2, . . . , 5}. Variables u, u x and u xx using the standard FEM and u, v and v x using the mixed FEM are investigated.

Accuracy analysis
Using DoFinder and the brute-force approach, E min 's are compared in Fig. 14. As can be seen, E min can be predicted correctly. Note that, for p = 1 using the mixed FEM, we cannot obtain E min of the brute-force approach because of the limited hardware. The associated N opt using DoFinder is shown in Fig. 15, and it is used to compute the solution of the optimal grid in Section 6.2.

Efficiency analysis
The CPU time required by DoFinder and the brute-force approach is shown in Fig. 16 and Fig. 17, respectively. In general, the CPU time associated with the brute-force approach decreases with increasing element degree. In comparison with the CPU time of the brute-force approach, the CPU time of DoFinder is negligible.
Furthermore, to obtain the solution on the optimal grid based on DoFinder, the computation time can be saved much in comparison with that using the brute-force approach. Fig. 18 gives the percentage of the CPU time saved by DoFinder. It shows a saving of the CPU time of around 70% for both the standard FEM and the mixed FEM.

Further development
For a given tolerance of a variable, denoted by tol var , we are able to quickly select the available FEM method, the associated minimal available p and minimal available number of DoFs. For example, when tol u = 1e-9, we can only achieve it using the mixed FEM, the associated minimal available p is 2, and minimal available number of DoFs is 638876.

Conclusions
A novel approach is presented to predict the highest attainable accuracy for second order, ordinary differential equations using the finite element method. In contrast to the brute-force approach, which uses successive h-refinements, this approach uses only a few coarse grid refinements to reach the region of asymptotic convergence. This approach is viable for the solution and its first and second derivatives, for both the standard FEM and the mixed FEM, for all element degrees. The algorithm for implementing the approach shows that the highest attainable accuracy can be accurately predicted using negligible CPU time and the CPU time to obtain this solution explicitly is significantly reduced: to compute the solution with the highest attainable accuracy using our approach results in a reduction of CPU time by around 70% for both the standard FEM and the mixed FEM. Future research will focus on the extension of this approach to 2D second order problems, where the influence of the linear system solver, local mesh refinement and boundary conditions might be significantly different from 1D problems. where C 1 = Cp p+1 . Therefore, the slope is β T = p + 1