The Numerov process over a non-uniform grid

The Numerov process is a solution method applicable to some classes of differential equations, that provides an error term of the fifth order in the grid size with a computational cost comparable to that of the finite-difference scheme. In the original formulation of the method, a uniform grid size is required; the paper shows a procedure for extending its applicability to a non-uniform grid in one dimension. The effectiveness of the procedure is tested on a model problem, and comparisons with other methods are carried out. Finally, it is shown how to extend the applicability of the method to a larger class of equations; among these, the mathematical model of semiconductor devices is important in view of its applications to the integrated-circuit technology.


Introduction
The Numerov Process (NP [1], cited in [2]), is applicable to the solution of some classes of differential equations; its error term is O(h 5 ) , with h the grid size, much superior to that of the finite-difference scheme, with a comparable computational cost. In the original formulation of NP, a uniform grid is necessary, which is a drawback for many practical applications; more recently, a method for extending the scheme to a variable step in one dimension has been shown [3]. This paper shows an approach, different from that of [3], for extending NP to a non-uniform grid in one dimension; the approach preserves the accuracy of the original NP by combining the latter with a numerical integration having the same order of error.
The interest in highly-accurate schemes derives, among others, from the need of solving the transport model of semiconductor devices, which is fundamental for the integrated-circuit technology. Unfortunately, the form of the semiconductor-device model as it stands is such that NP is not applicable. In the paper, a transformation is shown that makes the semiconductor model tractable with the onedimensional, variable-step version of NP.
The paper is organized as follows: the extension of NP to a one-dimensional, non-uniform grid is shown in Sect. 2, and the properties of the matrix derived from it are outlined in Sect. 3; the application to a model problem, along with the comparison with the finite-difference scheme, is shown in Sect. 4 using different grids. A comparison with the approach of [3] is carried out in Sect. 5. A more general class of differential equations is considered in Sect. 6; it is shown that, by a suitable transformation of the unknown function, also this class of equations becomes amenable to the application of NP. The example is important because it includes the matematical model of semiconductor devices. Finally, the conclusions are drawn in Sect. 7.

The Numerov process over a non-uniform grid
We start by considering the one-dimensional, second-order equation of the form (1) −u �� = F(u, x) , 1 3 where the dependence of F on u may be nonlinear. The standard Numerov Process (NP) requires a uniform grid (the derivation of NP is shown in Appendix I); uniformity is necessary because NP is based on series expansions, and the uniform spacing of the grid points makes the odd powers of the right and left expansions to cancel each other. After eliminating the odd powers, one determines the even powers in terms of u by recursively applying (1); this is also the reason why it is necessary to consider a form like (1) of the equation, where the first derivative of the unknown function does not appear.
To extend NP to a non-uniform grid one must preserve the possibility to cancel, albeit only locally, the odd powers of the expansions. An approach that fulfills this requirement is shown here, which combines a formally exact solution over a non-uniform grid with the improved precision of NP. Consider a one-dimensional, non-uniform grid with N internal nodes, namely x 0 < x 1 < ⋯ < x N < x N+1 , and let h i+1 = x i+1 − x i be the size of one of its elements. With reference to (1), for any Repeating on the next element x i ≤ x ≤ x i+1 , and integrating again, yields, after some calculations whose details are given in Appendix II, Fixing to x i the second integration limit of (2), and renaming from to x the integration variable transforms (2) into Similarly, fixing to x i+1 the second integration limit of (3), and renaming from to x the integration variable transforms The problem is thus shifted to the calculation of the two integrals Due to the form of F, integrals (6) depend on u, but not on u ′ . If G i and H i+1 are approximated with formulas where only the nodal values of u appear, namely, , the resulting expressions do not introduce extra unknowns with respect to those that explicitly appear in (4) and (5). As shown later, this goal can be accomplished in a manner that exploits the precision of NP. Assuming that integrals (6) have been calculated, the way in which the nodal values of u are found from (4), (5) depends on the boundary conditions. If the boundary conditions u 0 , u ′ 0 are given, the starting point is calculating u 1 from (5) after letting i = 0 there; to calculate u 2 one needs u ′ 1 , which is obtained by letting i = 1 in (4); then, one proceeds recursively until the Nth node is reached.
If, instead, the boundary conditions u 0 , u N+1 are given (which is the typical case in, e.g., the numerical analysis of semiconductor devices), one must eliminate the first derivatives from the system (4), (5). For this, one rearranges (4) to find shifting the index in the latter expression then yields u � i−1 = (u i − u i−1 + H i )∕h i which, combined with the former, provides then, one manipulates the term G i − H i ∕h i in order to give (7) the more symmetric form So far, no approximation has been introduced; to calculate the integrals in (8) one now applies Simpson's rule with c the midpoint of [a, b]. The error term of (9) is 5 (4) ] , with (4) the 4th derivative of calculated at some point of the integration interval [4].
Using index i − 1∕2 to denote the quantities at the midpoint of h i , index i + 1∕2 to denote those at the midpoint of h i+1 , and applying (9), transforms (8) into One notes that the left-hand side of (10) corresponds to the discretization of the second derivative obtained from a parabolic interpolation over a non-uniform grid; in turn, the central term of the right-hand side (without the 1/3 factor) is the discretized form of the right-hand side of −u �� = F when the finite-difference method is used (compare with (22)). The improvement inherent in (10) derives from Simpson's rule whose error term, as mentioned above, is O(h 5 i ) in the interval on the left of x i and O(h 5 i+1 ) in the interval on the right of it. The issue is to calculate the values F i−1∕2 and F i+1∕2 , that belong to the midpoint of each interval.

The system matrix
From the discussion carried out at the end of Appendix I, one may assume that the equation in hand has been linearized and that the solution is being carried out by iterations; therefore, within the single iteration the dependence of F on u is linear, and F i+1∕2 = c i+1∕2 u i+1∕2 + s i+1∕2 ; Eq. (10) then becomes In (11), the values of c and s at the nodes and at the midpoints of the elements are known, because the two functions are either prescribed or, like in the case considered here, are taken from the previous step of an iterative procedure. In contrast, u i−1∕2 and u i+1∕2 are not known; however, u is the solution of (1), for which the NP interpolation (41) holds because the three nodes considered are equally spaced. For u i−1∕2 , replicating (41) over the three nodes in hand yields Replacing the expressions of u i−1∕2 and u i+1∕2 into (11) and rearranging, yields with In conclusion, by combining the Simpson rule with an interpolation of the NP type, integrals (8) have been calculated in terms of the nodal values u i−1 , u i , u i+1 only, and the accuracy of NP has been kept. No constraint has been imposed on the elements; as a consequence, the scheme is applicable to a general, non-uniform grid. The outcome of the whole procedure is the N × N algebraic system (15). One notes that the system matrix of (15) is tridiagonal like in the finite-difference scheme; however, it is not symmetric because i−1 ≠ i due to the fact that a i−1 i ≠ a i i ; the asymmetry is not to be ascribed to the non-uniformity of the grid, nor to the use of NP, but to the presence of term c u in the equation to be solved. In fact, even if a uniform grid is used, the matrix resulting from it is still asymmetric (compared with (41)); conversely, if it happens that c = const , then (13) renders a i−1 i = a i i and makes the matrix symmetric. Another observation is that, if one takes the limit of a uniform grid, h i → h , interpolation (10) becomes Apart from the trivial case F = const , limit (19) is different from the interpolation found by directly applying NP to a (14) uniform grid, that is, (35  (35), one changes the coefficient of the right-hand side of (19) from h 2 ∕3 to h 2 ∕12 ; in this way, the weights of the F terms become (4,4,4), to be compared with (1, 10, 1) of (35). In both cases, the weights total to 12; however, the lateral weights in (35) are smaller because the lateral nodes of (35) are more distant from the central node x i than the lateral nodes of (19).

Model problem
In this section, we consider a model problem useful for an experimental validation of the procedure outlined in Sects. 2 and 3; we take the interval 0 ≤ x ≤ 1 and assume that the solution of −u �� = c u + s is with k, p two positive integers, so that the boundary conditions are u(0) = u(1) = 0 . One finds Function (20) is such that the frequency and amplitude of the oscillations increase from right to left; the case corresponding to k = 2 and p = 5 is considered (function (20) is shown as the green line in Figs. 1 and 2). Equation (20), with coefficients given by (21), is solved on different grids using the generalized NP method of Sects. 2 and 3, or the standard finite-difference scheme (FD); the latter yields Each outcome is compared with the true solution (20); more specifically, as error indicators we adopt .
where u FD i or u NP i is the value of the solution obtained at the ith node using the FD or NP scheme, and u i is the value of (20) at the same node.
The first comparison has been carried out by randomly generating the nodal positions within the integration domain; other comparisons have been carried out starting from a uniform grid, whose nodes have subsequently been shifted to obtain a prescribed density of nodes over the integration domain. Specifically, the shift was obtained by defining the new nodal positions y i , i = 1, … , N with the transformation where p is the positive integer appearing in the definition (20). Observing that the nodal density of the uniform grid is Q = 1∕(N + 1) , and that the density g(y) of the shifted grid fulfills the relation Q dx = g(y) dy , one finds In particular it is g(0) = (p + 1) g(1) , namely, the ratio g(0)/g(1) is the same as The error indicators (23) are shown in Table 1 for different numbers of grid nodes N; specifically, cols. 1, 2 of the table show the results obtained from the randomlygenerated grids, cols. 3, 4 those from the uniform grids, and cols. 5, 6 those from the left-shifted grids (although considering uniform grids is redundant with respect to the scope of this work, the results are shown for comparison with those of the other cases). In each set of simulations, the number of nodes has progressively been reduced until the error indicator (23) reached a value of the same order as the oscillation amplitude of u; after this threshold was reached, the number of nodes was not reduced any more and the simulation was discontinued, which explains the horizontal lines in Table 1.
Considering the general trend of the error indicator, the tables show that, as expected, the error improves from the case of the randomly-generated grids, to that of the uniform grids and, finally, to that of the left-shifted grids. As for the comparison between the two solution methods, it is found in all cases that the non-uniform NP is better by orders of magnitude than the corresponding FD.

Comparison with other approaches
An extension of NP to a non-uniform, one-dimensional grid has been proposed by Vigo-Aguiar and Ramos as a special case of the variable k-step Störmer-Cowell method [3]. One of the outcomes of [3] (indicated here as VR from the authors' initials) is determining a strategy for selecting the length of element h i+1 , given h 1 , … , h i , when solving (1) when the initial conditions u(x 0 ) = u 0 , u � (x 0 ) = u � 0 are prescribed; as shown by examples provided in [3], VR is applicable also to boundary-value problems. In the discretization scheme of VR, the left-hand side can be reduced to the form of (10); the right-hand side uses weights in which suitable combinations of the elements appear. Another work showing an extension of NP to a non-uniform, one-dimensional grid is [5], whose scheme is similar to that of VR: the left-hand side is still the same as in (10), whereas the weights at the right-hand side use the single ratio h i+1 ∕h i ; this yields an O(h 3 ) error term [6]. In the following, the comparison is carried out using again the model problem −u �� = c u + s with coefficients given by (20, 21); the method proposed here is compared with VR and with the boundary-value problem solver of MATLAB® [7] (the latter will be indicated with ML below). In contrast to VR, the scheme of ML considers a system of first-order equations, that are solved by a collocation method with a C 1 piecewise-cubic polynomial; the latter collocates at the ends of each element and at the midpoint, yielding an O(h 4 ) error term. It is also worth noting that ML constructs - a non-uniform grid by adding nodes until the specified tolerance is reached; this explains why the numbers of nodes of Table 2 are different from those of Table 1.
Remembering the description of the present method given in Sects. 2, 3, the equation to be solved is recast in integral form and the integral over each element is calculated by the Simpson rule; for the latter, the value of u at the element's midpoint is needed, which is obtained from the NP interpolation (41). It may be argued that the introduction of such a "ghost" point in the middle of each element is equivalent to doubling the number of nodes, which is of advantage when the error indicator (23) is used. For this reason, another set of simulations have been carried out with VR using a double number of nodes with respect to that used with NP; the results of such simulations (which are more expensive in terms of computing time, see Table 3) are labeled VR2. Instead, the doubling of the number of nodes was not used with ML because, as mentioned above, this method uses the elements' midpoint as well. Table 2 compares the error indicators of type (23) obtained with different solution methods. The table has been obtained by solving (1) using ML with six different values of its internal relative-tolerance parameter. This provided six non-uniform grids, whose numbers of nodes are given in column N; the corresponding error indicators are shown in column ML . Then, using the same grids, the solutions and the corresponding error indicators have been calculated again using the VR method with N nodes ( VR ), the VR method with 2 N nodes ( VR2 ), and the non-uniform-grid NP method ( NP ). The solutions corresponding to the last two lines of Table 2 are shown in Figs. 1 and 2, respectively, for ML and NP, that exhibit the smallest error indicators.
A last set of runs have been carried out to compare FD , VR , and NP over randomly-generated grids (Fig. 3), and over initially-uniform grids shifted by applying the (24) scheme (Fig. 4). Finally, the simulation times of the different methods are compared in Table 3 using a uniform grid by way of example; the calculations have been carried out on a Notebook with an AMD A4-4355N cpu, 1.9 GHz, 4.0 GB RAM.  1.2 × 10 0 2.2 × 10 +2 6.0 × 10 +1 2.1 × 10 −1 76 5.1 × 10 0 3.3 × 10 +2 2.9 × 10 +2 1.7 × 10 0  where b and Q are prescribed. The NP method cannot be applied to (26) as it stands; however, the form of (26) is readily reduced to that of (1) by introducing the auxiliary unknown with x 0 an arbitrary point; in fact, w fulfills An example of equation of the form (26) is given by the transport model of the semiconductor devices, which is important in view of its applications to the integrated-circuit technology; when the transport of electrons is considered, b stands for the normalized electric field and g for the electron concentration. Examples of applications of the NP method to the semiconductor model, over a uniform grid, are given in [8,9].

Conclusions
A method for extending the Numerov process to a non-uniform grid in one dimension has been shown. The result is achieved by acting on each grid element separately: first, the equation to be solved is locally recast in integral form, and the integral over the element is expressed with the Simpson rule; then, the extra unknown allocated at the midpoint of the element is expressed through the Numerov interpolation in terms of the two nodal unknowns belonging to the same element. The error term of the Simpson rule is of the same order as that of the Numerov process, so the accuracy of the latter is preserved. The effectiveness of the scheme has been tested on a model problem, in which both amplitude and oscillation frequency of the solution increase when one of the ends of the integration domain is approached. The method compares favorably with other methods that extend the Numerov Process over a non-uniform grid, and with MATLAB®. In the last part of the paper, it is shown how a suitable transformation of the unknown makes it possible to apply the Numerov process to a larger class of equations, that includes the mathematical model of semiconductor devices.

Appendix I: the Numerov process
An interpolation scheme more refined than the parabolic one is the so-called Numerov Process [1] (cited in [2]), which applies to second-order equations of the form The essence of the method is that the derivatives of the odd order are canceled, while those of even order are recursively reconstructed in terms of u by means of (29). Considering a uniform grid of size h, the nodal value at u i+1 is obtained by means of a Taylor expansion starting from u i ; truncating the expansion to the 5th order yields and, similarly, Adding (31) to (30), Now, the same form of the interpolation is obtained if the series expansion is applied to u ′′ instead of u; one finds in this case Neglecting the 6th-order term in (33) and eliminating u (4) i between (32) and (33) yields This interpolation is obtained by keeping terms up to the 5th order; as a consequence, its precision is better than that of the parabolic interpolation. The additional computational cost is due to the right-hand side of (35), whose calculation requires one more multiplication and two more additions with respect to that of the right-hand side of standard (29) −u �� = F(u, x) . (30) finite-difference scheme; the cost of inverting the matrix of the algebraic system (35) is, instead, the same. As anticipated, the method eliminates (without approximations) the odd-order derivatives and exploits the form of the equation to eliminate the even-order ones. This does not mean that the calculation of the odd-order derivatives is prevented; in fact, subtracting (31) from (30) one obtains Neglecting the 5th-order term in (36) and repeating the calculation for u ′′ , Neglecting again the 5th-order term and eliminating u (3) i between (36) and (37), which is rearranged as The calculation of the derivative keeps the terms up to the 4th order. It is also important to note that (39) provides the derivative at the ith node, not within one or the other neighboring elements like in the parabolic approximation.
Still considering (29), it is convenient to consider separately the case where the dependence of F on u is linear, namely, where c, s are two given functions of x, from the case where the dependence is nonlinear. In the linear case (35) and (39) become, respectively, In the nonlinear case, which is typical of, e.g., the mathematical model of semiconductor devices, one must preliminarily linearize (29) starting from a tentative solution ū ; specifically, one lets u =ū + u , this transforming (29) into which must then be solved by iterations. When dealing with (43) one may still use (41) and (42) provided the following replacements are used: At convergence, u = 0 and −ū �� = F(ū, x).
Availability of data and material (data transparency): Data and material are available.

Conflicts of interest Not applicable.
Code availability: Code is available.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.