Uniform and high-order discretization schemes for Sturm-Liouville problems via Fer streamers

The current paper concerns the uniform and high-order discretization of the novel approach to the computation of Sturm–Liouville problems via Fer streamers, put forth in Ramos and Iserles (Numer. Math. 131(3), 541—565 2015). In particular, the discretization schemes are shown to enjoy large step sizes uniform over the entire eigenvalue range and tight error estimates uniform for every eigenvalue. They are made explicit for global orders 4,7,10. In addition, the present paper provides total error estimates that quantify the interplay between the truncation and the discretization in the approach by Fer streamers.

The current paper concerns the discretization of the novel approach to the computation of Sturm-Liouville problems via Fer streamers, introduced in [14].
The motivation to discretize the new approach via Fer streamers stems from the local and global error estimates in [14] which guarantee large step sizes uniform over the entire eigenvalue range and tight error estimates uniform for every eigenvalue. The current paper shows how to retain these advantageous features under discretization.

Uniform versus asymptotic error estimates
It is well known that Eqs. 1-2 possess a unique countable family of solutions {λ j , y λ j } j ∈Z + 0 and that its eigenvalues are simple, bounded from below and accumulate only at infinity as well as that its eigenfunctions oscillate as functions of the eigenvalues [13], namely λ j < λ j +1 , lim j →+∞ (λ j /j 2 ) = (π/(b − a)) 2 , y λ j has exactly j zeros in (a, b).
(3) Moreover, a representation of the (λ, y λ ) is given by the solution of the initial value problem with initial condition which characterises the eigenvalues as roots and the eigenfunctions via the transition matrix To compute (λ, y λ ) one can then for instance set up a mesh with m ∈ Z + intervals c 0 := a < c 1 < · · · < c m := b, h k := c k+1 − c k , h min := min and approximate (λ, c k ) → Y λ (c k ) by integral series techniques, which require truncation of the series and discretization of the integrals. The Fer streamers approximation to (λ, t) → Y λ (t) in the truncation of the integral series in the lead paper [14] and in the discretization of the multivariate integrals in this paper, is virtually unique in the literature because it is based on the uniform regime h max → 0 + , uniformly w.r.t. λ ∈ q max − h −2 max , +∞ .
The only partial exception known to the author is the use of Magnus expansions in [11] which yields a method with global order four based on the bounded uniform regime h → 0 + , uniformly w.r.t. λ ∈ −h −2 , h −2 .
The theory based on the uniform regime in Eq. 8 is very different from the theory based on the asymptotic regimes common throughout the literature λ fixed and h → 0 + , (9) h fixed and λ → +∞, e.g., in the piecewise perturbation methods [7], in the right-correction Magnus series [2] and in the modified Magnus methods [9]. In addition, one can consider another regime. Namely, there are well known asymptotic formulas for λ j /j 2 in j → +∞, of which the leading term is given in Eq. 3, which can attain high order for 'large' eigenvalues [13].
The main difference of the uniform regime in Eq. 8 over the asymptotic regimes in Eqs. 9-10 is that the former yields approximations which hold uniformly in (λ, h), while the latter only controls errors in one variable at a time in either h (with fixed λ) or λ (with fixed h).
This has important consequences for the size of the constants in the big O notation in the approximation estimates as well as their sampling requirements and their ranges of validity. These manifest in four ways: Firstly, Eqs. 9-10 often lead to truncation estimates with 'large' constants in the big O notation for 'intermediary' eigenvalues. Secondly, Eqs. 9-10 usually result in more function evaluations for 'intermediary' eigenvalues. Thirdly, Eq. 9 can lead to quadrature estimates with 'large' constants in the big O notation, which are applicable only to 'small' eigenvalues. Fourthly, Eq. 10 is usually used with highly oscillatory quadrature, which is applicable only to 'large' eigenvalues and non-resonant integrals.
To illustrate the first point above, note the increase in global order from 6 to 8, legitimate under Eq. 9, in the truncation of the integral series in [2, Theorem 2] and [9,Theorem 4.2], is built upon Taylor expansions of oscillatory integrals. Thus, the O h 10 terms in [1, p. 35] and [9, p. 759] grow with λ. As a consequence, the constants in the truncation estimates grow with λ and are therefore applicable only to 'small' eigenvalues.
As for the second point above, the localized increase in function evaluations outside of Eqs. 9-10 can be traced to the quadrature estimates in [4] for the integral (which appears in approximations for global order greater than or equal to four) which vary for non-oscillatory | λ,h | 1, highly-oscillatory λ,h 1 and intermediate | λ,h | ≈ 1 regimes, where the first two arise from Eqs. 9-10, and the last arises from 'intermediary' (λ, h). With p quadrature points, quadrature estimates read [4]: . The ratio (p + 1)/(2p − 1) then quantifies a 50% increase in sampling required for 'intermediary' eigenvalues.
Ref. [4] also clarifies the third point above since it is precisely to account for 'large' constants in the big O notation that the aforementioned estimates for | λ,h | ≈ 1 are worse than the ones for | λ,h | 1. For the fourth point above, we note that quadrature estimates that rely on Eq. 10 often arise from asymptotic expansions in λ,h 1 [4]. However, since for fixed λ,h , smaller h leads to larger λ, these estimates are only valid for 'large' eigenvalues. In addition, quadrature estimates under λ,h 1 are often valid only in the absence of critical points and subject to a non-resonance condition [6]. In the context of the integral series in [2,7,9,14] and this paper, these considerations are key since the non-resonance condition is not satisfied in the integral (which appears in higher order approximations) Working with Eq. 8 rather than Eqs. 9-10, Fer streamers bypass such issues and control the truncation and discretization of Sturm-Liouville problems with error bounds which hold with large step sizes equally well for all 'small', 'intermediary' and 'large' eigenvalues.

Geometric integration
The Lie-group and Lie-algebra play a natural role in the qualitative approximation of the solution of Eqs. 4-5 given that We note that the Fer streamers approximations and the methods in [2,9] preserve this geometric property, unlike the methods in [7].

Computational complexity
We now discuss how the computational complexity of Fer streamers compares with other geometric integration techniques.

Number of steps in each numerical mesh
The truncation in [14] and the discretization in the current paper are based on the following: Equation 12 prevents 'large' constants in the big O notation in the error estimates. In detail, if λ < q min then the argument of certain hyperbolic cosines and sines is positive. If left unchecked, the argument becomes unbounded and the hyperbolic functions grow exponentially with the size of the argument. Equation 12 guarantees that the positive argument and the hyperbolic functions are bounded by a small constant. Since [2, p. 423] assumes that λ q max , this issue does not arise. Ledoux et al. [9] does not assume λ ≥ q min , and disregards this issue. As discussed in [14], there are Sturm-Liouville problems where Eq. 12 is automatically satisfied since there do not exist eigenvalues smaller than q min . Equation 11 enables an unhindered transition of the error estimates between the uniform regime in Eq. 8. In addition, it quantifies the impact of the magnitude of the potential to the Fer streamers approach to Sturm-Liouville problems. The fact that the scale of the potential influences the step size is noted, but not quantified, in [2, p. 416] and [9, p. 761].
We note that Eqs. 11-12 use the knowledge of a lower bound to the minimum of the potential and an upper bound to its maximum. In [2] the knowledge of q max is required since it focuses on λ q max . Ledoux et al. [9] does not use this information, since it does not control the positive argument of certain hyperbolic functions for λ < q min .
Equation 13 controls the non-uniformity of the mesh, which is intrinsically related to the size of the constants in the big O notation in the error estimates.

Error estimates and evaluations of the potential
The discretization in this paper with global order 4, 7, 10 with respect to Eq. 8 requires 3, 6, 9 (one at the right boundary and the rest in the interior) evaluations of the potential per mesh interval. Given the potential is continuous, this translates to 3m + 1, 6m + 1, 9m + 1 evaluations of the potential for a mesh with m intervals.
The discretization in [2, p. 422-429] with global order 4,8 with respect to the regime in Eq. 9 requires 2,4 (interior) potential evaluations per mesh interval, which corresponds to 2m, 4m potential evaluations, for a mesh with m intervals. It is also shown in [2,Eq. 48] that the global error in that work is bounded in the regime in Eq. 10.
The work by [9] extends the work by [2] from λ q max to λ ≤ q max and suggests different potential evaluations. For global order 4,8 it uses 3,5 (one at each boundary and the rest in the interior) potential evaluations per mesh interval. Since the potential is continuous, this results in 2m + 1, 4m + 1 potential evaluations, for a mesh with m intervals.
As discussed in Section 1.1, the 50% increase in function evaluations from Eqs. 9-10 to 8, is necessary for quadrature over 'intermediary' eigenvalues.

Amount of linear algebra
The discretization in this paper boils down to the quadrature of integrals of the form where t → Z λ (ht) possesses a plethora of behaviour that varies with (λ, h). The quadrature schemes in this paper are based on uniform approximations of Z λ (ht) in (λ, t) byZ λ,h (t) with the property that Eqs. 14-16, with Z λ (h·) replaced byZ λ,h (·), can be integrated exactly. Thus, the amount of linear algebra in the discretization schemes can be quantified by the number of terms in each integrand in Eqs. 14-16 withZ λ,h (·) instead of Z λ (h·), which grows exponentially with i) base equal to the product between the number of summands in each representation of Z λ (ht) times the number of quadrature points, and, ii) exponent equal to the number of commutators in each integrand plus one. Fortunately, the exponential growth of the number of terms in each integrand is heavily attenuated in the quadrature in this paper, since it requires less sampling points for the higher dimensional integrals than for the lower dimensional integrals, which represents a significant saving in linear algebra. In detail, in the sense of Eq. 8, global order: -four requiresZ λ,h (t) with 3 sampling points for Eq. 14, -seven requiresZ λ,h (t) with 6 sampling points for Eq. 14 and 3 for Eq. 15, -ten requiresZ λ,h (t) with 9 sampling points for Eq. 14, 6 for Eq. 15 and 3 for Eq. 16.
In particular, [2,9] do not enjoy the heavy attenuation of the exponential growth of the number of terms in each integrand described above for the quadrature schemes in this paper.

Sturm-Liouville problems via Fer streamers and truncation error estimates
The Fer streamers approach to Sturm-Liouville problems introduced in [14] is built around the Fer expansions representation of the solution of Eqs. 4-5 given by the integral series in Theorem 1 below, which relate to the eigensystem of Eqs. 1-2 through 6.

Theorem 1 ([3]) The solution of Eqs. 4-5 is given by the Fer expansions flow
In practice, integral series such as Magnus, Fer or Neumann expansions, require truncation of the series and discretization of the integrals, which with Sturm-Liouville problems are challenging since the solution of Eqs. 4-5 is often exponentially large or highly oscillatory.
One of the contributions in [14] is to provide a reinterpretation of Fer expansions, which bypasses such issues, and opens the door to approximate (λ, t) → Y λ (t) uniformly and to high order in Eq. 8, equally well throughout large, oscillatory and in-between cases.
The first insight in [14] is to use Lie-group/Lie-algebra techniques to rewrite the infinite series in Definition 1, which appear in Fer expansions and are hard to control, by amenable closed-form expressions, named 'Fer streamers', given below in Theorem 2.
Definition 2 Let X ∈ sl(2, R) and z ∈ C, and define

Theorem 2 ([14]) If l ∈ Z + and t ∈ [c k , c k+1 ], then the infinite series in Definition 1 for the Fer expansions in Theorem 1 are given in closed-form by the 'Fer streamers'
Moreover, the first Fer streamer takes the form The second insight in [14] is then to partition the λ-axis into and use Fer streamers to control the quantities in Theorem 1 via the bounds with respect to Eq. 8 in Theorem 3 below, which make possible to derive the truncation estimates with respect to Eq. 8 in Theorem 4 below, which forms the main result in [14].
In particular, Theorem 4 reduces the problem of approximating the infinite product of exponentials in the exact flow F λ (c k , c k+1 ) in Theorem 1 to approximating the finite product of exponentials in the truncated flowF λ, r (c k , c k+1 ) in Definition 3.
In addition, the fact that X → exp(X) has the simple form in Definition 1, further reduces the problem to approximating the finite number of exponents
In the present section, we will develop specialized quadrature that approximates Eq. 22 to high-order uniformly for all eigenvalues. This, in turn, will give rise to the discretized flowF λ,r (c k , c k+1 ) and discretized solutionỸ λ,r (c k+1 ) in Definition 10 below, which will themselves be computed exactly. Towards the end, we will prove in Theorem 12 below that these approximate the truncated flowF λ, r (c k , c k+1 ) and truncated solutionỸ λ, r (c k+1 ) to global order 4, 7, 10 for r = 1, log(3)/ log(2), 2, respectively.
To this end, we begin with the observation that each exponent in Eq. 22 can be approximated by multivariate integrals over polytopes, each of which is non-trivial, apart from the first. Indeed, the first term amounts to the computation of which can be carried out without concern, while the second term can be written as and the third term can be controlled by in (19) and (20), in (19) and (20), where the first equality is due to Definition 1 and Theorem 2, and the last two follow from Definition 2 and Theorem 3. Hence, for global order less than or equal to ten, the second and third terms boil down to quadrature of In detail, Fer streamers require: -for global order 4, the quadrature of Eq. 23 to local order 5, -for global order 7, the quadrature of Eqs. 23-24 to local order 8, -for global order 10, the quadrature of Eqs. 23-25 to local order 11.
We develop quadrature of Eqs. 23-25 via representations of their integrands that exploit their: -magnitude to reduce quadrature points, function evaluations and linear algebra, and, -behaviour to decrease quadrature error without using derivatives of the potential.
These quadrature schemes and representations vary across the three subsets in Eqs. 19-21 where B λ,1 (c k , t) exhibits different behaviour and magnitude.

Exposing the behaviour of the integrands
We start by making the behaviour of B λ,1 (c k , t) explicit in Theorems 5-7 below, which we exploit thereafter.

Definition 4 Let
To provide intuition before plunging into technicalities, the reader should be aware that the guiding principle that leads to the representations in Theorems 5-7 below is to rewrite the representation of B λ,1 (c k , t) in Eqs. 17-18 in terms of trigonometric functions with the argument ω λ,1 (c k , t). To this end, recall Eqs. 17-18 and invoke Definitions 2 and 4 to rewrite Since r λ,1 (c k , t) is close to 1 and λ,1 (c k , t) is close to 0 uniformly for every eigenvalue, the behaviour of B λ,1 (c k , t) will be encapsulated in terms of trigonometric functions with argument ω λ,1 (c k , t), provided some care is taken to make every singularity removable. As will become clear, this serves to make the behaviour and magnitude of B λ,1 (c k , t) explicit, which, in turn, serves to reduce the function evaluations and linear algebra in the quadrature and to decrease the quadrature error without using derivatives of the potential.
To make this guiding principle precise the next definition introduces the non-oscillatory parts f λ,1 (c k , t), ι λ,1 (c k , t) and g λ,1 (c k , t), which appear below in Theorems 5-7. Although important, it is technical in nature and the reader is encouraged to glance over it and return to it as required.
With Definition 5 in hand, we encapsulate the behaviour of B λ,1 (c k , t) in the following three theorems concisely via the Hadamard/entrywise product denoted by . Eq. 19 then   π(B λ,1 (c k , t)

Theorem 5 If λ lies in
Furthermore, the derivatives f Furthermore, the derivatives ι (j ) λ,1 (c k , ·) can be bounded independently of λ.
Proof The representation follows from Eqs. 17-18 and Definitions 2 and 5 together with Eqs. 28-30. The terms are arranged in order to make explicit and render every singularity removable. Finally, g λ,1 (c k , t) can be bounded independently of λ because the derivatives of Eq. 27 can be bounded independently of λ.

Definition 6 below introduces B fine λ,1 (c k , c k + h k t) as a means to decompose the fine and coarse parts of B λ,1 (c k , c k + h k t). This fine and coarse decomposition is made precise with Corollaries 1-3 below which show that
uniformly over the entire eigenvalue range. It is then in the next subsection, in Theorems 8-10, that this fine and coarse decomposition is shown to bear fruit in the form of reduced requirements for quadrature.

Definition 6 Let B fine
λ,1 (c k , c k + h k t) be the unique element in sl(2, R) such that π(B λ,1 (c k , c k + h k t)) (19) and (20), The next three Corollaries follow immediately from Theorems 5-7 and Definition 6.

Reduced requirements for quadrature
Theorems 8-10 below highlight the synergy between the Lie bracket and the representations in Theorems 5-7, Definition 6 and Corollaries 1-3: they act together to decrease the magnitude of each multivariate integral, making it smaller than expected! In particular, as discussed in the subsequent text, they illustrate that it is possible to use these representations to develop a quadrature which exploits the magnitude of each integrand in order to reduce the number of function evaluations and amount of linear algebra.  ) and (20), Proof Follows by straightforward computation from Definition 6.  (19) and (20), Proof Follows by straightforward computation from Definition 6.  (19) and (20), Proof Follows by straightforward computation from Definition 6.
As a result of Theorems 8-10, quadrature of Eqs. 23-25 should be replaced with that of since this results in fewer function evaluations and amount of linear algebra, given that: -for global order 4, as a consequence of Theorem 8, quadrature of Eq. 23 to local order 5 is equivalent to quadrature of Eq. 34 to local order 3, -for global order 7, according to Theorem 8, quadrature of Eq. 23 to local order 8 is equivalent to quadrature of Eq. 34 to local order 6, and, following Theorem 9, quadrature of Eq. 24 to local order 8 is equivalent to quadrature of Eq. 35 to local order 3, -for global order 10, Theorem 8 guarantees quadrature of Eq. 23 to local order 11 is equivalent to quadrature of Eq. 34 to local order 9, Theorem 9 ensures quadrature of Eq. 24 to local order 11 is equivalent to quadrature of Eq. 35 to local order 6, and, Theorem 10 ensures quadrature of Eq. 25 to local order 11 is equivalent to quadrature of Eq. 36 to local order 3.

Interpolation
We develop quadrature for Eq. 34-36 by polynomial interpolation of the nonoscillatory parts of B fine λ,1 (c k , c k + h k t), given by f λ,1 (c k , c k + h k t), ι λ,1 (c k , c k + h k t) and g λ,1 (c k , c k + h k t), as exposed in Corollaries 1-3. This gives rise to the approximationsB fine λ,1,T j (c k , c k + h k t) in: in each of Eqs. 19, 20 and 21 by the right hand side of, respectively: -Eq. 31 with t → [f λ,1 (c k , c k + h k t)] i,1 replaced by polynomial interpolation at T j , -Eq. 32 with t → [ι λ,1 (c k , c k + h k t)] i,1 replaced by polynomial interpolation at T j , -Eq. 33 with t → [g λ,1 (c k , c k + h k t)] i,1 replaced by polynomial interpolation at T j .
The question then becomes which (j, T j ) should be used for quadrature of Eq. 34-36 to achieve prescribed local order. A question we address in the next three subsubsections:
In particular, fewer interpolation points are needed for higher dimensional integrals than for lower dimensional integrals, which represents a significant saving in linear algebra.

Interpolation points that reduce quadrature error without derivatives of the potential
As in the analysis in [5,6,10] for highly oscillatory Fourier-type integrands, the mildly to highly oscillatory behaviour in Eq. 21 made explicit in Corollary 3 also suggests polynomial interpolation at the endpoints since this yields smaller quadrature error. This would lead to the evaluation of g λ,1 (c k , c k ), which depends on q c + k : something that would be best to avoid since the derivative of the potential might not be available. Fortunately, since in Corollary 3 there is a 't' term in front of every ' g λ,1 (c k , c k + h k t) i,1 ' term, this is automatically achieved at the left boundary point. Hence, there is no need to interpolate at the left boundary point. As for t ∈ (0, 1], g λ,1 (c k , c k + h k t) does not depend on the derivative of the potential and should be interpolated at the right boundary point. Given this, choose τ 1 : = 0 and τ j := 1, in which case, t → tp f λ,1 (c k ,c k +h k ·),j −1 (t) is the unique j degree interpolation polynomial which interpolates t → tf λ,1 (c k , c k + h k t) at the j + 1 points 0, τ 1 , · · · , τ j −1 , 1 .

Data
As discussed in Sections 3.4.1-3.4.2, it is beneficial to choose T j := S j −1 ∪ {1} where S j −1 :⊆ (0, 1) has j − 1 distinct points where j = 3, 6, 9 vary with Eqs. 34-36 and local order. To reduce potential evaluations, take S 2 :⊆ S 5 :⊆ S 8 , in which case, one lets S := S l−2 for global order l, in which case, the polynomial interpolation in this section requires the following data If the antiderivative of the potential is not available, then it is possible to approximate, up to local order, the antiderivative data and the exact integration of the result.

Discretization error estimates
Definition 9 and Theorem 11 below make explicit the quadrature error in Section 3.5, while Definition 10 and Theorem 12 below clarify the manner in which the quadrature error, which lives in the Lie algebra sl(2, R), affects the quantities in the Lie group SL(2, R).
Proof Follows from Theorems 8-10 and the discussion in Section 3.4.

Definition 10
Let r ∈ {1, log(3)/ log(2), 2}, and define the discretized flow, the discretized solution, the discretization local error, and the discretization global error by, respectively:  (19) and (20), (19) and (20), Regarding local error, for r = 1 the local discretization error can be written as where the first and second equalities are due to Definitions 3 and 10, the third equality is due to Eq. 38 and the fourth equality is due to Eq. 40. More generally, for arbitrary r, one can show which, with Theorems 3 and 11, yields the desired estimate. Regarding global error, for arbitrary r, the global discretization error obeys the recursion relation with initial condition and general rule = log e L disc. λ,r (c k ,c k+1 ) exp Ad exp(D λ,0 (c k ,c k+1 )) G disc. λ,r (c k ) + h.o.t.
where the first four equalities are due to Definitions 3 and 10, the fifth to Eq. 39, the sixth to Eq. 40, and the last to Eq. 38. The global discretization error expressions in Eqs. 41-42 then result in which, with Assumption 1, Theorem 3 and Theorem 11, result in the desired estimate.

Total error estimates
Having controlled in Section 2 the truncation error incurred from approximating F λ (c k , c k+1 ) and Y λ (c k+1 ) byF λ, r (c k , c k+1 ) andỸ λ, r (c k+1 ) for r ∈ Z + , and controlled in Section 3 the discretization error sustained from approximat-ingF λ, r (c k , c k+1 ) andỸ λ, r (c k+1 ) byF λ,r (c k , c k+1 ) andỸ λ,r (c k+1 ) for r ∈ The previous theorem links the truncation estimates in [14] with the discretization estimates in this paper in that their sum controls the total error in the Fer streamers approach. As can be seen from Theorems 4 and 12, the bounds on the discretization errors λ,r (c k , c k+1 ), G disc. λ,r (c k+1 ), are larger than the bounds on the truncation errors L trun.

Conclusions
It has been shown that to preserve the advantageous properties of the Fer streamers approach to Sturm-Liouville problems in [14] under discretization, while simultaneously minimizing the computational complexity by reducing function evaluations and amount of linear algebra in the discretization, quadrature requires not the original representation of Fer streamers, but rather relies on a characterization designed to expose their magnitude and behaviour. Tight total error estimates, uniform for every eigenvalue, have also been established that quantify the interplay between the truncation and discretization in the approach by Fer streamers which have shown that the discretization errors outweigh the truncation errors.
The principal advantage of the Fer streamers approach to Sturm-Liouville problems is that its truncation and discretization error estimates hold uniformly for all eigenvalues, which guarantee large step sizes over the eigenvalue range. This is especially significant given that the error estimates in alternative techniques apply only to 'small' or 'large' eigenvalues (c.f. Section 1.1). Compared with the alternative geometric integration techniques in the right-correction Magnus series [2] and in the modified Magnus methods [9], which do not possess error estimates uniform over the entire eigenvalue range, the Fer streamers approach presents an interesting trade-off in computational complexity: although it requires a 50% increase in function evaluations for each univariate integral in order to control all 'small', 'intermediary' and 'large' eigenvalues, it also enjoys a significant decrease in linear algebra for each multivariate integral (see Section 1.

3).
These results open new questions which will be the subject of future work and include: -Free Lie algebras, and, -Step size strategies.

Free Lie algebras
As discussed in this paper, in order to approximate Y λ (c k+1 ) byỸ λ,r (c k+1 ), it is necessary to computeF λ,r (c k , c k+1 ), which itself relies on the computation of I fine λ,i,T j (c k , c k+1 ). Given that each integralĨ fine λ,i,T j (c k , c k+1 ) can be integrated exactly by first expanding the individual terms of the commutators of the finite sums iñ B fine λ,1,T j (c k , c k + h k t), then integrating exactly every summand, and finally adding up all the integrated terms, a useful question is how to perform such operations while incurring the smallest possible volume of linear algebra and, by extension, running time. Since this volume of linear algebra relates to the number of commutators that result from integrating the individual terms, there exist three mechanisms that can be used to decrease this volume of linear algebra, which are: -Firstly, free Lie algebra techniques and Hall basis, which lead to fewer commutators via a systematic use of commutator identities such as: skew symmetry, Jacobi's identity, etc, -Secondly, when collected in a Hall basis, certain linear combinations between different integrands are then identically zero, and, -Thirdly, when collected in a Hall basis, certain linear combinations between different integrands integrate exactly to zero.
Hence, these efficiency mechanisms will be explored to optimise running time.
For these problems, Fig. 1 depicts the absolute and relative errors between reference eigenvalues and eigenvalues computed via Fer streamers with the procedure described above. Figure 1 highlights two issues related to the mesh strategy described at the start of this subsection. On the one hand, it yields step size h = 0.03 for Eq. 43 that gives relative error around 10 −16 for λ ≥ 0.5 × 10 5 which may be more accuracy than necessary, and thus may consume more resources (e.g., running time or pointwise evaluations of the potential) than necessary. On the other hand, for Eq. 44 the same strategy yields relative error around 10 −8 for every λ with step size h = 0.57.
These examples illustrate the fact that the behaviour of the solutions of Sturm-Liouville problems depend on the local behaviour of the potential q. With this in mind, one can: -Relax Assumption 1 which currently uses global information about q, namely, its minimum q min and maximum q max , to instead use local information on each [c k , c k+1 ] based on the pointwise evaluations of q at the interpolation points of Section 3.4.3, and, -Explore non-equidistant meshes, and investigate automatic mesh selection algorithms, along the lines of [8].