Spectral Rectangular Collocation Formula: An Approach for Solving Oscillatory Initial Value Problems and/or Boundary Value Problems in Ordinary Differential Equations

The idea of rectangularization of a N N × differentiation matrix through collocation was recently suggested to be more useful in discretization process, especially when the traditional row replacement approach fails. In this regard, we employ the state-of-the-art technique of rectangularization to discretize some oscillatory initial value problems (IVPs) and also extended the new approach to a non-linear boundary value problem (BVP). The numerical implementation was composed using some few lines of executable Python codes. Our findings are instructive and quite revealing and supported by evidence from our numerical experiments and simulations.


Introduction
The genesis of the word spectral is not absolutely clear but presumably, arises from the genuine use of Fourier sines and cosines as basis functions [1] particularly in connection with scale between two extreme points, namely the spectrum. Spectral methods should not be confused with finite element methods (FEMs) based on the fact that they are closely related. The spectral methods use basis functions which are non-zero over the whole region, while the use of basis functions which are non-zero only on small sub-regions is particular to FEMs. Also, spectral methods are known to be global methods, and seek one solution for the entire region, rather than finite difference or finite element methods which seek a solution only in a local area to the element or node N.
Spectral collocation methods (SCMs) mainly use Chebyshev properties for its computations, reason why they are sometimes referred to as Chebyshev spectral collocation or pseudospectral methods. These SCMs have been firmly established in the numerical solution of differential equations (DEs) see [2,3,4] and other references therein. SCMs are based on the representation of a real, continuous, well-behaved function on some interval [a, b] using fast and stable algorithms such as barycentric interpolation formula in [5] and finite Fourier transform in [2]. Compared with other formalizations of spectral methods such as the Galerkin approximations (see [6]), they offer convenience in the presence of variable coefficients and non-linearities [7] which make these algorithms more accurate, for a given computational power. SCMs have become increasingly popular for solving DEs due to the fact that they enjoy rapid convergence. See for example the solution of unsteady boundary layer flows that are described by systems of coupled non-linear DEs in [8]. The unknown solution of the DE is approximated by a global interpolant, such as polynomial of high degree. This global interpolant is then differentiated exactly, and the expansion coefficients are determined by requiring the equations to be satisfied at an appropriate number of collocation points [9]. Since interpolation, differentiation, and evaluation are all linear operations, the process of obtaining approximations to the values of the derivative of a function at the collocation points can be expressed as a matrix-vector multiplication [9].
Spectral discretization of DEs have been of great use in computational mathematics, engineering and other branches of science. This idea become more interesting with the emergence of the traditional row replacement strategy see ( [3,10,11,12,13]) but not quite efficient for solving DEs with boundary conditions other than -1 and 1. In this case, the choice of row to enforce the boundary conditions become ad hoc [7]. Due to this snag, [7] introduced an efficient means of rectangularizing the N N × differentiation matrix to enable the enforcement of the boundary conditions. In this paper, we discretize the problem posed in the sense of the work in [7]. However, the explicit form of the work in [7] can be found in [11]. Interestingly, emerging computational approach for reformulating all kinds of DEs to systems of first order equations as been established in [14].

Barycentric Lagrange Interpolation Formulae
The Lagrange formula is in most cases the method of choice for dealing with polynomial interpolants contrary to the general belief that it is not good for approximation because of the identity property. The point is that the Lagrange formula must be modified through the formulae of barycentric interpolation. Consider the Lagrange formula as 0 ( ) ( ) , with property 1, ( ) , 0,..., . 0, , As in [5], we define the barycentric weight as 1 , 0,..., , and the barycentric formula of the first kind as 0 ( ) ( ) .
Equation (5), like the other authors is referred to as the improved Lagrange formula of (1). The formula (5) now  formula depends on the order which increases the computational cost as N increases, this many ordering leads to numerical instability [5]. However, these does not imply that the Newton formula have no relevance again. There are many cases where Newton interpolation is preferable, see for example, the construction of the multi-step methods for the numerical solution of ordinary differential equations (ODEs) in [15,16,17,18,19]. In addition, the barycentric formula of the first kind (5) can be further modified into a more elegant method, which is usually used in practice (see [5]). The barycentric formula of the second kind is given by 0 0 ( ) .
The barycentric Lagrange interpolation formula (6) is the one with special and beautiful symmetry [5]. If weights i w (4) have common factor, they may be cancelled without any impact on ( ) P x [20]. Methods (6) and (5) (6) can be found in [5]. Fast reduction of generalized companion matrix pairs for the barycentric Lagrange interpolants have also been consider in the excellent works in [21,22]. The formulae defined in this paper are being collocated at specifically chosen set of points 0 ,... . , N x x As in [3], the Chebyshev point of the first kind given by and the Chebyshev point of the second kind is given by The interpolation based on these points (7,8) converges as . N → ∞ Although, (8) has been used in practice due to the fact that it takes values at the boundary. Recently, equation (7) is becoming of necessity because of the crucial role it plays in the rectangularization of the N N × matrix, see the excellent work in [7,10,23]. The stability of barycentric Lagrange interpolation formulae (5) and (6) have been investigated in [24,25,26]. Barycentric formula of the second kind (6) may diverge when interpolated globally but can still enjoy the spectral accuracy provided that the given data if outside the interval [-1, 1] is scale back to the interval [-1, 1]. This can be achieved by using the transformation:  (6) is that it generalizes to rational functions which has been a great tool for building Chebfun software package [25].

Spectral Rectangularization through Collocation
Spectral rectangularization approach as demonstrated recently in the work in [7] is proved to be a convenient means of solving linear and non-linear DEs with general boundary conditions. This new approach in [7] was used in an elegant manner to construct a re-sampling matrix through the use of the barycentric formula (6). This barycentric re-sampling is an efficient, fast, and stable method to convert a polynomial interpolant, as earlier indicated by a vector of its value at a set of nodes, to its representation as values on other set of nodes. The special formulae given in this section are in the sense of the formulae discussed in the preceding section. Given a set of points vector where the re-sampling matrix ,− ∈ ℝ − × is given by .
The vector  [7]. When 0 m > , the re-sampling outcomes in an unchanged loss of information as a new vector of function values is shorter than the original. This process [7] referred to as down-sampling effect but the boundary conditions are appended to the rectangular system to yield a square matrix again. From the rectangular point of view, an m th-order differential operator is naturally discretized by an ( ) N N m × − matrix, allowing the m boundary constraints to be appended to form an invertible ( ) ( ) N m N m − × − system [7]. Naturally, the barycentric for these points are known as Consider a differentiation matrix which the generalized form is given as Then the differentiation matrix N D is pre-multiply by the re-sampling matrix (12). In particular, the case of a homogeneous differential equation of the form The simple DE (14) can easily be discretized by setting it as a system The second order differential operator (14) (2) D is rectangularized by a 2 N N − × re-sampling matrix (12) on the left hand side and similarly for f on the right hand side to give room for enforcement of the boundary conditions which makes the system a square matrix again. Hence, the system can be solved easily.

Numerical Experiment
In this section, we employ the relation ,  (15) to obtain the N for the discretization process of each problem in this section. We have successfully introduced this transformation for obtaining the node in the course of the research as it can rarely be found in the literature. Consequently, this is to synergize the relationship between the step size and node during the implementation process.
The implementation was carried out using numpy, numpy.linalg, pylab while plottings were obtained using matplotlib.pyplot, which are all packages on Python. The DE (17) represents motion on a perturbed circular orbit in complex plane in which the point u(x) spirals slowly outward such its distance from the origin at any given time x is 2 2 ( ) (Re( ( ))) (Im( ( ))) .

Spectral Rectangularization Approach for Non-linear DE
The spectral rectangularization by collocation method for non-linear DEs is another area of exploit. In this non-linear case, the scalar boundary value problem (BVP) brings no significant complications for the approach. We consider an mth order non-linear BVP where h is the non-linear differential operator and i g is the non-linear differential functional. To achieve this task, we employ the Newton-iterative scheme for polynomial approximation as illustrated in [31]. At 1 m = in (19), we have the vector form as On perturbing the solution u, say δ + u and expanding ( ) We multiply the re-sampling matrix , N m P − by (24).
, , where ′ h is a N square matrix. Now, for the case of the system We can express ( ) where  is the Hadamard point-wise or element-wise product. We then differentiate ( ) h u + ϒ u (32) We are interested in the element in the main diagonal that is the reason why we apply diag(). We then substitute our result into equation (24)

Conclusion
We have properly formulated and implemented our spectral method based on the state-of-the-art rectangularization approach for the solution of oscillatory IVPs and nonlinear BVP in DEs. This approach requires the discretization of polynomials, forcing DEs into polynomial through the introduction of the re-sampling matrix Emphasis are now being placed on spectral methods when approximating for which we have introduced the concept of transforming the N into H using the relation (15) before discretizing the IVPs and BVP through rectangularization. In doing this, we observe spectral convergence unlike other conventional approaches. Our results of striking accuracy with efficient use of computer resources outperform the existing results (see Table 1, Table 2, Table 3). In particular, result as in Figure 3, Figure 4, Figure 5 shows the spectral accuracy and efficiency of the new rectangularization approach. It is our hope that the algorithm of emerging numerical codes for the discretization of DEs will be dominated by the used of the rectangularization approach just as it has already started in [7].