Fourier Operational Matrices of Differentiation and Transmission : Introduction and Applications

and Applied Analysis 3 The vector representation of the above-mentioned truncated series can be written in the following form: y N (t)


Introduction
In the last four decades, numerical methods which are based on the operational matrices of integration (especially for orthogonal polynomials and functions) have received considerable attention for dealing with a huge size of applied mathematics problems such as system identification, state space analysis, optimal control, and senstivity analysis.The key idea of these methods is based on the integral expression where Φ() = [Φ 1 (), Φ 2 (), . . ., Φ  ()] is an arbitrary basis vector and  is a ( + 1) × ( + 1) constant matrix, called the operational matrix of integration.The matrix  has already been determined for many types of orthogonal (or nonorthogonal) bases such as Walsh functions [5][6][7], block-pulse functions [8], Laguerre polynomials [9], Chebyshev polynomials [10], Legendre polynomials [11], Hermite polynomials [12], Fourier series [13], Bernstein polynomials [14], and Bessel series [15].As a primary research work which was based on the operational matrices of integration, one can refer to the work of Corrington [5].In [5], the author proposed a method of solving nonlinear differential and integral equations using a set of Walsh functions as the basis.His method is aimed at obtaining piecewise constant solutions of dynamic equations and requires previously prepared tables of coefficients for integrating Walsh functions.
To alleviate the need for such tables, Chen and Hsiao [6,7] introduced an operational matrix to perform integration of Walsh functions.This operational matrix approach has been applied to various problems such as time-domain analysis and synthesis of linear systems, piecewise constant-feedbackgain determination for optimal control of linear systems and for inverting irrational Laplace transforms.On the other hand, since the beginning of 1994, the Bernoulli, Chebyshev, Laguerre, Bernstein, Legendre, Taylor, Hermite, and Bessel matrix methods have been used in the works [4,[16][17][18][19][20][21][22][23][24][25][26][27][28] to solve high-order linear and nonlinear Fredholm Volterra integro-differential difference equations and their systems.The main characteristic of these approaches is based on the operational matrices of differentiation instead of integration.The best advantage of these techniques with respect to the integration methods is that, in the fundamental matrix relations, there is not any approximation symbol, meanwhile in the integration forms such as (1), the approximation symbol could be seen obviously.In other words, where  is the operational matrix of differentiation for any selected basis such as the above-mentioned polynomials, functions, and truncated series.The readers can see that there is no approximation symbol in (2), meanwhile this can be seen in ( 1) by using operational matrices of integration.For justifying this expression, one can refer to this subject that after the differentiating of an th degree polynomial we usually reach to a polynomial which has less than th degree.However, in the integration processes the degree of polynomials would be increased.
To the best of our knowledge, this is the first work concerning the Fourier matrix method for solving high-order linear ordinary differential equations (ODEs) in a differentiation form of view.This partially motivated our interest in such methods.In this paper, in the light of the methods [4, 16-20, 24, 25] and by means of the matrix relations between the truncated Fourier series and their derivatives, we propose a new method, namely, the Fourier matrix method for solving th order linear ODEs with constant coefficients in the form where () is a known function and the initial conditions  () (−) =   , for  = 0, 1, . . .,  − 1 are given.Also, we present a new method by using matrix relations between the truncated Fourier series and their transmissions for solving the following th order difference equation: Moreover, by means of the matrix relations between the truncated Fourier series and their derivatives together with their transmissions and by using suitable Gaussian collocation nodes, we obtain the numerical solution of generalized Pantograph equation with variable coefficients in the form subject to the initial conditions  () (−) =   , for  = 0, 1, . . .,  − 1.
It should be noted that delay differential equations have a wide range of application in science and engineering.Functional-differential equations with proportional delays are usually refereed as Pantograph equations or generalized Pantograph equations.The applications of Pantograph equation are in different fields such as number theory, economy, biology, control, electrodynamics, nonlinear dynamical systems, quantum mechanics, probability theory, astrophysics, cell growth, and other industrial applications.For some applications of this equation, we refer the interested reader to [21].Pantograph equation was studied by many authors numerically and analytically.A complete survey about numerical and analytical methods about generalized Pantograph equations has been provided in [21].
The rest of the paper is organized as follows.In Section 2, we introduce some mathematical preliminaries of Fourier series together with their operational matrices of differentiation and transmission.Section 3 is devoted to apply the Fourier matrix methods for solving (3), (4), and (5) using the Fourier operational matrices of differentiation and transmission.An error estimation of the collocation scheme for solving the generalized Pantograph equation ( 5) is provided in Section 4. In Section 5, the proposed method is applied to several numerical examples.Also conclusions and future works are given in Section 6.

Some Properties of Fourier Series
Fourier series decomposes periodic functions or periodic signals into the sum of a (possibly infinite) set of simple oscillating functions, namely, sines and cosines.The Fourier series has many applications in electrical engineering, vibration analysis, acoustics, optics, signal processing, image processing, quantum mechanics, econometrics, thin-walled shell theory, and so forth.A square integrable function () can be expanded in terms of Fourier series as follows: where the coefficients   and f are given by In practice, we use only the first (2 + 1) terms of the Fourier series.Therefore, our aim is to approximate the solution of (3), (4), and (5) as follows: Abstract and Applied Analysis 3 The vector representation of the above-mentioned truncated series can be written in the following form: According to the property   () = (), we have    () =   (), where    () denotes the first derivative of   ().By repeating this procedure, we have  ()   () =  () (), where  ()   () denotes the th derivative of   ().On the other hand, the matrix relation between   () and (  ())  can be constructed by using the elementary calculus as follows: where  is the (2+1)×(2+1) Fourier operational matrix of differentiation.The above relation can be written in the form   () = ()  .Accordingly, the th derivative of () can be given by Therefore,  ()  () has the following form: We now construct the operational matrices of transmission in both cases of forward and backward by using the following simple relations from the elementary calculus: Thus, the relation between (( + ))  and   () is as follows: Moreover, the relation between (( − ))  and   () can be written in the following form: In other words, ( + ) = ()(  + )  and ( − ) = ()(  − )  .Therefore, It should be noted that ( 12) and ( 16) are the fundamental matrix relations of this section.In other words, the operational matrices of differentiation  and transmission   + would be applied in our methods in the next section abundantly.

Matrix and Collocation Methods Based on Fourier Operational Matrices
In this section, we use the matrix relations ( 12) and ( 16) for solving (3), (4), and (5).In the solution approximation procedure of (3), we use only the fundamental relation (12), and for approximating the solution of (4), we use only the fundamental relation (16), meanwhile both of relations ( 12) and ( 16) together with the collocation method based on Legendre Gauss points are applied to solve numerically the generalized Pantograph equation (5).

Fourier Matrix Method for Solving High-Order Linear
ODEs.We now derive an algorithm for solving (3).To do this, let the solution of (3) be approximated by the first (2 + 1)-terms of the truncated Fourier series.Hence, if we write where the unknown Fourier coefficient vector  and the Fourier vector () are given in (9), then the th derivative of   () can be expressed in the matrix form by  ()  () =  () () = ()(  )   for  = 0, 1, . . ..Moreover, we need to approximate the given function () which exists in (3) by using (7).In other words, we assume that where is a known vector with the aid of (7).By substituting the approximations ( 12) and ( 18) into (3), we get According to the completeness property of Fourier series, the matrix equation corresponding to the basic equation ( 3) can be written as  =  or in the form of augmented matrix [; ], where Therefore, the basic equation ( 3) is changed into a system of (2+1) linear algebraic equations with unknown coefficients  0 ,  1 , . . .,   , ŷ1 , . . ., ŷ which can be written in the form or in augmented matrix form Now we should impose the initial conditions  () (−) =   for  = 0, 1, . . .,  − 1 to the associated system of algebraic According to (12), we can write  () (−) = (−)(  )   and hence the matrix-vector representation of the initial conditions is as follows: Therefore, the matrix form of the initial conditions can be written as where Finally, from ( 23) and ( 25), (3) subject to the considered initial conditions reduces to the following system of algebraic equations: where The matrix coefficients W of the above system is sparse, and for solving system of linear equations W = F, we can use some iterative Krylov subspace methods and determine the vector  and hence the approximate solution   () = () is obtained.

Fourier Matrix Method for Solving High-Order Linear
Difference Equations.In this subsection, we propose an idea for approximating the solution of (4) by using Fourier operational matrix of transmission.For this purpose, we assume that the solution of ( 4) can be approximated by truncated Fourier series in the form () ≈   () = ().
According to (16), we have ( +   ) ≈   ( +   ) = ( +   ) = ()(   + )   for  = 0, 1, . . ., .Similar to the previous method, the known function () can be approximated by truncated Fourier series in the form () ≈   () = () where the entries of  would be calculated from (7).By substituting these matrix forms in (4), we get Hence, the matrix equation corresponding to the basic equation ( 4) can be written as  =  or in the form of augmented matrix [; ], where Therefore, the basic equation ( 4) is changed into a system of (2+1) linear algebraic equations with unknown coefficients  0 ,  1 , . . .,   , ŷ1 , . . ., ŷ which can be written in the form or in augmented matrix form According to the structure of Fourier operational matrices of transmissions, the matrix coefficients  of the above system is tridiagonal (see for instance ( 14)), and for solving system of linear equations  = , we can use some iterative krylov subspace methods and determine the vector  and hence the approximated solution   () = () is obtained.

Fourier Collocation Method for Solving Generalized Pantograph Equations.
In this subsection, we use the collocation method (via Legendre Gauss points) based on Fourier operational matrices of differentiation and transmission to solve numerically the generalized Pantograph equation (5).
Hence, the matrix equation (36) corresponding to the basic equation ( 5) can be written as  =  or in the form of augmented matrix [; ], where Therefore, the basic equation ( 5) is changed into a system of (2 −  + 1) linear algebraic equations with unknown coefficients  0 ,  1 , . . .,   , ŷ1 , . . ., ŷ , which can be written in the form or in augmented matrix form Now we should impose the initial conditions  () (−) =   to the associated algebraic system of equations in a proper manner.These initial conditions can be written as Therefore, the matrix form of the initial conditions can be written as where Finally, from (40) and ( 42), the generalized Pantograph equation ( 5) subject to the initial conditions  () (−) =   reduces to the following system of algebraic equations: where The matrix coefficients W of the above system may be sparse, and for solving system of linear equations W = F, we can use some iterative krylov subspace methods and determine the vector  and hence the approximated solution   () = () is obtained.
Remark 1.For more information about iterative krylov subspace methods, one can point out to [29].In this book, several iterative methods have been introduced for solving large sparse linear systems.

Error Estimation
In this section, we will give an error estimation (which was previously proposed in the works [27,28]) for the collocation scheme that is proposed for solving Pantograph equation ( 5) with the aid of residual function.This idea may help us to estimate the error especially if the exact solution of (5) does not exist.For this purpose, we can rewrite (5) with the initial conditions  () (−) = 0,  = 0, 1, . . .,  − 1.
Similar to the previous collocation method, we now solve the above error problem and approximate   () by the aid of truncated Fourier series (with the 2 M + 1 terms) in the following form: If the exact solution of ( 5) is unknown, then the absolute errors |  (  )| = |(  ) −   (  )| are not found.However, the absolute errors can be computed approximately with the aid of the estimated absolute error function | , M()|.

Numerical Examples
In this section, several numerical examples are given to illustrate the accuracy and effectiveness of the proposed methods, and all of them are performed on a computer using programs written in MAPLE 13.In this regard, we have reported in the tables and figures the values of the exact solution (), the approximate solution   (), the absolute error function   () = |() −   ()|, and estimate absolute error function  , M() (or the absolute residual functions associated to (3), (4), and ( 5)) at any selected points of the given computational interval.It should be noted that in the first two examples we provide equations in which our results are exact, meanwhile other approaches which were based on the operational matrices of integration and differentiation such as [5,11,14,[16][17][18] could not obtain the exact solutions.
Example 2. As the first example we consider the following ODE with the exact solution () = sin() subject to the initial condition (−) = 0.
According to (8), one can approximate the solution of the above ODE with assumption  = 1 as follows: Our aim is to determine the vector components  with the aid of Fourier operational matrix of differentiation.By using (10) (i.e.,   () = ()  ), we have Moreover, we can approximate the known function () = sin() + cos() by the aid of the truncated Fourier series in the form By substituting the matrix forms of (51)-( 53) into (50) we have Since (−) = 0, then (−) = [1/2 − 1 0] = 0. Thus, the above matrix equation transforms into the following matrix equation: [ The solution of the above system is  = [0 0 1]  and hence which is the exact solution.
According to (8), one can approximate the solution of the above difference equation with assumption  = 1 as follows: Our aim is to determine the vector components  with the aid of Fourier operational matrix of transmission.By using (16) (i.e., ( +   ) = ()( Moreover, we can approximate the known function () = − sin() − cos() by the aid of the truncated Fourier series in the form By substituting the matrix forms of ( 59)-( 61) into (57), we have (65) In Table 1, we provide the absolute values of errors for  = 2 and 3.Moreover, we estimate the mentioned errors by solving the associated error problem which is constructed in Section 4 by using M = 3, 4, and 5.This table shows the agreement of the estimate errors with the actual errors.(66) Since the problem (66) has a complicated exact solution, we compare the results of collocation method based on presented method (PM) with the other four methods in Table 2.These methods that are used to obtain the approximate solutions of this Example are Walsh series approach [1], Laguerre series technique [2], Taylor series method [3], and Hermit series collocation approach [4].Moreover, it seems that the Walsh series approach [1] and Laguerre series technique [2]
equations  =  in a proper way.For this reason and for clarity of presentation, we can remove the last  rows of  and the last  entries of  without loss of generality.However, we can do this work by the first  rows of the  and the first  entries of .Thus, we have 2− ) cos (2   2− ) . . .sin (   2− )