A matrix formulation of the Tau method for the numerical solution of non-linear problems
Kourosh Parand, Amin Ghaderi, Mehdi Delkhosh, Reza Pourgholi
aa r X i v : . [ m a t h . NA ] A ug Noname manuscript No. (will be inserted by the editor)
A matrix formulation of the Tau method for thenumerical solution of non-linear problems
Kourosh Parand · Amin Ghaderi · Mehdi Delkhosh · Reza Pourgholi
Received: date / Accepted: date
Abstract
The purpose of this research is to propose a new approach namedthe shifted Bessel Tau (SBT) method for solving higher-order ordinary dif-ferential equations (ODE). The operational matrices of derivative, integraland product of shifted Bessel polynomials on the interval [ a, b ] are calculated.These matrices together with the Tau method are utilized to reduce the solu-tion of the higher-order ODE to the solution of a system of algebraic equationswith unknown Bessel coefficients. The comparisons between the results of thepresent work and other the numerical method are shown that the present workis computationally simple and highly accurate.
Keywords
Bessel polynomials · Nonlinear ODE · Tau method · Operationalmatrices · Boundary value problems
Mathematics Subject Classification (2000) · · · K. ParandDepartment of Computer Sciences, Shahid Beheshti University, G.C., Tehran, Iran.Department of Cognitive Modelling, Institute for Cognitive and Brain Sciences, ShahidBeheshti University, G.C, Tehran, Iran.E-mail: k [email protected]. GhaderiDepartment of Cognitive Modelling, Institute for Cognitive and Brain Sciences, ShahidBeheshti University, G.C, Tehran, Iran.E-mail: [email protected]. DelkhoshDepartment of Mathematics and Computer Science, Islamic Azad University, BardaskanBranch, Bardaskan, Iran.E-mail: [email protected]. PourgholiSchool of Mathematics and Computer Sciences, Damghan University, Damghan, Iran.E-mail: [email protected] Kourosh Parand et al.
Differential equations and their solutions play a major role in science andengineering. A physical event can be modeled by the differential equation,an integral equation or an integro-differential equation or a system of theseequations. Since few of these equations cannot be solved explicitly, it is of-ten necessary to resort to the numerical techniques which are appropriatecombinations of numerical integration and interpolation. The solution of thisequation occurring in physics, biology and engineering are based on numeri-cal methods. One of the powerful method for solving the ordinary differentialequations, integral and integro-differential equations are spectral methods inboth terms of accuracy and simplicity. Spectral methods, in the context of nu-merical schemes for differential equations, generically belong to the family ofweighted residual methods (WRMs) [1]. WRMs represent a particular groupof approximation techniques, in which the residuals (or errors) are minimizedin a certain way and thereby leading to specific methods, including Galerkin,Petrov-Galerkin, collocation and Tau formulations [2,3,4,5]. The basis of spec-tral method is always polynomials (orthogonal polynomials). These types ofbasis for approximation of enough smooth functions in the finite domain leadsto the exponential convergence.In Tau method, the operational matrices of derivative, integral and productof basis functions are obtained and then these matrices together to reduce thesolution of these problems to the solution of a system of algebraic equations.Tau method is based on expanding the required approximate solution as theelements of a complete set of basis functions. Polynomial series and orthog-onal functions have received considerable attention in dealing with variousproblems of differential, integral and integro-differential equations. Polynomi-als are incredibly useful mathematical tools as they are simply defined andcan be calculated quickly on computer systems and represent a tremendousvariety of functions. They can be differentiated and integrated easily, and canbe pieced together to form spline curves that can approximate any functionto any accuracy desired. Also, the main characteristic of this technique is thatit reduces these problems to those of solving a system of algebraic equations,thus greatly simplifies the problems.In recent years, the differential, integral and integro-differential equations havebeen solved using the homotopy perturbation method [8,9], the fractional or-der of rational Bessel functions collocation (FRBC) method [54], the radialbasis functions method [10], the collocation method [11], the Homotopy anal-ysis method [12,13], the Tau method [14,15], the Variational iteration method[16,17], the Legendre matrix method [7,18], the Haar Wavelets operationalmatrix [19], the Shifted Chebyshev direct method [20], the generalized frac-tional order of the Chebyshev collocation approach [42], the Legendre Waveletsoperational matrix [21], Sine-Cosine Wavelets operational matrix [22], the op-erational matrices of Bernstein polynomials [23]. Razzaghi et al. in [24,25,26]presented the integral and product operational matrices based on the Fourier,the Taylor and the shifted-Jacobi series. Parand et al. in [27,28,29,30] ap- au method for the numerical solution of non-linear problems 3 plied the spectral method to solve nonlinear ordinary differential equations onsemi-infinite intervals. Their approach was based on a rational Tau method.They obtained the operational matrices of derivative and product of rationalChebyshev and Legendre functions and then applied these matrices togetherwith the Tau method to reduce the solution of these problems to the solutionof a system of algebraic equations. Doha [31] has derived the shifted Jacobioperational matrix of fractional derivatives which is applied together with thespectral Tau method for the numerical solution of dynamical systems.Yousefiet al. in [32,33,34,35] have presented Legendre wavelets and Bernstein opera-tional matrices and used them to solve miscellaneous systems. Also, recently,in [36] a new shifted-Chebyshev operational matrix of fractional integrationof arbitrary order is introduced and applied together with the spectral Taumethod for solving linear fractional differential equations (FDEs). In [37,38]the operational matrix form in a Hybrid of block-pulse functions and anotherset of functions like Taylor series, Legendre and Chebyshev has been usedfor finding the solution of the various classes of dynamical systems. Lakestani[39] constructed the operational matrix of fractional derivative of order α inthe Caputo sense using the linear B-spline functions. In [40], a general for-mulation for the d dimensional orthogonal functions and their derivative andproduct matrices has been presented. Ordinary operational matrices have beenutilized together with the Tau method to reduce the solution of partial dif-ferential equations (PDEs) to the solution of a system of algebraic equations.Recently in [41], a class of two dimensional nonlinear Volterra’s integral equa-tions has been solved using operational matrices of Legendre polynomials. Theoperational matrices of integration and product together with the collocationpoints have been utilized to reduce the solution of the integral equation to thesolution of a nonlinear algebraic equation system.The remainder of the paper is organized as follows: In section 2, the func-tion approximation and convergence of the method by using Bessel polyno-mials are illustrated . In sections 3, the general procedure for formulation ofoperational matrices of integration, differentiation, and product are explainedand are proved. In section 4, our numerical findings and demonstrate the va-lidity, accuracy and applicability of the operational matrices by consideringnumerical examples are reported. Also a conclusion is given in the last sec-tion. The Bessel functions arise in many problems in physics possessing cylindricalsymmetry, such as the vibrations of circular drumheads and the radial modesin optical fibers. Bessel functions are usually defined as a particular solutionof a linear differential equation of the second order which known as Bessel’sequation.
Kourosh Parand et al.
Bessel differential equation of order n ∈ R is: x d y ( x ) dx + x dy ( x ) dx + ( x − n ) y ( x ) = 0 , x ∈ ( −∞ , ∞ ) . (1)One of the solutions of Eq. (1) by applying the method of Frobenius is asfollows [43]: J n ( x ) = ∞ X r =0 ( − r r !( n + r )! ( x r + n , (2)where series (2) is convergent for all x ∈ ( −∞ , ∞ ).Bessel functions and polynomials are used to solve the more number ofproblems in physics, engineering, mathematics, and etc., such as Blasius equa-tion, Lane-Emden equations, integro-differential equations of the fractionalorder, unsteady gas equation, systems of linear Volterra integral equations,high-order linear complex differential equations in circular domains, systemsof high-order linear Fredholm integro-differential equations, nolinear Thomas-Fermi on semi-infinite domain, etc. [44,45,46,47,48,49,50,51,52,53,54].Bessel polynomials has been introduced as follows [49]: B n ( x ) = [ N − n ] X r =0 ( − r r !( n + r )! ( x r + n , x ∈ [0 , . where n ∈ N , and N is the number of basis of Bessel polynomials.Now, We introduce new shifted Bessel equation in the interval [ a, b ] by: Q n ( x ) = B n ( x − ab − a ) , x ∈ [ a, b ] , n = 0 , , , ..., N. where a, b ∈ R .2.1 Approximation of functionsLet us define Γ = { x | a ≤ x ≤ b } and L w ( Γ ) = { v : Γ → R | v is measurable and k v k w < ∞} , where k v k w = Z ba | v ( x ) | w ( x ) dx ! / , with w ( x ) = b − a , is the norm induced by inner product of the space L w ( Γ )as follows: h v, u i w = Z ba v ( x ) u ( x ) w ( x ) dx. Now, suppose that Q N = span { Q ( x ) , Q ( x ) , . . . , Q N ( x ) } , au method for the numerical solution of non-linear problems 5 where Q N is a finite-dimensional subspace of L w ( Γ ), dim ( Q N ) = N + 1, so Q N is a closed subspace of L w ( Γ ). Therefore, Q N is a complete subspace of L w ( Γ ). Assume that f ( x ) be an arbitrary element in L w ( Γ ). Thus f ( x ) has aunique best approximation in Q N subspace, say ˆ b ∈ Q N , that is ∃ ˆ b ( x ) ∈ Q N , ∀ b ( x ) ∈ Q N , k f ( x ) − ˆ b ( x ) k w ≤k f ( x ) − b ( x ) k w . Notice that we can write b ( x ) function as a linear combination of the basisvectors of Q N subspace.We know function of f ( x ) can be expanded by N + 1 terms of Besselpolynomials as: f ( x ) = f N ( x ) + R ( x ) , that is f N = N X n =0 a n Q n ( x ) = A T Q ( x ) , (3)where Q ( x ) = [ Q ( x ) , Q ( x ) , · · · , Q N ( x )] T and R ( x ) ∈ Q ⊥ N that Q ⊥ N is theorthogonal complement. So f ( x ) − f N ( x ) ∈ Q ⊥ N and b ( x ) ∈ Q N are orthogonalwhich we denote it by: (cid:0) f ( x ) − f N ( x ) (cid:1) ⊥ b ( x ) , thus f ( x ) − f N ( x ) vector is orthogonal over all of basis vectors of Q N subspaceas: h f ( x ) − f N ( x ) , Q i ( x ) i w = h f ( x ) − A T Q ( x ) , Q i ( x ) i w = 0 , i = 0 , , · · · , N, hence h f ( x ) − A T Q ( x ) , Q T ( x ) i w = 0 , therefore A can be obtained by h f ( x ) , Q T ( x ) i w = h A T Q ( x ) , Q T ( x ) i w ,A T = h f ( x ) , Q T ( x ) i w h Q ( x ) , Q T ( x ) i − w , where K = h Q ( x ) , Q T ( x ) i w = Z ba Q ( x ) Q ( x ) T w ( x ) dx, (4)is an ( N + 1)( N + 1) matrix and is said dual matrix of Q ( x ).2.2 Convergence of the methodThe following theorem shows that by increasing N , the approximation solu-tion f N ( x ) is convergent to f ( x ) exponentially. Kourosh Parand et al.
Theorem 1. ( Taylor’s formula ) Suppose that f ( x ) ∈ C [0 , b ] and D k f ( x ) ∈ C [0 , b ], where k = 0 , , ..., N . Then we have f ( x ) = N X i =0 x i i ! D i f (0 + ) + x N +1 ( N + 1)! D N +1 f ( ξ ) , (5)with 0 ξ x, ∀ x ∈ [0 , b ]. And thus | f ( x ) − N X i =0 x i i ! D i f (0 + ) | M x N +1 ( N + 1)! , (6)where M > | D N +1 f ( ξ ) | . Proof:
See Ref. [55].
Theorem 2
Suppose that D k f ( x ) ∈ C [0 , b ] for k = 0 , , ..., N, and Q N isthe subspace generated by { Q ( x ) , Q ( x ) , ..., Q N ( x ) } . If f N ( x ) = F T Q ( x ) isthe best approximation to f from Q N , then the error bound is presented asfollows k f ( x ) − f N ( x ) k w M ( N + 1)! (cid:0) b N +2 (2 N + 3) (cid:1) , (7)where M > | D N +1 f ( x ) | , x ∈ [0 , b ] . Proof.
By theorem 1, we have y ( x ) = P Ni =0 x i i ! D i f (0 + ) and | f ( x ) − y ( x ) | M x N +1 N + 1! . Since the best approximation to f ( x ) from Q N is F T Q ( x ), and y ( x ) ∈ Q N ,thus k f ( x ) − f N ( x ) k w k f ( x ) − y ( x ) k w M ( N + 1)! Z b x N +2 w ( x ) dx, where the weight function is w ( x ) = b , thus by integration of above equation,Eq. (7) can be proved (cid:4) . Now, we can show Q ( x ) as metrix product of Y and X s ( x ) as follows: Q ( x ) = Y X s ( x ) , (8)where Q ( x ) = [ Q ( x ) , Q ( x ) , ..., Q N ( x )] T and X s ( x ) = [1 , (cid:0) x − ab − a (cid:1) , (cid:0) x − ab − a (cid:1) , ..., (cid:0) x − ax − b (cid:1) N ] T .There are two different forms for Y . If N is odd [47,48,49,50,51,52]: au method for the numerical solution of non-linear problems 7 Y = − · · · ( − N − ( N − )!( N − )!2 N − · · · ( − N − ( N − )!( N +12 )!2 N · · · ( − N − ( N − )!( N +12 )!2 N − · · · ... ...0 0 0 · · · N − N −
00 0 0 · · · N !2 N And if N is even: Y = − · · · ( − N ( N )!( N )!2 N · · · ( − N − ( N − )!( N )!2 N −
00 0 · · · ( − N − ( N − )!( N +22 )!2 N ... ... ... · · · ... ...0 0 0 · · · N − N −
00 0 0 · · · N !2 N In this step, X s ( x ) as metrix product of S and X ( x ) are expaned as follows: X s ( x ) = S X ( x ) , (9)which X ( x ) = [1 , x, ..., x N ] T and S = · · · − ab − a b − a · · · a ( b − a ) − a ( b − a ) b − a ) · · · (cid:18) N (cid:19)(cid:0) b − a (cid:1) N (cid:0) − a (cid:1) − N (cid:18) N (cid:19)(cid:0) b − a (cid:1) N (cid:0) − a (cid:1) − N +1 (cid:18) N (cid:19)(cid:0) b − a (cid:1) N (cid:0) − a (cid:1) − N +2 · · · (cid:0) b − a (cid:1) N Therfore, Q ( x ) will be rewritten as follows: Q ( x ) = M X ( x ) , (10)where M = YS .We know a triangular matrix is invertible if and only if all its diagonal entriesare nonzero, and the product of two invertible matrices also invertible, so X ( x ) = M − Q ( x ) . (11)This shows that we can represent X ( x ) as metrix product of M − and Q ( x ). Kourosh Parand et al. X ( x ) vector as follows: ddx X ( x ) = P X ( x ) , (12)where P i,j = (cid:26) i, i − j = 10 , otherwise f or i, j = 0 , , ..., N , The matrix D is named as the (N + 1) × (N + 1) operational matrix of deriva-tive if and only if ddx Q ( x ) = D Q ( x ) . (13) Lemma 1.
Suppose that M , M − and D are (N + 1) × (N + 1) matriceswhich satisfy respectively Eqs. (10), (11) and (13). Then D = M P M − . (14) Proof:
By using expansion of ddx Q ( x ), we have ddx Q ( x ) = ddx M X ( x ) = M ddx X ( x ) = M P X ( x ) = M P M − Q ( x ) ⇒ D = M P M − (cid:4) . Remark: If D be a first order operational matrix of drivative of Q ( x ) then k th order operational matrix of drivative of Q ( x ) is as follows: d k dx k Q ( x ) = D k Q ( x ) . (15)3.2 Operational matrix of integralThe integral operational matrix for basis vector X ( x ) operates as: Z xa X ( t ) dt ≃ LX ( x ) , (16)where L i,j = i +1 j − i = 1 and i = N − a i +1 i +1 j = 0 and i = Nr j i = N otherwise f or i, j = 0 , , ..., N , where R = [ r , r , ..., r m ] T .We just need to approximate R xa t N dt = x N +1 N +1 − a N +1 N +1 . au method for the numerical solution of non-linear problems 9 According to the relation (3), we have x N +1 N + 1 − a N +1 N + 1 ≃ N X n =0 a n Q n ( x ) = R T Q ( x ) , (17)thus R T ≃ h x N +1 N + 1 , Q T ( x ) i w h Q ( x ) , Q T ( x ) i − w , f or w ( x ) = 1 b − a . The matrix I is named as the (N + 1) × (N + 1) operational matrix of thederivative of the vector Q(x) if and only if Z x Q ( t ) dt ≃ I Q ( x ) . (18) Lemma 2.
Suppose that M , M − and L are the (N + 1) × (N + 1) matriceswhich satisfy respectively Eqs (10), (11) and (18), then I = M L M − . (19) Proof:
By using expansion of R xa Q ( t ) dt, we have Z x Q ( t ) dt = Z x M X ( t ) dt = M Z x X ( t ) dt ≃ M LX ( x ) = M L M − Q ( x ) ⇒ I = M L M − (cid:4) . X ( x ) vectors will also be appliedas follows: X ( x ) X ( x ) T V ≃ e V X ( x ) T , (20)where V = [ v , v , ..., v N ]. According to the relation (3), we have e V ≃ h X ( x ) , X T ( x ) i w h X ( x ) , X T ( x ) i − w , f or w = 1 b − a . The matrix e C is named as (N + 1) × (N + 1) operational matrix of product ofQ(x) if and only if Q ( x ) Q ( x ) T C ≃ e C Q ( x ) T , where C = [ c , c , ..., c N ]. Lemma 3.
Suppose that M , M − and e V are (N + 1) × (N + 1) matriceswhich satisfy respectively Eqs. (10), (11) and (20). Then e C = M e V M − . (21) Proof:
Assume that V = M T C so by using expansion of Q ( x ) Q ( x ) T C, wehave Q ( x ) Q ( x ) T C = M X ( x ) M X ( x ) C = M X ( x ) X ( x ) T M T C = M X ( x ) X ( x ) V ≃ M e V X ( x ) = M e V M − Q ( x ) ⇒ e C = M e V M − (cid:4) . In the present paper, to illustrate the effectiveness of Tau method, several testexamples are carried out. A comparison of our results with those obtainedby other methods reveals that our methods are very effective and convenient.Now we will demonstrate the use of the SBT method in several engineeringproblems.
Example 1:
We consider the governing equation for unsteady two-dimensionalflow and heat transfer of a viscous fluid [56] which the equation is convertibleto the following system of the nonlinear equations: f ′′′′ − S ( xf ′′′ + 3 f ′′ − f f ′′ ) − G f ′′ = 0 , (22) θ ′′ + S P r (2 f θ ′ − xθ ′ ) + P r Ec ( f ′′ + 12 δ f ′ ) = 0 , (23)with the boundary conditions f (0) = A, f ′ (0) = 0 , θ (0) = 1 ,f (1) = 12 , f ′ (1) = 0 , θ (1) = 0 . (24)To solve this example, we approximate f ( x ) and θ ( x ) functions by N + 1 basisof the shifted Bessel polynomials on the interval [0 ,
1] as bellow: f ( x ) ≃ N X i =0 a i Q i ( x ) = F T Q ( x ) , (25) θ ( x ) ≃ N X i =0 c i Q i ( x ) = C T Q ( x ) , (26)where F = [ a , a , ..., a N ] T and C = [ c , c , ..., c N ] T .By using Eqs. (13) and (20) we have: f ( j ) ( x ) = F T D j Q ( x ) ,xf ( j ) = ( U T Q ( x ))( F T D j Q ( x )) = U T Q ( x ) Q ( x ) T ( F T D j ) T | {z } Z ≃ U T ˜ ZQ ( x ) f f ′′ ( x ) = F T BF T D J ( X ) = F T Q ( x ) Q ( x ) T ( x ) ( D ) T F | {z } O ≃ F T ˜ OQ ( x ) ,f ′ ( x ) = ( F T D Q ( x )) = F T D Q ( x ) Q ( x ) T D T F | {z } W ≃ F T D ˜ W Q ( x ) ,f ′′ ( x ) = ( F T D Q ( x )) = F T D Q ( x ) Q ( x ) T ( x ) ( D ) T F | {z } E ≃ F T D ˜ EQ ( x ) ,θ ( j ) ( x ) = C T D j Q ( x ) ,xθ ′ ( x ) = ( U T Q ( x ))( C T D Q ( x )) = U T Q ( x ) Q ( x ) T ( x ) D T C | {z } V ≃ U T ˜ V Q ( x ) ,f θ ′ = F T Q ( x ) C T D Q ( x ) = F T x T D T C | {z } H ≃ F T ˜ HQ ( x ) , au method for the numerical solution of non-linear problems 11 where Q ( x ) Q ( x ) T Z ≃ ˜ ZQ ( x ), Q ( x ) Q ( x ) T O ≃ ˜ OQ ( x ), Q ( x ) Q ( x ) T E ≃ ˜ EQ ( x ), Q ( x ) Q ( x ) T V ≃ ˜ V Q ( x ) , Q ( x ) Q ( x ) T H ≃ ˜ HQ ( x ), and D j is the jth power ofthe matrix D and we also have: x = X ( x ) T u = Q ( x ) T ( M − ) T u | {z } U = Q ( x ) T U, where u = [0 , , , , ..., T . Therefore, Eqs. (22) and (23) can be rewritten as: F T D Q ( x ) − S (cid:16) U T ˜ ZQ ( x ) + 3 F T D Q ( x ) − F T ˜ OQ ( x ) (cid:17) − G F T D Q ( x ) = 0 , (27) C T D j Q ( x ) + S P r (cid:16) F T ˜ HQ ( x ) − U T ˜ V Q ( x ) (cid:17) − P r Ec (cid:16) F T D ˜ EQ ( x ) + 12 δ F T D ˜ W Q ( x ) (cid:17) = 0 , (28) Now, by using Tau method [57] and Eq. (4), we have: Z (cid:18) F T D − S (cid:16) U T ˜ Z + 3 F T D − F T ˜ O (cid:17) − G F T D (cid:19) Q ( x ) Q ( x ) T dx = (cid:18) F T D − S (cid:16) U T ˜ Z + 3 F T D − F T ˜ O (cid:17) − G F T D (cid:19) Z Q ( x ) Q ( x ) T dx = (cid:18) F T D − S (cid:16) U T ˜ Z + 3 F T D − F T ˜ O (cid:17) − G F T D (cid:19) K = 0 , (29) Z (cid:18) C T D j + S P r (cid:16) F T ˜ H − U T ˜ V (cid:17) − P r Ec (cid:16) F T D ˜ E + 12 δ F T D ˜ W (cid:17)(cid:19) Q ( x ) Q ( x ) T dx = (cid:18) C T D j + S P r (cid:16) F T ˜ H − U T ˜ V (cid:17) − P r Ec (cid:16) F T D ˜ E + 12 δ F T D ˜ W (cid:17)(cid:19) K = 0 . (30)Since K is an invertible matrix, thus one has: F T D − S (cid:16) U T ˜ Z + 3 F T D − F T ˜ O (cid:17) − G F T D = 0 , (31) C T D j + S P r (cid:16) F T ˜ H − U T ˜ V (cid:17) − P r Ec (cid:16) F T D ˜ E + 12 δ F T D ˜ W (cid:17) = 0 , (32)In this method, to satisfy the boundary conditions Eq. (24), we have substi-tuted the four following equations in four rows of Eq. (31): f (0) = F T Q (0) = A, f ′ (0) = F T D Q (0) = 0 ,f (1) = F T Q (1) = 12 , f ′ (1) = F T D Q (1) = 0 , and we have also substituted two following equations in two rows of Eq. (32): θ (0) = C T Q (0) = 1 ,θ (1) = C T D Q (1) = 0 , Finally, we generate 2 N + 2 algebraic equations, therefore by solving aboveequations the unknown vectors F and C are achieved.Effects of deformation parameter S and porosity parameter A are displayed inFigs. (1) and (2). Fig. (3) shows the influences of flow parameters (in suction case) on the tem-perature profile. Fig. (4) shows the influences of physical parameters on veloc-ity and temperature distributions in the case of injection.Fig.(5) shows residual error by the SBT method with N = 15 on f ′ ( x ) and θ ( x ) for the suction cases.Tables (1) and (2) show a comparison between the given values f ′ ( x ) and θ ( x )by the variational iteration method (VIM)[56] and the SBT methods with N = 15 for A = 0 . , S = 0 . , M = 0 . , P r = 0 . , Ec = 0 . δ = 0 .
1. Itis evident from the table that the present method is efficiency and accuracyof VIM method. Table 3 displays the numerical values of Nusselt number ofdifferent values of flow parameters.
Example 2:
In this example, we consider the following Lane-Emden typeequation on the interval [0 ,
3] as follows [50]: y ′′ ( x ) + 2 x y ′ ( x ) − x + 3) y ( x ) = 0 , x , (33)with the initial conditions y (0) = 1 , y ′ (0) = 0 , (34)The exact solution of this problem is y ( x ) = e x . (35)To solve this example, we approximate y ′′ ( x ) functions by N + 1 basis of theBessel polynomials on the interval [0 ,
3] as bellows: y ′′ ( x ) ≃ N X i =0 a i Q i ( x ) = A T Q ( x ) (36)where A = [ a , a , ..., a N ] T . By using Eqs. (13) and (20) we have: x = X ( x ) T w = Q ( x ) T ( M − ) T w | {z } U = Q ( x ) T U, (2 x + 3 x ) = X ( x ) T r = Q ( x ) T ( M − ) T r | {z } Z = Q ( x ) T Z, X ( x ) T p = Q ( x ) T ( M − ) T p | {z } V = Q ( x ) T V, For approximating the solution y ′ ( x ) and y ( x ), we apply one and two-timeintegration on both sides Eq. (48): au method for the numerical solution of non-linear problems 13 Z x y ′′ ( x ) dx = y ′ ( x ) − y ′ (0) ≃ Z x A T Q ( x ) dx ⇒ y ′ ( x ) ≃ A T I Q ( x ) ,y ′′ ( x ) x = A T Q ( x ) Q ( x ) T U ≃ A T ˜ U Q ( x ) , Z x y ′ ( x ) dx = y ( x ) − y (0) ≃ Z x A T I Q ( x ) dx − ⇒ y ( x ) = A T I Q ( x ) + V T Q ( x ) = ( A T I + V T ) | {z } L Q ( x ) = L T Q ( x ) ,y ( x )(2 x + 3 x ) = L T Q ( x ) Q ( x ) T Z ≃ L T ˜ ZQ ( x ) , where w = [0 , , , , ..., T , r = [0 , , , , , ..., T and p = [1 , , , ..., T .Therefore, Eq. (33) can be rewritten as: A T ˜ U Q ( x ) + 2 A T I Q ( x ) − L T ˜ ZQ ( x ) = 0 , (37)As in a typical Tau method [57] , we have: A T ˜ U + 2 A T I − L T ˜ Z = 0 . (38)Finally, we generate N + 1 algebraic equations, therefore by solving aboveequations the unknown vector A is achieved.Table (4) reports the comparison of y ( x ) given by exact value, the Hermitefunctions collocation (HFC) method applied by Parand et al. [61] and theproposed method in this paper with N = 40 . The absolute error graphs ofLane-Emden equation obtained by present method for N = 20 ,
30 and 40 areshown in Figure (6).
Example 3:
We consider the following nonlinear Abel equation on the inter-val [0 ,
1] as follows: y ′ ( x ) = sin ( x ) y ( x ) − xy ( x ) + x y ( x ) − x , x , (39)with the boundary condition y (0) = 0 . (40)To solve this example, we approximate y ′ ( x ) functions by N + 1 basis ofshifted Bessel polynomials on the interval [0 ,
1] as bellows: y ′ ( x ) ≃ N X i =0 a i Q i ( x ) = A T Q ( x ) , (41) where A = [ a , a , ..., a N ] T . By using Eqs. (13) and (20) we have: x = X ( x ) T w = Q ( x ) T ( M − ) T w | {z } U = Q ( x ) T U,x = X ( x ) T r = Q ( x ) T ( M − ) T r | {z } Z = Q ( x ) T Z,x = X ( x ) T p = Q ( x ) T ( M − ) T p | {z } V = Q ( x ) T V,sin ( x ) = Q ( x ) T S,y ( x ) = A T I Q ( x ) ,y ( x ) x = (cid:16) A T I Q ( x ) Q ( x ) T I T A |{z} W (cid:17) Q ( x ) T U ≃ A T I ˜ W ˜ U Q ( x ) ,y ( x ) x = A T I Q ( x ) Q ( x ) T Z ≃ A T I ˜ ZQ ( x ) ,y ( x ) sin ( x ) = (cid:16) A T I ˜ W Q ( x ) Q ( x ) T I T A |{z} W (cid:17) Q ( x ) T S ≃ A T I ˜ W ˜ SQ ( x ) , where w = [0 , , , , ..., T , r = [0 , , , , ..., T , p = [0 , , , , , ..., T and S T = h sin ( x ) , Q T ( x ) i w h Q ( x ) , Q T ( x ) i − w , .Therefore, Eq. (39) can be rewritten as: A T Q ( x ) − A T I ˜ W ˜ SQ ( x ) + A T I ˜ ZQ ( x ) − A T I ˜ ZQ ( x ) + V T Q ( x ) = 0 , (42)As in a typical Tau method [57] , we have: A T − A T I ˜ W ˜ S + A T I ˜ Z − A T I ˜ Z + V T = 0 , (43)In this step, we generate N + 1 algebraic equations, therefore by solving aboveequations the unknown vector A is achieved.The residual error graphs of the nonlinear Abel equation obtained by presentmethod for N = 5 and 10 are shown in Figure (7). Table (5) displays givenvalues y ( x ) by the fractional order of the Chebyshev orthogonal functions(FCFs) collocation method used by parand et al. [58], the present method andthe residual error with N = 10. Example 4:
Consider the nonlinear the standard Lane-Emden problem y ′′ ( x ) + 1 x y ′ ( x ) + y = 0 , x y (0) = 1 , y ′ (0) = 0 . (45)Thus, the solution of the Lane-Emden equation is obtained by the SBTmethod with N = 8 and N = 12.The residual error graphs of the nonlinear standard Lane-Emden obtainedby the present method for N = 8 and 12 are shown in Figure (8). Table (6)reports the comparison of the obtained values y ( x ) by the present methodwith N = 12, Horedt [59], the Bessel orthogonal functions collocation (BFC) au method for the numerical solution of non-linear problems 15 method utilized by Parand et al. [45] and the Chebyshev orthogonal functions(GFCFs) collocation method proposed by Parand et al. [60]. Example 5:
Finally, consider the boundary-value the Troesch’s problem y ′′ ( x ) = γ sinh( γy ( x )) , x , (46)with the boundary conditions y (0) = 0 , y (1) = 1 , (47)We approximate y ( x ) functions by N + 1 basis of shifted Bessel polynomialson the interval [0 ,
1] as bellows: y ( x ) ≃ N X i =0 a i Q i ( x ) = A T Q ( x )where A = [ a , a , ..., a N ] T .By using the SBT technique described in the paper, we can generate a setof linear algebraic equations as follows: A T D − γ ( γA T + γ A T ˜ A + γ A T ˜ A ) = 0 , (48)where Q ( x ) Q ( x ) T A ≃ ˜ AQ ( x ).To satisfy the boundary conditions Eq. (47), we have substituted the fol-lowing equations in two rows of Eq. (48): y (0) = A T Q (0) = 0 , y (1) = A T Q (0) = 1 . (49)Table (7) displays the comparison of the obtained values y ( x ) by the modi-fied non-linear shooting method (MNLSM) applied by Alias [62], the numericalscheme based on the modified homotopy perturbation technique used by [63]and the present method with N = 10 and γ = 0 . The shifted Bessel Tau (SBT) method is a useful technique for solving ordi-nary differential equations. Numerical programs using this method are oftenconsiderably faster with greater accuracy than other standard techniques suchas the variational iteration method, shooting method and homotopy perturba-tion technique. Interest in the Tau method, for a long time regarded only as atool for the construction of accurate approximations of a very restricted classof functions, has been enhanced. In this investigation, the Bessel polynomialsoperational matrices of integration, differentiation, and product are derived.A general procedure for the formulation of these matrices is given. By utilizingthe operational matrices and the Tau method, we reduce the solution of theoriginal equations to the solution of a nonlinear algebraic equation system.Therefore, all of these equations can be solved by Newton method for theunknown coefficients. Some problems with initial or boundary conditions are considered in one dimensional contexts, although the proposed technique canbe easily implemented for two and three dimensional problems. The main ad-vantage of the present method is its simplicity and convenience for computeralgorithms. Illustrative examples demonstrate the validity and applicability ofthe present method.
References
1. B.A. Finlayson, The Method of Weighted Residuals and Variational Principles. AcademicPress, New York (1972).2. S. Kazem, M. Shaban, J.A. Rad, Solution of Coupled Burgers equation based on opera-tional matrices of d-dimensional orthogonal functions, Z. Naturforsch. A, 67 (5), 267-274.3. K. Parand, M. Dehghan, F. Baharifard, Solving a laminar boundary layer equation withthe rational gegenbauer functions. Appl. Math. Model., 37 (2013) 851-863 .4. E. Doha, A. Bhrawy, Efficient spectral-galerkin algorithms for direct solution for second-order differential equations using jacobi polynomials. Numer. Algorithm, 42 (2006) 137-164.5. S. Shateyi, S.S. Motsa, Y. Khan, A new piecewise spectral homotopy analysisof theMichaelis-Menten enzymatic reactions model. Numer. Algorithm, 66 (2014) 495-510.6. C. Lanczos, Trigonometric interpolation of empirical and analytical functions. J. Math.Phys., 17 (1938) 123-199.7. A. Saadatmandi, M. Dehghan, A Tau approach for solution of the space fractional diffu-sion equation. Comput. Math. Appl., 62 (2011) 1135-1142.8. A. Yildirim, S. Sezer, Y. Kaplan, Analytical approach to Boussinesq equation with spaceand timefractional derivatives, Int. J. Numer. Meth. Fluids, 66 (2011) 1315-1324.9. J. H. He, Application of homotopy perturbation method to nonlinear wave equations,Chaos Soliton Fract, 26 (2005) 695-700.10. S. Kazem, J. Rad, K. Parand, Radial basis functions methods for solving Fokker-Planckequation, Eng. Anal. Bound. Elem., 36 (2012)181-189.11. K. Parand, M. Shahini, M. Dehghan, Rational Legendre pseudospectral approach forsolving nonlinear differential equations of Lane-Emden type, J. Comput. Phys., 228 (2009)8830-8840.12. S. Abbasbandy, E. Shivanian, A new analytical technique to solve Fredholms integralequations, Numer. Algorithms, 56 (2011) 27-43.13. S. J. Liao, Series solution of nonlinear eigenvalue problems by means of the homotopyanalysis method, Nonlinear Anal. Real, 10 (2009) 2455-2470.14. K. Parand, M. Razzaghi, Rational Legendre approximation for solving some physicalproblems on semi-infinite intervals, Phys. Scripta, 69 (2004) 353-357.15. K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving higher-order or-dinary differential equations, Int. J. Comput. Math., 81 (2004) 73-80.16. S. Yousefi, M. Dehghan, The use of Hes variational iteration method for solving varia-tional problems, Int. J. Comput. Math., 87 (2010) 1299-1314.17. A. Wazwaz, A reliable treatment of singular Emden-Fowler initial value problems andboundary value problems, Appl. Math. Comput., 217 (2011) 10387-10395.18. A. Saadatmandi, M. Dehghan, A new operational matrix for solving fractional-orderdifferential equations, Comput. Math. Appl., 59 (2010) 1326-1336.19. J. S. Gu, W. S. Jiang, The haar wavelets operational matrix of integration, Int. J. Syst.Sci., 27 (1996) 623-628.20. I. Horng, J. Chou, Shifted Chebyshev direct method for solving variational problems,Int. J. Syst. Sci., 16 (1985) 855-861.21. M. Razzaghi, S. Yousefi, The legendre Wavelets operational matrix of integration, Int.J. Syst. Sci., 32 (2001) 495-502.22. M. Razzaghi, S. Yousefi, Sine-Cosine Wavelets operational matrix of integration and itsapplications in the calculus of variations, Int. J. Syst. Sci., 33 (2002) 805-810.23. S. A. Yousefi, M. Behroozifar, Operational matrices of bernstein polynomials and theirapplications, Int. J. Syst. Sci., 41 (2010) 709-716.au method for the numerical solution of non-linear problems 1724. M. Razzaghi, M. Razzaghi, Fourier series direct method for variational problems, Int.J. Control, 48 (1988) 887-895.25. M. Razzaghi, M. Razzaghi, Taylor series analysis of time-varying multi-delay systems,Int. J. Control, 50 (1989) 183-192.26. M. Razzaghi, M. Razzaghi, Shifted-Jacobi series direct method for variational problems,Int. J. Syst. Sci., 20 (1989) 1119-1129.27. K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving Volterras popu-lation model, Appl. Math. Comput., 149 (2004) 893-900.28. K. Parand, M. Razzaghi, Rational Chebyshev tau method for solving higher-order or-dinary differential equations, Int. J. Comput. Math., 81 (2004) 73-80.29. K. Parand, H. Yousefi, M. Delkhosh, A. Ghaderi,A novel numerical technique to obtainan accurate solution to the Thomas-Fermi equation, Eur. Phys. J. Plus, 131 (7), 1-16.30. K Parand, P Mazaheri, M Delkhosh, A Ghaderi, New numerical solutions forsolving Kidder equation by using the rational Jacobi functions, SeMA J., (2016)doi:10.1007/s40324-016-0103-z.31. E. H. Doha, A. H. Bhrawy and S. S. Ezz-Eldien, A new Jacobi operational matrix:An application for solving fractional differential equations, Appl. Math. Model., 36 (2012)4931-4943.32. S. A. Yousefi, H. Jafari and M. A. Firoozjaeeand S. Momani and C. M. Khaliqued,Application of Legendre wavelets for solving fractional differential equations, Comput.Math. Appl., 62 (2011) 1038-1045.33. S. A. Yousefi, F. Khellat,The linear Legendre mother wavelets operational matrix ofintegration and its application, J. Franklin Inst., 343 (2006) 181-190.34. S. A. Yousefi, M. Behroozifar, M. Dehghan, Numerical solution of the nonlinear age-structured population models by using the operational matrices of Bernstein polynomials,Appl. Math. Modell., 36 (2012) 945-963.35. S. A. Yousefi, M. Behroozifar, M. Dehghan, The operational matrices of Bernstein poly-nomials for solving the parabolic equation subject to specification of the mass, J. Comput.Appl. Math., 335 (2011) 5272-5283.36. A. H. Bhrawy, The operational matrix of fractional integration for shifted Chebyshevpolynomials, Appl. Math. Lett., 26 (2013) 25-31.37. H. R. Marzban, M. Razzaghi, Optimal control of linear delay systems via hybrid ofblock-pulse and Legendre polynomials, J. Franklin Inst., 341 (2004) 279-293.38. H. R. Marzban, M. Shahsiah, Solution of piecewise constant delay systems using hybridof block-pulse and Chebyshev polynomials, Optim. Contr. Appl. Met., 32 (2011) 647-659.39. M. Lakestani, M. Dehghan, S. I. Pakchin, The construction of operational matrix offractional derivatives using B-spline functions, Commun. Nonlinear. Sci. Numer. Simulat.,17 (2012) 1149-1162.40. S. Kazem, M. Shaban, J. A. Rad, Solution of the coupled burgers equation based onoperational matrices of d-dimensional orthogonal functions, Z. Naturforsch. A., 67 (2012)267-274.41. S. Nemati, P. M. Lima, Y. Ordokhani, Numerical solution of a class of two-dimensionalnonlinear Volterra integral equations using Legendre polynomials, J. Comput. Appl.Math., 242 (2013) 53-69.42. K. Parand, M.i Delkhosh, Solving Volterras population growth model of arbitrary orderusing the generalized fractional order of the Chebyshev functions, Ricerch. Mat., 65 (2016)307-328.43. W.W. Bell, Special functions for scientists and engineers, D. Van Nostrand Company,CEf Canada, 1967.44. K. Parand, M. Nikarya, J.A. Rad, F. Baharifard, A new reliable numerical algorithmbased on the first kind of Bessel functions to solve Prandtl-Blasius laminar viscous flowover a semi-infinite flat plate, Z. Naturforsch. A, 67 (2012) 665-673.45. K. Parand, M. Nikarya, J.A. Rad, Solving non-linear Lane-Emden type equations usingBessel orthogonal functions collocation method, Celest. Mech. Dyn. Astr., 116 (2013) 97-107.46. M. Delkhosh, The conversion a Bessel’s equation to a self-adjoint equation and appli-cations, World Appl. Sci. J., 15, 1687-1691 (2011).8 Kourosh Parand et al.47. N. Sahin, S. Yuzbasi, M. Gulsu, A collocation approach for solving systems of linearVolterra integral equations with variable coefficients, Comput. Math. Appl., 62 (2011)755-769.48. S. Yuzbasi, M. Sezer, A numerical method to solve a class of linear integro-differentialequations with weakly singular kernel, Math. Method Appl. Sci., 35 (2012) 621-632.49. S. Yuzbasi, N. Sahin, M. Sezer, Bessel polynomial solutions of high-order linear Volterraintegro-differential equations, Comput. Math. Appl., 62 (2011) 1940-1956.50. S. Yuzbasi, A numerical approach for solving a class of the nonlinear Lane-Emden typeequations arising in astrophysics, Math. Method Appl. Sci., 34 (2011) 2218-2230.51. S. Yuzbasi, Bessel collocation approach for solving continuous population models forsingle and interacting species, Appl. Math. Model., 36 (2012) 3787-3802.52. S. Yuzbas, A collocation method based on the Bessel functions of the first kind forsingular perturbated differential equations and residual correction, Math. Method Appl.Sci., 38 (2015) 3033-3042.53. E. Tohidi, H. Nik, A Bessel collocation method for solving fractional optimal controlproblems, Appl. Math. Model., 39 (2015) 455-465.54. K. Parand, A. Ghaderi, H. Yousefi, M. Delkhosh, A new approach for solving nonlinearThomas-Fermi equation based on fractional order of rational Bessel functions, Electron.J. Differential Equations, 2016 (2016), No. 331, pp. 1-18.55. Z. Odibat, S. Momani, An algorithm for the numerical solution of differential equationsof fractional order, J. Appl. Math. Info., 26 (2008) 15-27.56. S.I. Khan, N. Ahmed, U. Khan, S. Ullah Jan, S.T. Mohyud-Din, Heat transfer analysisfor squeezing flow between parallel disks, J. Egypt. Math. Society, 23 (2015) 445-450.57. D. Gottlieb, M. Hussaini, S. Orszg, Theory and applications of spectral methods inspectral methods for partial differential equations, SIAM, Philadelphia, 1984.58. K. Parand, M. Delkhosh, M. Nikarya, Novel Orthogonal Functions for Solving Differ-ential Equations of Arbitrary Order, Tbilisi Math. J., 10 (2017) 31-55.59. G. P. Horedt, Polytropes: Applications in Astrophysics and Related Fileds. Kluwer,Dordecht (2004).60. K. Parand, M. Delkhosh, An effective numerical method for solving the nonlinear sin-gular Lane-Emden type equations of various orders, J. Teknologi, 79 (2017) 25-36 .61. K. Parand, M. Dehghan, A. R. Rezaei, S. M. Ghaderi, An approximation algorithmfor the solution of the nonlinear Lane-Emden type equations arising in astrophysics usingHermite functions collocation method, Comp. Phys. Commun., 181 (2010) 1096-1108.62. N. Alias, A. Manaf, A. Ali, M. Habib, Solving Troeschs problem by using modifiednonlinear Shooting method, J. Teknologi, 78 (2016) 45-52.63. X. Feng, L. Mei, G. He, An efficient algorithm for solving Troesch’s problem, Appl.Math. Comp., 189 (2007) 500-507.(a) Effect of A on f ′ ( x ) (b) Effect of S on f ′ ( x ) Fig. 1
Effect of A and S on f ′ ( x ) for the suction cases, Example 1.au method for the numerical solution of non-linear problems 19(a) Effect of A on f ′ ( x ) (b) Effect of S on f ′ ( x ) Fig. 2
Effect of A and S on f ′ ( x ) for the injection cases, Example 1.(a) Effect of Pr on θ ( x ) (b) Effect of A on θ ( x ) Fig. 3
Effect of Pr and A on θ ( x ) for the suction cases, Example 1.(a) Effect of Pr on θ ( x ) (b) Effect of S on θ ( x ) Fig. 4
Effect of Pr and S on θ ( x ) for the injection cases, Example 1.0 Kourosh Parand et al.(a) Effect of A on f ′ ( x ) (b) Effect of A on θ ( x ) Fig. 5
Logarithm graphs of residual error by effects of A on f ′ ( x ) and θ ( x ) for the suctioncases with N = 15, Example 1. Table 1
Comparison of the obtained values f ′ ( x ) by the variational iteration method(VIM)[56] and the present method, Example 1. x Bessel Tau method VIM[56]
Res ( x )0.2 0.384801280290557 0.384801 1.39535e-120.4 0.575554143054330 0.575554 8.21911e-140.6 0.575174675564395 0.575174 1.90287e-120.8 0.384040432902588 0.384040 1.91562e-12 Table 2
Comparison of the obtained values θ ( x ) by the variational iteration method(VIM)[56] and the present method, Example 1. x Bessel Tau method VIM[56]
Res ( x )0.2 0.806144282850332 0.806144 8.25970e-140.4 0.607022034290869 0.607022 7.26187e-140.6 0.407003862839112 0.407004 7.24240e-140.8 0.206120555470330 0.206121 8.19659e-14 Table 3
Comparison of values of Nusselt number (1 − at ) / Nu for different values of P r, Ec , and δ when A = 0 . , S = 0 . , M = 0 . N = 15 andthe variational iteration method (VIM) [56], Example 1.Pr Ec δ (1 − at ) / Nu VIM [56]0.0 0.2 0.1 1.00000000000000000 1.000000.1 1.01893685003010788 1.018940.2 1.03787156909761297 1.037870.3 0.0 0.99859957004178043 0.998600.2 1.05680415677248143 1.056800.4 1.11500874350318244 1.115010.3 0.0 1.08487154864359339 1.084870.5 1.11074408599955706 1.110741.0 1.18836169806744806 1.18836au method for the numerical solution of non-linear problems 21
Fig. 6
Logarithm graphs of absolute error function, Example 2.
Table 4
Comparison of the obtained values y ( x ) by exact value, HFC method [61] and thepresent method with N = 40, Example 2. x Exact value HFC [61] Peresent method0.01 1.00010000500016667083 1.0000999826 1.000100005000166657220.02 1.00040008001066773341 1.0004000642 1.000400080010667748810.05 1.00250312760579508497 1.0025031064 1.002503127605795073090.10 1.01005016708416805754 1.0100501492 1.010050167084168055460.20 1.04081077419238822675 1.0408107527 1.040810774192388223990.50 1.28402541668774148407 1.2840253862 1.284025416687741478180.70 1.63231621995537897012 1.6323161777 1.632316219955378962940.80 1.89648087930495135334 1.8964808279 1.896480879304951360370.90 2.24790798667647141917 2.2479078937 2.247907986676471412321.00 2.71828182845904523536 2.7182819166 2.718281828459045241841.5 9.48773583635852572055 — 9.487735836358525723002.0 54.5981500331442390781 — 54.59815003314423907192.5 518.012824668342025939 — 518.0128246683420259473.0 8103.08392757538400770 — 8103.083927575384007652 Kourosh Parand et al.
Fig. 7
Logarithm graphs of residual error function, Example 3.
Table 5
Numerical given values y ( x ) by FCFs collocation method [58] and the presentmethod with N = 10, Example 3. x FCFs [58] Present method Residual error0.1 -2.500e-5 -2.541201381e-5 2.13471e-60.2 -4.004e-4 -4.000957821e-4 1.42107e-60.3 -2.032e-3 -2.033184449e-3 3.48798e-60.4 -6.459e-3 -6.459447204e-3 7.35610e-60.5 -1.591e-2 -1.591376982e-2 6.64911e-60.6 -3.346e-2 -3.346334651e-2 1.36557e-60.7 -6.327e-2 6.327405918e-2 4.86862e-60.8 -1.111e-1 -1.111347323e-1 9.41866e-60.1 -1.855e-1 1.855105463e-1 1.23827e-51.0 -2.999e-1 -2.999554921e-1 1.98917e-4
Table 6
Comparison of the obtained values y ( x ) by the present method with N = 12,Horedt [59], BFC method [45] and GFCFs collocation method [60], Example 4. x Present method Horedt [59] BFC [45] GFCFs [60]0.1 0.998334998549872 0.9983350 0.99833499854 0.998334999860.3 0.985133946938390 – – –0.5 0.959352715810926 0.9593527 0.95935271580 0.959352715850.7 0.922170348514590 – – –1.0 0.848654111411546 0.8486541 0.84865411140 0.848654096031.5 0.695367147241325 – – –2.0 0.529836429310169 – – –au method for the numerical solution of non-linear problems 23
Fig. 8
Logarithm graphs of residual error function, Example 4.
Table 7
Comparison of the obtained values y ( x ) by Alias [62], Feng [63] and the presentmethod with N = 10 and γ = 0 .
5, Example 5. xx