New Spectral Second Kind Chebyshev Wavelets Algorithm for Solving Linear and Nonlinear Second-Order Differential Equations Involving Singular and Bratu Type Equations

Abstract

A new spectral algorithm based on shifted second kind Chebyshev wavelets operational matrices of derivatives is introduced and used for solving linear and nonlinear second-order two-point boundary value problems. The main idea for obtaining spectral numerical solutions for these equations is essentially developed by reducing the linear or nonlinear equations with their initial and/or boundary conditions to a system of linear or nonlinear algebraic equations in the unknown expansion coefficients. Convergence analysis and some efficient specific illustrative examples including singular and Bratu type equations are considered to demonstrate the validity and the applicability of the method. Numerical results obtained are compared favorably with the analytical known solutions.

1. Introduction

Spectral methods are one of the principal methods of discretization for the numerical solution of differential equations. The main advantage of these methods lies in their accuracy for a given number of unknowns (see, e.g., [1–3]). For smooth problems in simple geometries, they offer exponential rates of convergence/spectral accuracy. In contrast, finite difference and finite-element methods yield only algebraic convergence rates. The three most widely used spectral versions are the Galerkin, collocation, and tau methods. Collocation methods have become increasingly popular for solving differential equations, also they are very useful in providing highly accurate solutions to nonlinear differential equations (see, e.g., [4, 5]).

The subject of wavelets has recently drawn a great deal of attention from mathematical scientists in various disciplines. It is creating a common link between mathematicians, physicists, and electrical engineers. Wavelets theory is a relatively new and an emerging area in mathematical research. It has been applied to a wide range of engineering disciplines; particularly, wavelets are very successfully used in signal analysis for wave form representation and segmentations, time frequency analysis, and fast algorithms for easy implementation. Wavelets permit the accurate representation of a variety of functions and operators. Moreover, wavelets establish a connection with fast numerical algorithms (see [6, 7]).

The subject of nonlinear differential equations is a well-established part of mathematics, and its systematic development goes back to the early days of the development of calculus. Many recent advances in mathematics, paralleled by a renewed and flourishing interaction between mathematics, sciences, and engineering, have again shown that many phenomena in applied sciences, modelled by differential equations, will yield some mathematical explanation of these phenomena [8–10].

Even-order differential equations have been extensively discussed by a large number of authors due to their great importance in various applications in many fields. For example, in the sequence of papers [1, 2, 11], the Galerkin method is applied for handling these problems. They constructed suitable basis functions which satisfy the boundary conditions of the given differential equation. For this purpose, they used compact combinations of various orthogonal polynomials. The suggested algorithms in these papers are suitable for handling one- and two-dimensional linear high even-order boundary value problems.

In this paper, we aim to give some algorithms for handling both linear and nonlinear second-order boundary value problems based on introducing a new matrix of derivatives, then applying Petrov-Galerkin method on linear equations and collocation method on nonlinear equations. High accurate spectral wavelets solutions are achieved with a small number of retained modes compared with the usual spectral methods, also we can handle singular differential equations with discontinuous variable coefficients.

Among the second-order boundary value problems is the one-dimensional Bratu problem which has a long history. Bratu’s own paper appeared in 1914 [12]; generalizations are sometimes called the Liouville-Gelfand or Liouville-Gelfand-Bratu problem in honor of Gel’fand [13] and the nineteenth century work of the great French mathematician Liouville. In recent years, it has been a popular testbed for numerical and perturbation methods [14–16].

Simplification of the solid fuel ignition model in thermal combustion theory yields an elliptic nonlinear partial differential equation, namely, the Bratu problem. Also due to its use in a large variety of applications, many authors have contributed to the study of such problem. Some applications of Bratu problem are the model of thermal reaction process, the Chandrasekhar model of the expansion of the Universe, chemical reaction theory, nanotechnology, and radiative heat transfer (see [17–21]).

The application of Legendre wavelets for solving differential and integral equations is thoroughly considered by many authors (see, [22–26]). Also, Chebyshev wavelets are used for solving some differential, fractional and integral equations (see [5, 27–31]).

In [32, 33], a cosine and sine (CAS) wavelets operational matrix of fractional order integration has been derived and used to solve integrodifferential equations of fractional order.

One approach for solving differential equations is based on converting the differential equations into integral equations through integration, approximating various signals involved in the equation by truncated orthogonal series, and using the operational matrix of integration, to eliminate the integral operations.

Special attentions have been given to applications of block pulse functions [34], Legendre polynomials [35], Chebyshev polynomials [36], Haar wavelets [37], Legendre wavelets [24, 25, 38], and Chebyshev wavelets [27]. Another approach is based on using operational matrix of derivatives in order to reduce the problem into solving a system of algebraic equations (see [23]).

The main aim of this paper is to develop a new spectral algorithm for solving second-order two-point boundary value problems based on shifted second kind Chebyshev wavelets operational matrix of derivatives. The method reduces the differential equation with its initial and/or boundary conditions to a system of algebraic equations in the unknown expansion coefficients. Large systems of algebraic equations may lead to greater computational complexity and large storage requirements. However, the second kind Chebyshev wavelets are structurally sparse; this reduces drastically the computational complexity of solving the resulting algebraic system.

The structure of the paper is as follows. In Section 2, we give some relevant properties of second kind Chebyshev polynomials and their shifted forms. In Section 3, we develop a new shifted second kind Chebyshev wavelets operational matrices of derivatives (SCWOMD); also we ascertain the convergence analysis of the proposed algorithm. As an application of SCWOMD, numerical solutions of second-order linear and nonlinear initial value problems or two-point boundary value problems are implemented and presented in Section 4. In Section 5, some numerical examples are presented to show the efficiency and the applicability of the presented algorithm. Some concluding remarks are given in Section 6.

2. Some Properties of Second Kind Chebyshev Polynomials and Their Shifted Forms

In the present section, we discuss some relevant properties of the second kind Chebyshev polynomials and their shifted forms.

2.1. Second Kind Chebyshev Polynomials

It is well known that the second kind Chebyshev polynomials are defined on [−1,1] by 𝑈_𝑛(𝑥)=sin(𝑛+1)𝜃sin𝜃,𝑥=cos𝜃. (1) These polynomials are orthogonal on [−1,1]; that is, ¹₋₁√1−𝑥²𝑈_𝑚(𝑥)𝑈_𝑛(𝑥)𝑑𝑥=𝜋0,𝑚≠𝑛,2,𝑚=𝑛. (2) The following properties of second kind Chebyshev polynomials (see, e.g., [39]) are of fundamental importance in the sequel. They are eigenfunctions of the following singular Sturm-Liouville equation: 1−𝑥²𝐷²𝜙_𝑘(𝑥)−3𝑥𝐷𝜙_𝑘(𝑥)+𝑘(𝑘+2)𝜙_𝑘(𝑥)=0, (3) where 𝐷≡𝑑/𝑑𝑥 and may be generated by using the recurrence relation 𝑈_𝑘+1(𝑥)=2𝑥𝑈_𝑘(𝑥)−𝑈_𝑘−1(𝑥),𝑘=1,2,3,…, (4) starting from 𝑈₀(𝑥)=1and 𝑈₁(𝑥)=2𝑥, or from Rodrigues formula 𝑈_𝑛(𝑥)=(−2)^𝑛(𝑛+1)!(2𝑛+1)!√1−𝑥²𝐷^𝑛1−𝑥^2𝑛+(1/2). (5) The following theorem is needed hereafter.

Theorem 1. The first derivative of second kind Chebyshev polynomials is given by 𝐷𝑈_𝑛(𝑥)=2^𝑛−1_{()𝑘=0𝑘+𝑛𝑜𝑑𝑑}(𝑘+1)𝑈_𝑘(𝑥). (6)

For a proof of Theorem 1, see [39].

2.2. Shifted Second Kind Chebyshev Polynomials

The shifted second kind Chebyshev polynomials are defined on [0,1] by 𝑈^∗_𝑛(𝑥)=𝑈_𝑛(2𝑥−1). All results of second kind Chebyshev polynomials can be easily transformed to give the corresponding results for their shifted forms. The orthogonality relation with respect to the weight function √𝑥−𝑥² is given by ¹₀√𝑥−𝑥²𝑈^∗_𝑛(𝑥)𝑈^∗_𝑚(𝑥)𝑑𝑥=𝜋0,𝑚≠𝑛,8,𝑚=𝑛. (7) The first derivative of 𝑈^∗_𝑛(𝑥) is given in the following corollary.

Corollary 2. The first derivative of the shifted second kind Chebyshev polynomial is given by 𝐷𝑈^∗_𝑛(𝑥)=4^𝑛−1_{()𝑘=0𝑘+𝑛𝑜𝑑𝑑}(𝑘+1)𝑈^∗_𝑘(𝑥). (8)

3. Shifted Second Kind Chebyshev Operational Matrix of Derivatives

Wavelets constitute of a family of functions constructed from dilation and translation of single function called the mother wavelet. When the dilation parameter 𝑎 and the translation parameter 𝑏 vary continuously, then we have the following family of continuous wavelets: 𝜓_𝑎,𝑏(𝑡)=|𝑎|^−1/2𝜓𝑡−𝑏𝑎,𝑎,𝑏∈ℝ,𝑎≠0. (9) Second kind Chebyshev wavelets 𝜓_𝑛𝑚(𝑡)=𝜓(𝑘,𝑛,𝑚,𝑡) have four arguments: 𝑘, 𝑛 can assume any positive integer, 𝑚 is the order of second kind Chebyshev polynomials, and 𝑡 is the normalized time. They are defined on the interval [0,1] by 𝜓_𝑛𝑚(𝑡)=2^(𝑘+3)/2√𝜋𝑈^∗_𝑚2^𝑘𝑛𝑡−𝑛,𝑡∈2^𝑘,𝑛+12^𝑘,0,otherwise,𝑚=0,1,…,𝑀,𝑛=0,1,…,2^𝑘−1. (10)

3.1. Function Approximation

A function 𝑓(𝑡) defined over [0,1] may be expanded in terms of second kind Chebyshev wavelets as 𝑓(𝑡)=^∞∞_𝑛=0𝑚=0𝑐_𝑛𝑚𝜓_𝑛𝑚(𝑡), (11) where 𝑐_𝑛𝑚=𝑓(𝑡),𝜓_𝑛𝑚(𝑡)_𝜔=¹₀𝜔(𝑡)𝑓(𝑡)𝜓_𝑛𝑚(𝑡)𝑑𝑡, (12) and 𝑤(𝑡)=√𝑡−𝑡². If the infinite series is truncated, then 𝑓(𝑡) can be approximated as 𝑓(𝑡)≃^{2𝑘−1𝑀}_𝑛=0𝑚=0𝑐_𝑛𝑚𝜓_𝑛𝑚(𝑡)=𝐶^𝑇Ψ(𝑡), (13) where 𝐶 and Ψ(𝑡) are 2^𝑘(𝑀+1)×1 matrices defined by 𝐶=𝑐_0,0,𝑐_0,1,…,𝑐_0,𝑀,…,𝑐₂^𝑘_−1,𝑀,…,𝑐₂^𝑘_−1,1,…,𝑐₂^𝑘_−1,𝑀^𝑇,Ψ(𝑡)=𝜓_0,0,𝜓_0,1,…,𝜓_0,𝑀𝜓,…,_{2𝑘−1,𝑀},…,𝜓_2𝑘−1,1,…,𝜓_{2𝑘−1,𝑀}^𝑇. (14)(15)

3.2. Convergence Analysis

We state and prove a theorem ascertaining that the second kind Chebyshev wavelet expansion of a function 𝑓(𝑥), with bounded second derivative, converges uniformly to 𝑓(𝑥).

Theorem 3. A function 𝑓(𝑥)∈𝐿²_𝜔[0,1], with |𝑓^′′(𝑥)|⩽𝐿, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to 𝑓(𝑥). Explicitly, the expansion coefficients in (12) satisfy the following inequality: 𝑐_𝑛𝑚<8√2𝜋𝐿(𝑛+1)^5/2(𝑚+1)²,∀𝑚>1,𝑛⩾0. (16)

Proof. From (12), it follows that 𝑐_𝑛𝑚=2^(𝑘+3)/2√𝜋^{(𝑛+1)/2𝑘}_𝑛/2^𝑘𝑓(𝑥)𝑈^∗_𝑚2^𝑘𝑥−𝑛𝜔2^𝑘𝑥−𝑛𝑑𝑥. (17) If we set 2^𝑘𝑥−𝑛=cos𝜃 in (17), then we get 𝑐_𝑛𝑚=2^{(−𝑘+3)/2}√𝜋×^𝜋₀𝑓cos𝜃+𝑛2^𝑘()=2sin𝑚+1𝜃sin𝜃𝑑𝜃_{(−𝑘+1)/2}√𝜋×^𝜋₀𝑓cos𝜃+𝑛2^𝑘[cos𝑚𝜃−cos(𝑚+2)𝜃]𝑑𝜃, (18) which gives after integration by parts two times 𝑐_𝑛𝑚=22^5𝑘/2√2𝜋^𝜋0𝑓^′′cos𝜃+𝑛2^𝑘𝛾_𝑚(𝜃)𝑑𝜃, (19) where 𝛾_𝑚(𝜃)=sin𝜃𝑚sin(𝑚−1)𝜃−()𝜃𝑚−1sin𝑚+1−𝑚+1sin𝜃()𝜃𝑚+2sin𝑚+1−()𝜃𝑚+1sin𝑚+3.𝑚+3 (20) Therefore, we have 𝑐_𝑛𝑚=12^{(5𝑘−1)/2}√𝜋^𝜋₀𝑓^′′cos𝜃+𝑛2^𝑘𝛾_𝑚(𝜃)=1𝑑𝜃2^{(5𝑘−1)/2}√𝜋^𝜋₀𝑓^′′cos𝜃+𝑛2^𝑘𝛾_𝑚(𝜃)⩽𝐿𝑑𝜃2^{(5𝑘−1)/2}√𝜋^𝜋₀𝛾_𝑚(𝜃⩽𝐿√)𝑑𝜃𝜋2^{(5𝑘−1)/2}×1𝑚1+1𝑚−1+1𝑚+11𝑚+2+1𝑚+1=𝐿√𝑚+3𝜋2^{(5𝑘−5)/2}𝑚²<√+2𝑚−32𝐿𝜋2^{(5𝑘−5)/2}(𝑚+1)². (21) Since 𝑛⩽2^𝑘−1, we have 𝑐_𝑛𝑚<8√2𝜋𝐿(𝑛+1)^5/2(𝑚+1)², (22) and this completes the proof of Theorem 3.

A shifted second kind Chebyshev wavelets operational matrix of the first derivative is stated and proved in the following theorem.

Theorem 4. Let Ψ(𝑡) be the second kind Chebyshev wavelets vector defined in (15). Then, the first derivative of the vector Ψ(𝑡) can be expressed as 𝑑Ψ(𝑡)(𝑡),𝑑𝑡=𝐃Ψ (23) where 𝐃 is 2^𝑘(𝑀+1) square operational matrix of derivative and is defined by 𝐃=𝐹0...00𝐹...0............00...𝐹, (24) in which 𝐹 is an (𝑀+1) square matrix whose (𝑟,𝑠)th element is defined by 𝐹_𝑟,𝑠=2^𝑘+2𝑠,𝑟⩾2,𝑟>𝑠,(𝑟+𝑠)𝑜𝑑𝑑,0,𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒. (25)

Proof. If we make use of the shifted second kind Chebyshev polynomials, then the 𝑟th element of the vector Ψ(𝑡) defined in (15) can be written in the following form: Ψ_𝑟(𝑡)=𝜓_𝑛𝑚(𝑡)=2^(𝑘+3)/2√𝜋𝑈^∗_𝑚2^𝑘𝑡−𝑛𝜒_{[𝑛/2𝑘,(𝑛+1)/2𝑘]},𝑖=1,2,…,2^𝑘(𝑀+1), (26) where 𝑟=𝑛(𝑀+1)+(𝑚+1), 𝑚=0,1,…,𝑀, 𝑛=0,1,…,(2^𝑘−1), and 𝜒_{[𝑛/2𝑘,(𝑛+1)/2𝑘]}=𝑛1,𝑡∈2^𝑘,𝑛+12^𝑘,0,otherwise. (27) If we differentiate (26) with respect to 𝑡, then we get 𝑑Ψ_𝑟(𝑡)=2𝑑𝑡√𝜋^(2𝑘+3)/2𝑑𝑈𝑑𝑡^∗𝑚2^𝑘𝑡−𝑛𝜒_{[𝑛/2𝑘,(𝑛+1)/2𝑘]}. (28) It is to be noted here that the R.H.S. of (28) is zero outside the interval [𝑛/2^𝑘,(𝑛+1)/2^𝑘]; hence, its second kind Chebyshev wavelets expansion only has those elements in Ψ(𝑡) that are nonzero in the interval [𝑛/2^𝑘,(𝑛+1)/2^𝑘]; that is, Ψ_𝑖(𝑡), 𝑖=𝑛(𝑀+1)+1,𝑛(𝑀+1)+2,…,𝑛(𝑀+1)+𝑀+1, and as a consequence, this enables one to expand (𝑑Ψ_𝑟(𝑡)/𝑑𝑡) in terms of shifted second kind Chebyshev wavelets Ψ_𝑖(𝑡), 𝑖=𝑛(𝑀+1)+1,𝑛(𝑀+1)+2,…,𝑛(𝑀+1)+𝑀+1, in the form 𝑑Ψ_𝑟(𝑡)=𝑑𝑡^{(𝑛+1)(𝑀+1)}_{()𝑖=𝑛𝑀+1+1}𝑎_𝑖Ψ_𝑖(𝑡), (29) and accordingly (28) implies that the operational matrix 𝐃 is a block matrix defined as in (24). Moreover, we have 𝑑𝑈^∗₀(𝑡)𝑑𝑡=0, (30) and this implies that 𝑑Ψ𝑑𝑡^𝑟(𝑡)()()2=0for𝑟=1,𝑀+1+1,2𝑀+1+1,…,_𝑘−1(𝑀+1)+1. (31) Consequently, the elements of the first row of matrix 𝐹 defined in (24) are zeros. Now and with the aid of Corollary 2, the first derivative of Ψ_𝑟(𝑡) may be expressed in the form 𝑑Ψ_𝑟(𝑡)=2𝑑𝑡_(2𝑘+7)/2√𝜋^𝑚−1_{𝑗=0(𝑗+𝑚)odd}𝑗+1𝑈^∗_𝑗2^𝑘𝑡−𝑛𝜒_{[𝑛/2𝑘,(𝑛+1)/2𝑘]}. (32) Expanding this equation in terms of second kind Chebyshev wavelets basis, we get 𝑑Ψ_𝑟(𝑡)𝑑𝑡=2^{𝑘+2𝑚−1}_{()𝑠=0𝑠+𝑚odd}𝑠Ψ_{𝑛(𝑀+1)+𝑠}(𝑡), (33) so if 𝐹_𝑟,𝑠 are chosen as 𝐹_𝑟,𝑠=2^𝑘+2()𝑠,𝑟=2,…,𝑀+1,𝑠=1,…,𝑟−1,𝑟+𝑠odd,0,otherwise, (34) then (23) is obtained, and the proof of the theorem is completed.

Corollary 5. The operational matrix for the 𝑛th derivative can be obtained from 𝑑^𝑛Ψ(𝑡)𝑑𝑡^𝑛=𝐃^𝑛Ψ(𝑡),𝑛=1,2,…, (35) where 𝐃^𝑛 is the 𝑛th power of matrix 𝐃.

4. Second-Order Two-Point Boundary Value Problems

In this section, we are interested in solving linear and nonlinear two-point boundary value problems subject to homogenous or nonhomogenous initial or boundary conditions based on the wavelets operational matrices constructed in Section 3.

4.1. Linear Second-Order Two-Point Boundary Value Problems

Consider the linear second-order differential equation 𝑦^′′(𝑥)+𝑓₁(𝑥)𝑦^′(𝑥)+𝑓₂(𝑥)𝑦(𝑥)=𝑔(𝑥),𝑥∈[0,1], (36) subject to the initial conditions 𝑦(0)=𝛼,𝑦^′(0)=𝛽, (37) or the boundary conditions 𝑦(0)=𝛼,𝑦(1)=𝛽, (38) or the most general mixed boundary conditions 𝑎₁𝑦(0)+𝑎₂𝑦^′(0)=𝛼,𝑏₁𝑦(1)+𝑏₂𝑦^′(1)=𝛽. (39) If we approximate 𝑦(𝑥), 𝑓₁(𝑥), 𝑓₂(𝑥), and 𝑔(𝑥) in terms of the second kind Chebyshev wavelets basis, then one can write 𝑦(𝑥)≃^{2𝑘−1𝑀}_𝑛=0𝑚=0𝑐_𝑛𝑚𝜓_𝑛𝑚(𝑥)=𝐶^𝑇Ψ(𝑥),𝑓₁(𝑥)≃^{2𝑘−1𝑀}_𝑛=0𝑚=0𝑓_𝑛𝑚𝜓_𝑛𝑚(𝑥)=𝐹^𝑇₁Ψ(𝑥),𝑓₂(𝑥)≃^{2𝑘−1𝑀}_𝑛=0𝑚=0𝑓_𝑛𝑚𝜓_𝑛𝑚(𝑥)=𝐹^𝑇₂Ψ(𝑡),𝑔(𝑥)≃^{2𝑘−1𝑀}_𝑛=0𝑚=0𝑔_𝑛𝑚𝜓_𝑛𝑚(𝑥)=𝐺^𝑇Ψ(𝑥), (40) where 𝐶^𝑇, 𝐹^𝑇₁, 𝐹^𝑇₂, and 𝐺^𝑇 are defined similarly as in (14). Relations (23) and (35) enable one to approximate 𝑦^′(𝑥) and 𝑦^′′(𝑥) as 𝑦^′(𝑥)≃𝐶^𝑇𝐃Ψ(𝑥),𝑦^′′(𝑥)≃𝐶^𝑇𝐃²Ψ(𝑥). (41) Now, substitution of relations (40) and (41) into (36) enables us to define the residual, 𝑅(𝑥), of this equation as 𝑅(𝑥)=𝐶^𝑇𝐃²Ψ(𝑥)+𝐹^𝑇₁Ψ(𝑥)(Ψ(𝑥))^𝑇𝐃^𝑇𝐶+𝐹^𝑇₂Ψ(𝑥)(Ψ(𝑥))^𝑇𝐶−𝐺^𝑇Ψ(𝑥), (42) and application of the Petrov-Galerkin method (see [40]) yields the following (2^𝑘(𝑀+1)−2) linear equations in the unknown expansion coefficients, 𝑐_𝑛𝑚; namely, ¹₀√𝑥−𝑥²Ψ_𝑗(𝑥)𝑅(𝑥)𝑑𝑥=0,𝑗=1,2,…,2^𝑘(𝑀+1)−2. (43) Moreover, the initial conditions (37), the boundary conditions (38), and the mixed boundary conditions (39) lead, respectively, to the following equations: 𝐶^𝑇Ψ(0)=𝛼,𝐶^𝑇(0)𝐶𝐃Ψ=𝛽,_𝑇Ψ(0)=𝛼,𝐶^𝑇Ψ(1)𝑎=𝛽,₁𝐶^𝑇Ψ(0)+𝑎₂𝐶^𝑇(0)𝑏𝐃Ψ=𝛼,₁𝐶^𝑇Ψ(1)+𝑏₂𝐶^𝑇𝐃Ψ(1)=𝛽. (44)(45)(46) Thus, (43) with the two equations of (44) or (45) or (46) generates 2^𝑘(𝑀+1) a set of linear equations which can be solved for the unknown components of the vector 𝐶, and hence an approximate spectral wavelets solution to 𝑦(𝑥) given in (40) can be obtained.

4.2. Nonlinear Second-Order Two-Point Boundary Value Problems

Consider the nonlinear differential equation 𝑦^′′(𝑥)=𝐹𝑥,𝑦(𝑥),𝑦^′(𝑥), (47) subject to the initial conditions 𝑦(0)=𝛼,𝑦^′(0)=𝛽, (48) or the boundary conditions 𝑦(0)=𝛼,𝑦(1)=𝛽, (49) or the most general mixed boundary conditions 𝑎₁𝑦(0)+𝑎₂𝑦^′(0)=𝛼,𝑏₁𝑦(1)+𝑏₂𝑦^′(1)=𝛽. (50) If we follow the same procedure of Section 4.1, approximate 𝑦(𝑥) as in (40) and make use of (23) and (35), then we obtain 𝐶^𝑇𝐃²Ψ(𝑥)=𝐹𝑥,𝑦(𝑥),𝐶^𝑇𝐃Ψ(𝑥). (51)

To find an approximate solution to 𝑦(𝑥), we compute (51) at the first 2^𝑘(𝑀+1)−2 roots of 𝑈^∗_{2𝑘(𝑀+1)}(𝑥). These equations with (44) or (45) or (46) generate 2^𝑘(𝑀+1) nonlinear equations in the expansion coefficients, 𝑐_𝑛𝑚, which can be solved with the aid of Newton’s iterative method.

5. Numerical Results and Discussions

In this section, the presented algorithms in Section 4 are applied to solve linear and nonlinear second-order initial and boundary value problems as well as Bratu’s equations. Some examples subject to different initial, boundary, and mixed boundary conditions are considered to illustrate the efficiency and the applicability of the proposed algorithms. All computations are performed by using Mathematica 8.

Example 1. Consider the following linear second-order boundary value problem: 𝑦^′′(𝑥)+(1−𝑥)𝑦^′(𝑥)+2𝑦(𝑥)=1+2𝑥−𝑥²(0)(1)sin𝑥,0<𝑥<1,𝑦=1,𝑦=0, (52) with the exact solution 𝑦(𝑥)=(1−𝑥)cos𝑥.
In Table 1, the maximum absolute error 𝐸 is listed for various values of 𝑘 and 𝑀.

Table 1 

Maximum absolute error 𝐸 for Example 1.

Example 2. Consider the following linear second-order boundary value problem: 𝑦^′′(𝑥)−sin𝑥𝑦^′(𝑥)+𝑦(𝑥)=sin²𝑥;0<𝑥<1, (53) subject to the mixed boundary conditions 𝑦(0)+𝑦^′(0)=1,sin(1)𝑦(1)+cos(1)𝑦^′(1)=0, (54) with the exact solution 𝑦(𝑥)=cos𝑥.
In Table 2, the maximum absolute error 𝐸 is listed for various values of 𝑘 and 𝑀.

Table 2 

Maximum absolute error 𝐸 for Example 2.

Example 3. Consider the following nonlinear initial value problem (see [41]): 4𝑦^′′(𝑥)−2𝑦^′(𝑥)²(𝑥)(0)+𝑦=0,0<𝑥<1,𝑦=−1,𝑦^′(0)=0, (55) with the exact solution 𝑦(𝑥)=(𝑥²/8)−1.
We solve (55) using the algorithm described in Section 4.2 for the case corresponds to 𝑀=2, 𝑘=0 to obtain an approximate solution of 𝑦(𝑥). First, if we make use of (24) and (35), then the two operational matrices 𝐃 and 𝐃² are given, respectively, by 𝐃=000400080,𝐃²=0000003200. (56) Moreover, Ψ(𝑥) can be evaluated to give Ψ(𝑥)=2𝜋28𝑥−432𝑥²−32𝑥+6. (57) If we set 𝐶=𝑐_0,0,𝑐_0,1,𝑐_0,2^𝑇=𝜋2𝑐₀,𝑐₁,𝑐₂^𝑇, (58) then (55) takes the form 4𝐶^𝑇𝐃²Ψ(𝑥)−2𝐶^𝑇𝐃Ψ(𝑥)²+𝐶^𝑇Ψ(𝑥)=0, (59) which is equivalent to 𝑐₀−2𝑐₁+131𝑐₂−64𝑐²₁+512𝑐₁𝑐₂−1024𝑐²₂+4𝑐₁−4𝑐₂−256𝑐₁𝑐₂+1024𝑐²₂𝑥+16𝑐₂−256𝑐²₂𝑥²=0. (60) We only need to satisfy this equation at the first root of 𝑈^∗₃(𝑥), that is, at 𝑥₁=(2−√2)/4, to get 𝑐₀−64𝑐²₁+√2𝑐₁256𝑐₂−1+129−512𝑐₂𝑐₂=0. (61) Furthermore, the use of initial conditions in (55) leads to the two equations 2𝑐₀−4𝑐₁+6𝑐₂𝑐=−1,₁−4𝑐₂=0. (62) The solution of the nonlinear system of (61)-(62) gives 𝑐₀=−123256,𝑐₁=164,𝑐₂=1,256 (63) and consequently 𝑦(𝑥)=−123,1256,16422568𝑥−432𝑥²=𝑥−32𝑥+6₂8−1, (64) which is the exact solution.

Example 4. Consider the following nonlinear second-order boundary value problem: 𝑦^′′(𝑥)+𝑦^′(𝑥)²−𝑦(𝑥)=sinh²𝑦(0)(1)𝑥,=1,𝑦=cosh1, (65) with the exact solution 𝑦(𝑥)=cosh𝑥.
In Table 3, we list the maximum absolute error 𝐸, for various values of 𝑘 and 𝑀.

Table 3 

Maximum absolute error 𝐸 for Example 4.

Example 5. Consider the following linear singular initial value problem (see [23]): 𝑦^′′(𝑥)+𝑓(𝑥)𝑦^′(𝑥)(𝑥)(𝑥),𝑦(0)+𝑦=𝑔=0,𝑦^′(0)=1, (66) where 𝑓(𝑥)=4−4𝑥+31,0⩽𝑥<4,1𝑥−3,14⩽𝑥⩽1, (67) and 𝑔(𝑥) is chosen such that the exact solution of (66) is 𝑦(𝑥)=𝑒^𝑥sin𝑥/(1+𝑥²). We apply the algorithm described in Section 4.1 with 𝑘=2. In Table 4, we list the maximum absolute error for various values of 𝑀, while in Figure 1, the exact solution is compared with the three numerical solutions corresponding to 𝑘=2, 𝑀=2, 𝑘=2, 𝑀=4, 𝑘=2,  and 𝑀=6 and also in Table 5 we give a comparison between the error obtained in [23] and the present method, 𝑘=2, 𝑀=8.

Table 4 

Maximum absolute error 𝐸 for Example 5.

Table 5 

Comparison between different errors for Example 5.

Figure 1 

Numerical and exact solutions of Example 5.

Example 6. Consider the following Bratu equation (see [17–20]): 𝑦^′′(𝑥)+𝜆𝑒^𝑦(𝑥)=0,𝑦(0)=𝑦(1)=0,0⩽𝑥⩽1, (68) with the analytical solution 𝑦(𝑥)=−2lncosh((𝜃/4)(2𝑥−1))(),cosh𝜃/4 (69) where 𝜃 is the solution of the nonlinear equation 𝜃=√2𝜆cosh𝜃. The presented algorithm in Section 4.2 is applied to numerically solve (68), for the three cases corresponding to 𝜆=1, 2 and 3.51 which yield 𝜃=1.51716, 2.35755, and 4.66781, respectively. In Table 6, the maximum absolute error 𝐸 is listed for 𝑘=0 and various values of 𝑀, while in Table 7, we give a comparison between the best errors obtained by our method and some other numerical methods used to solve Example 6. This table shows that our method is more accurate if compared with the methods developed in [17–20]. In addition, Figure 2 illustrates a comparison between different solutions obtained by our algorithm in case of 𝜆=1, 𝑘=0, and 𝑀=1,2,3.

Table 6 

Maximum absolute error 𝐸 for Example 6.

Table 7 

Comparison between the best errors for Example 6 for 𝑘=0.

Figure 2 

Different solutions of Example 6.

Remark 6. It is worth noting here that the obtained numerical results in the previous six solved examples are very accurate, although the number of retained modes in the spectral expansion is very few, and again the numerical results are compared favorably with the known analytical solutions.

6. Concluding Remarks

In this paper, an algorithm for obtaining a numerical spectral solution for second-order linear and nonlinear boundary value problems is discussed. The derivation of this algorithm is essentially based on constructing the shifted second kind Chebyshev wavelets operational matrix of differentiation. One of the main advantages of the presented algorithm is its availability for application on both linear and nonlinear second-order boundary value problems including some singular equations and also Bratu type equations. Another advantage of the developed algorithm is that high accurate approximate solutions are achieved using a small number of the second kind Chebyshev wavelets. The obtained numerical results are compared favorably with the analytical ones.

Acknowledgments

This paper was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah. The authors, therefore, acknowledge with thanks DSR technical and financial support. Also, they would like to thank the referee for his valuable comments and suggestions which improved the paper into its present form.