Abstract

The idea of rectangularization of a 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.

AMS Mathematics Subject Classification 2010: 65D05, 65D25, 65L05, 65L10.

1. 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 ¹ 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 ⁵ and finite Fourier transform in ². Compared with other formalizations of spectral methods such as the Galerkin approximations (see ⁶), they offer convenience in the presence of variable coefficients and non-linearities ⁷ 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 ⁸. 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 ⁹. 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 ⁹.

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 ⁷. Due to this snag, ⁷ introduced an efficient means of rectangularizing the 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 ⁷. However, the explicit form of the work in ⁷ can be found in ¹¹. Interestingly, emerging computational approach for reformulating all kinds of DEs to systems of first order equations as been established in ¹⁴.

2. 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

(1)

where

(2)

with property

(3)

As in ⁵, we define the barycentric weight as

(4)

and the barycentric formula of the first kind as

(5)

Equation (5), like the other authors is referred to as the improved Lagrange formula of (1). The formula (5) now demands floating point operations (flops) to evaluate the weights (4) and does not depend on data This feature allows for the evaluation of several functions, and it only requires operations each provided the (4) are known while Newton formula requires the computation of a new tableau for each new function. The formula (5) in addition, does not necessarily depend on the order in which the nodes are arranged, it can easily be updated recursively given new data set while the Newton formula depends on the order which increases the computational cost as increases, this many ordering leads to numerical instability ⁵. 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 ⁵). The barycentric formula of the second kind is given by

(6)

The barycentric Lagrange interpolation formula (6) is the one with special and beautiful symmetry ⁵. If weights (4) have common factor, they may be cancelled without any impact on ²⁰. Methods (6) and (5) both enjoy the advantage of updating weight (4) in operations to subsume a new set of The trigonometric form of formulas (4), (5), and (6) can be found in ⁵. 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 As in ³, the Chebyshev point of the first kind given by

(7)

and the Chebyshev point of the second kind is given by

(8)

The interpolation based on these points (7,8) converges as 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 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:

given two maps and Here, Otherwise one should use the barycentric formula of the first kind (5) instead. However, one attractive ingredient in the barycentric formula of the second kind (6) is that it generalizes to rational functions which has been a great tool for building Chebfun software package ²⁵.

3. Spectral Rectangularization through Collocation

Spectral rectangularization approach as demonstrated recently in the work in ⁷ is proved to be a convenient means of solving linear and non-linear DEs with general boundary conditions. This new approach in ⁷ 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 and the barycentric weights formula as

(9)

We consider a given data with a unique interpolating polynomial of degree express in the sense of formula (6) as

(10)

Since the interpolant depends linearly on the values , we can express the evaluation of the unique polynomial at a set of points vector in a matrix-vector multiplication form as

(11)

where the re-sampling matrix is given by

(12)

The vector now defines a unique polynomial interpolant of the data if and in equations (9) and (12) are replaced by and respectively ⁷. When , the re-sampling outcomes in an unchanged loss of information as a new vector of function values is shorter than the original. This process ⁷ 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 th-order differential operator is naturally discretized by an matrix, allowing the boundary constraints to be appended to form an invertible system ⁷. Naturally, the barycentric weights for these points are known as

(13)

Consider a differentiation matrix which the generalized form is given as

Then the differentiation matrix is pre-multiply by the re-sampling matrix (12). In particular, the case of a homogeneous differential equation of the form

(14)

The simple DE (14) can easily be discretized by setting it as a system

The second order differential operator (14) is rectangularized by a 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.

4. Numerical Experiment

In this section, we employ the relation where H, N are step-size and number of nodes respectively. The relation can be rearranged as

(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 N 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.

Example 4.1 Oscillatory DE, Source: ²⁷

(16)

where and . The exact solution is

Example 4.2 Complex DE, Source: ³⁰

(17)

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

(18)

The theoretical solution is given by We analyze problem (17) using our rectangular discretization approach in the interval [0,40] which corresponds to 20 orbits of the point at .

5. 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

(19)

(20)

where is the non-linear differential operator and is the non-linear differential functional. To achieve this task, we employ the Newton-iterative scheme for polynomial approximation as illustrated in ³¹. At in (19), we have the vector form as

(21)

On perturbing the solution u, say and expanding in Taylor's series about we have

(22)

(23)

(24)

We multiply the re-sampling matrix by (24).

(25)

where is a square matrix. Now, for the case of the system

(26)

We can express and and further express

(27)

where is the Hadamard point-wise or element-wise product. We then differentiate

(28)

(29)

(30)

(31)

(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) to have

(33)

where the Jacobian matrix and . More clearly,

(34)

Hence, we can solve the system (34).

Example 5.1 Consider a non-linear problem. Source: ³²

(35)

6. 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 non-linear BVP in DEs. This approach requires the discretization of polynomials, forcing DEs into polynomial through the introduction of the re-sampling matrix which re-samples the differentiated polynomials into a lower-degree subspace. The solutions of the DEs are then represented as finite vectors of their values at certain points. Each such vector implicitly defines a global interpolant, and operations on the DEs are replaced by the same operations applied to the corresponding interpolants.

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 ⁷.

Acknowledgements

The authors are thankful to Professor Neil Turok for his Next Einstein Initiative in Africa which has led to the establishment of the African Institute for Mathematical Sciences, Kigali City, Rwanda where the research work was carried out. Special thanks to instructor Mr. Edward Korveh for his guidance in the course of this research work.

References