# Fourier Operational Matrices of Differentiation and Transmission: Introduction and Applications

Authors: F. Toutounian; Emran Tohidi; A. Kilicman
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper introduces Fourier operational matrices of differentiation and transmission for solving high-order linear differential and difference equations with constant coefficients. Moreover, we extend our methods for generalized Pantograph equations with variable coefficients by using Legendre Gauss collocation nodes. In the case of numerical solution of Pantograph equation, an error problem is constructed by means of the residual function and this error problem is solved by using the mentioned collocation scheme. When the exact solution of the problem is not known, the absolute errors can be computed approximately by the numerical solution of the error problem. The reliability and efficiency of the presented approaches are demonstrated by several numerical examples, and also the results are compared with different methods.


---

## Body

## 1. 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 assystem identification, state space analysis, optimal control, and senstivity analysis. The key idea of these methods is based on the integral expression
(1)∫0t‍Φ(τ)dτ≈Φ(t)P,
where Φ(t)=[Φ1(t),Φ2(t),…,ΦN(t)] is an arbitrary basis vector and P is a (N+1)×(N+1) constant matrix, called the operational matrix of integration. The matrix P has already been determined for many types of orthogonal (or nonorthogonal) bases such as Walsh functions [5–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-feedback-gain 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–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,
(2)Φ′(t)=Φ(t)B,
where B 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 Nth degree polynomial we usually reach to a polynomial which has less than Nth 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 mth order linear ODEs with constant coefficients in the form
(3)y(m)(t)=∑k=0m-1‍βky(k)(t)+f(t),t∈[-π,π],
where f(t) is a known function and the initial conditions y(k)(-π)=ηk, for k=0,1,…,m-1 are given.Also, we present a new method by using matrix relations between the truncated Fourier series and their transmissions for solving the followingmth order difference equation:
(4)∑k=0m‍γky(t+αk)=f(t),t∈[-π,π].
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
(5)y(m)(t)=∑j=0J‍∑k=0m-1‍pjk(t)y(k)(λjkt+μjk)+g(t),t∈[-π,π],
subject to the initial conditions y(k)(-π)=θk, for k=0,1,…,m-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 Section2, 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.


## 2. 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 functionf(t) can be expanded in terms of Fourier series as follows:
(6)f(t)=f02+∑n=1∞‍(fncos(nt)+f^nsin(nt)),t∈(-π,π),
where the coefficients fn and f^n are given by
(7)fn=1π∫-ππ‍f(t)cos(nt)dt(n≥0),f^n=1π∫-ππ‍f(t)sin(nt)dt(n≥1).
In practice, we use only the first (2N+1) terms of the Fourier series. Therefore, our aim is to approximate the solution of (3), (4), and (5) as follows:
(8)y(t)≈yN(t)=y02+∑n=1N‍(yncos(nt)+y^nsin(nt)),t∈[-π,π],2N+1>m.
The vector representation of the above-mentioned truncated series can be written in the following form:
(9)yN(t)=[12cos(t)cos(2t)⋯cos(Nt)sin(t)sin(2t)⋯sin(Nt)]nnn×[y0y1y2⋮yNy^1y^2⋮y^N]=F(t)Y,t∈(-π,π).
According to the property yN(t)=F(t)Y, we have yN′(t)=F′(t)Y, where yN′(t) denotes the first derivative of yN(t). By repeating this procedure, we have yN(k)(t)=F(k)(t)Y, where yN(k)(t) denotes the kth derivative of yN(t). On the other hand, the matrix relation between FT(t) and (F′(t))T can be constructed by using the elementary calculus as follows:
(10)[12cos(t)cos(2t)⋮cos(Nt)sin(t)sin(2t)⋮sin(Nt)]′︷(F′(t))T=[000⋯000⋯0000⋯0-10⋯0000⋯00-2⋯0⋮⋮⋮⋱⋮⋮⋮⋱⋮000⋯000⋯-N010⋯000⋯0002⋯000⋯0⋮⋮⋮⋱⋮⋮⋮⋱⋮000⋯N00⋯0]︷M[12cos(t)cos(2t)⋮cos(Nt)sin(t)sin(2t)⋮sin(Nt)]︷FT(t),
where M is the (2N+1)×(2N+1) Fourier operational matrix of differentiation. The above relation can be written in the form F′(t)=F(t)MT. Accordingly, the kth derivative of F(t) can be given by
(11)F(k)(t)=F(t)(MT)k,k=0,1,…,m.
Therefore, yN(k)(t) has the following form:
(12)yN(k)(t)=F(k)(t)Y=F(t)(MT)kY,k=0,1,…,m.
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:
(13)cos(nt±α)=cos(nt)cos(α)∓sin(nt)sin(α),sin(nt±α)=sin(nt)cos(α)±cos(nt)cos(α),nnnnnnnnnnnnnnnnnnnnnnnnnnn1≤n≤N.
Thus, the relation between (F(t+α))T and FT(t) is as follows:
(14)[12cos(t+α)cos(2t+α)⋮cos(Nt+α)sin(t+α)sin(2t+α)⋮sin(Nt+α)]︷(F(t+α))T=[100⋯000⋯00cos(α)0⋯0-sin(α)0⋯000cos(α)⋯00-sin(α)⋯0⋮⋮⋮⋱⋮⋮⋮⋱⋮000⋯cos(α)00⋯-sin(α)0sin(α)0⋯0cos(α)0⋯000sin(α)⋯00cos(α)⋯0⋮⋮⋮⋱⋮⋮⋮⋱⋮000⋯sin(α)00⋯cos(α)]︷D+α×[12cos(t)cos(2t)⋮cos(Nt)sin(t)sin(2t)⋮sin(Nt)]︷FT(t).
Moreover, the relation between (F(t-α))T and FT(t) can be written in the following form:
(15)[12cos(t-α)cos(2t-α)⋮cos(Nt-α)sin(t-α)sin(2t-α)⋮sin(Nt-α)]︷(F(t-α))T=[100⋯000⋯00cos(α)0⋯0sin(α)0⋯000cos(α)⋯00sin(α)⋯0⋮⋮⋮⋱⋮⋮⋮⋱⋮000⋯cos(α)00⋯sin(α)0-sin(α)0⋯0cos(α)0⋯000-sin(α)⋯00cos(α)⋯0⋮⋮⋮⋱⋮⋮⋮⋱⋮000⋯-sin(α)00⋯cos(α)]︷D-α×[12cos(t)cos(2t)⋮cos(Nt)sin(t)sin(2t)⋮sin(Nt)]︷FT(t).
In other words, F(t+α)=F(t)(D+α)T and F(t-α)=F(t)(D-α)T. Therefore,
(16)y(t+αk)=F(t+αk)Y=F(t)(D+αk)TY,y(t-αk)=F(t-αk)Y=F(t)(D-αk)TY,0≤k≤m.
It should be noted that (12) and (16) are the fundamental matrix relations of this section. In other words, the operational matrices of differentiation M and transmission D+α would be applied in our methods in the next section abundantly.


## 3. 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).


### 3.1. 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 (2N+1)-terms of the truncated Fourier series. Hence, if we write
(17)y(t)≈yN(t)=y02+∑n=1N‍(yncos(nt)+y^nsin(nt))=F(t)Y,t∈[-π,π],2N+1>m,
where the unknown Fourier coefficient vector Y and the Fourier vector F(t) are given in (9), then the kth derivative of yN(t) can be expressed in the matrix form by yN(k)(t)=F(k)(t)Y=F(t)(MT)kY for k=0,1,…m. Moreover, we need to approximate the given function f(t) which exists in (3) by using (7). In other words, we assume that
(18)f(t)≈fN(t)=f02+∑n=1N‍(fncos(nt)+f^nsin(nt))=F(t)F,t∈[-π,π],
where F=[f0f1f2⋯fNf^1f^2⋯f^N]T is a known vector with the aid of (7). By substituting the approximations (12) and (18) into (3), we get
(19)F(t)(MT)mY=∑k=0m-1‍βkF(t)(MT)kY+F(t)F.
According to the completeness property of Fourier series, the matrix equation corresponding to the basic equation (3) can be written as WY=F or in the form of augmented matrix [W;F], where
(20)W=[wij]={(MT)m-∑k=0m-1‍βk(MT)k}.
Therefore, the basic equation (3) is changed into a system of (2N+1) linear algebraic equations with unknown coefficients y0,y1,…,yN,y^1,…,y^N which can be written in the form
(21)WY=F,
or in augmented matrix form
(22)[W;F]=[w00w01w02…w0(2N);f0w10w11w12…w1(2N);f1w20w21w22…w2(2N);f2⋮⋮⋮⋱⋮;⋮w(2N)0w(2N)1w(2N)2…w(2N)(2N);f^N].
Now we should impose the initial conditions y(k)(-π)=ηk for k=0,1,…,m-1 to the associated system of algebraic equations WY=F in a proper way. For this reason and for clarity of presentation, we can remove the last m rows of W and the last m entries of F without loss of generality. However, we can do this work by the first m rows of the W and the first m entries of F. Thus, we have
(23)[W^;F^]=[w00w01w02…w0(2N);f0w10w11w12…w1(2N);f1w20w21w22…w2(2N);f2⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(2N-m)(2N);f^N-m].
According to (12), we can write y(k)(-π)=F(-π)(MT)kY and hence the matrix-vector representation of the initial conditions is as follows:
(24)F(-π)(MT)kY=[ηk],k=0,1,…,m-1.
Therefore, the matrix form of the initial conditions can be written as
(25)UkY=[ηk]or[Uk;ηk],k=0,1,…,m-1,
where
(26)Uk=F(-π)(MT)k=[uk0,uk1,…,uk(2N)],k=0,1,…,m-1.
Finally, from (23) and (25), (3) subject to the considered initial conditions reduces to the following system of algebraic equations:
(27)W~Y=F~,
where
(28)[W~;F~]=[w00w01w02…w0(2N);f0w10w11w12…w1(2N);f1w20w21w22…w2(2N);f2⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(2N-m)(2N);f^N-mu00u01u02…u0(2N);η0u10u11u12…u1(2N);η1u20u21u22…u2(2N);η2⋮⋮⋮⋱⋮;⋮u(m-1)0u(m-1)1u(m-1)2…u(m-1)(2N);η(m-1)].
The matrix coefficients W~ of the above system is sparse, and for solving system of linear equations W~Y=F~, we can use some iterative Krylov subspace methods and determine the vector Y and hence the approximate solution yN(t)=F(t)Y is obtained.


### 3.2. 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 y(t)≈yN(t)=F(t)Y. According to (16), we have y(t+αk)≈yN(t+αk)=F(t+αk)Y=F(t)(D+αk)TY for k=0,1,…,m. Similar to the previous method, the known function f(t) can be approximated by truncated Fourier series in the form f(t)≈fN(t)=F(t)F where the entries of F would be calculated from (7). By substituting these matrix forms in (4), we get
(29)∑k=0m‍γkF(t)(D+αk)TY=F(t)F.
Hence, the matrix equation corresponding to the basic equation (4) can be written as AY=F or in the form of augmented matrix [A;F], where
(30)A=[aij]={∑k=0m‍γk(D+αk)T}.
Therefore, the basic equation (4) is changed into a system of (2N+1) linear algebraic equations with unknown coefficients y0,y1,…,yN,y^1,…,y^N which can be written in the form
(31)AY=F,
or in augmented matrix form
(32)[A;F]=[a00a01a02…a0(2N);f0a10a11a12…a1(2N);f1a20a21a22…a2(2N);f2⋮⋮⋮⋱⋮;⋮a(2N)0a(2N)1a(2N)2…a(2N)(2N);f^N].
According to the structure of Fourier operational matrices of transmissions, the matrix coefficients A of the above system is tridiagonal (see for instance (14)), and for solving system of linear equations AY=F, we can use some iterative krylov subspace methods and determine the vector Y and hence the approximated solution yN(t)=F(t)Y is obtained.


### 3.3. 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). To do this, let the solution of (5) be approximated by the first (2N+1)-terms truncated Fourier series in the form y(t)≈yN(t)=F(t)Y, where the Fourier coefficient vector Y is unknown (should be determined) and the Fourier vector F(t) is given in (9). According to (12), the kth derivative of yN(t) can be expressed in the matrix form yN(k)(t)=F(k)(t)Y=F(t)(MT)kY.By substituting the above approximations into (5), we get
(33)F(t)(MT)mY=∑j=0J‍∑k=0m-1‍pjk(t)F(λjkt+μjk)(MT)kY+g(t).
To find the solution yN(t), we first collocate (33) at the nodes ti,i=0,1,…,2N-m as the roots of P^2N-m+1(t) (where P^2N-m+1(t) denotes the shifted Legendre polynomial of (2N-m+1) degree in the interval [-π,π]) yields
(34)F(ti)(MT)mY=∑j=0J‍∑k=0m-1‍pjk(ti)F(λjkti+μjk)(MT)kY+g(ti)nnnnnnnnnnnnnnnnnnnnni=0,1,…,2N-m,
where
(35)F(λjkt+μjk)=[12cos(λjkt+μjk)cos(2(λjkt+μjk))⋯sin(N(λjkt+μjk))].
The matrix system (34) can be rewritten as follows:
(36){D(MT)m-∑j=0J‍∑k=0m-1‍PjkE(λjk,μjk)(MT)k}Y=G,
where
(37)D=[F(t0)F(t1)⋮F(t2N-m)]=[12cos(t0)cos(2t0)…sin(Nt0)12cos(t1)cos(2t1)…sin(Nt1)⋮⋮⋮⋱⋮12cos(t2N-m)cos(2t2N-m)…sin(Nt2N-m)],E(λjk,μjk)=[F(λjkt0+μjk)F(λjkt1+μjk)⋮F(λjkt2N-m+μjk)]=[F(λjkt0)D+μjkF(λjkt1)D+μjk⋮F(λjkt2N-m)D+μjk]=[F(λjkt0)F(λjkt1)⋮F(λjkt2N-m)]D+μjk=[12cos(λjkt0)cos(2λjkt0)…sin(Nλjkt0)12cos(λjkt1)cos(2λjkt1)…sin(Nλjkt1)⋮⋮⋮⋱⋮12cos(λjkt2N-m)cos(2λjkt2N-m)…sin(Nλjkt2N-m)]D+μjk,Pjk=[pjk(t0)00…00pjk(t1)0…000pjk(t2)…0⋮⋮⋮⋱⋮000…pjk(t2N-m)],G=[g(t0)g(t1)g(t2)⋮g(t2N-m)].
Hence, the matrix equation (36) corresponding to the basic equation (5) can be written as WY=G or in the form of augmented matrix [W;G], where
(38)W=[wij]={D(MT)m-∑j=0J‍∑k=0m-1‍PjkE(λjk,μjk)(MT)k}.
Therefore, the basic equation (5) is changed into a system of (2N-m+1) linear algebraic equations with unknown coefficients y0,y1,…,yN,y^1,…,y^N, which can be written in the form
(39)WY=G,
or in augmented matrix form
(40)[W;G]=[w00w01w02…w0(2N);g(t0)w10w11w12…w1(2N);g(t1)w20w21w22…w2(2N);g(t2)⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(N-m)(2N);g(t2N-m)].Now we should impose the initial conditionsy(k)(-π)=θk to the associated algebraic system of equations in a proper manner. These initial conditions can be written as
(41)F(-π)(MT)kY=[θk],k=0,1,…,m-1.
Therefore, the matrix form of the initial conditions can be written as
(42)UkY=[θk]or[Uk;θk],k=0,1,…,m-1,
where
(43)Uk=F(-π)(MT)k=[uk0,uk1,…,uk(2N)],k=0,1,…,m-1.
Finally, from (40) and (42), the generalized Pantograph equation (5) subject to the initial conditions y(k)(-π)=θk reduces to the following system of algebraic equations:
(44)W~Y=G~,
where
(45)[W~;G~]=[w00w01w02…w0(2N);g(t0)w10w11w12…w1(2N);g(t1)w20w21w22…w2(2N);g(t2)⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(2N-m)N;g(t2N-m)u00u01u02…u0(2N);θ0u10u11u12…u1(2N);θ1u20u21u22…u2(2N);θ2⋮⋮⋮⋱⋮;⋮u(m-1)0u(m-1)1u(m-1)2…u(m-1)(2N);θm-1].
The matrix coefficients W~ of the above system may be sparse, and for solving system of linear equations W~Y=F~, we can use some iterative krylov subspace methods and determine the vector Y and hence the approximated solution yN(t)=F(t)Y 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.


## 3.1. 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 (2N+1)-terms of the truncated Fourier series. Hence, if we write
(17)y(t)≈yN(t)=y02+∑n=1N‍(yncos(nt)+y^nsin(nt))=F(t)Y,t∈[-π,π],2N+1>m,
where the unknown Fourier coefficient vector Y and the Fourier vector F(t) are given in (9), then the kth derivative of yN(t) can be expressed in the matrix form by yN(k)(t)=F(k)(t)Y=F(t)(MT)kY for k=0,1,…m. Moreover, we need to approximate the given function f(t) which exists in (3) by using (7). In other words, we assume that
(18)f(t)≈fN(t)=f02+∑n=1N‍(fncos(nt)+f^nsin(nt))=F(t)F,t∈[-π,π],
where F=[f0f1f2⋯fNf^1f^2⋯f^N]T is a known vector with the aid of (7). By substituting the approximations (12) and (18) into (3), we get
(19)F(t)(MT)mY=∑k=0m-1‍βkF(t)(MT)kY+F(t)F.
According to the completeness property of Fourier series, the matrix equation corresponding to the basic equation (3) can be written as WY=F or in the form of augmented matrix [W;F], where
(20)W=[wij]={(MT)m-∑k=0m-1‍βk(MT)k}.
Therefore, the basic equation (3) is changed into a system of (2N+1) linear algebraic equations with unknown coefficients y0,y1,…,yN,y^1,…,y^N which can be written in the form
(21)WY=F,
or in augmented matrix form
(22)[W;F]=[w00w01w02…w0(2N);f0w10w11w12…w1(2N);f1w20w21w22…w2(2N);f2⋮⋮⋮⋱⋮;⋮w(2N)0w(2N)1w(2N)2…w(2N)(2N);f^N].
Now we should impose the initial conditions y(k)(-π)=ηk for k=0,1,…,m-1 to the associated system of algebraic equations WY=F in a proper way. For this reason and for clarity of presentation, we can remove the last m rows of W and the last m entries of F without loss of generality. However, we can do this work by the first m rows of the W and the first m entries of F. Thus, we have
(23)[W^;F^]=[w00w01w02…w0(2N);f0w10w11w12…w1(2N);f1w20w21w22…w2(2N);f2⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(2N-m)(2N);f^N-m].
According to (12), we can write y(k)(-π)=F(-π)(MT)kY and hence the matrix-vector representation of the initial conditions is as follows:
(24)F(-π)(MT)kY=[ηk],k=0,1,…,m-1.
Therefore, the matrix form of the initial conditions can be written as
(25)UkY=[ηk]or[Uk;ηk],k=0,1,…,m-1,
where
(26)Uk=F(-π)(MT)k=[uk0,uk1,…,uk(2N)],k=0,1,…,m-1.
Finally, from (23) and (25), (3) subject to the considered initial conditions reduces to the following system of algebraic equations:
(27)W~Y=F~,
where
(28)[W~;F~]=[w00w01w02…w0(2N);f0w10w11w12…w1(2N);f1w20w21w22…w2(2N);f2⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(2N-m)(2N);f^N-mu00u01u02…u0(2N);η0u10u11u12…u1(2N);η1u20u21u22…u2(2N);η2⋮⋮⋮⋱⋮;⋮u(m-1)0u(m-1)1u(m-1)2…u(m-1)(2N);η(m-1)].
The matrix coefficients W~ of the above system is sparse, and for solving system of linear equations W~Y=F~, we can use some iterative Krylov subspace methods and determine the vector Y and hence the approximate solution yN(t)=F(t)Y is obtained.


## 3.2. 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 y(t)≈yN(t)=F(t)Y. According to (16), we have y(t+αk)≈yN(t+αk)=F(t+αk)Y=F(t)(D+αk)TY for k=0,1,…,m. Similar to the previous method, the known function f(t) can be approximated by truncated Fourier series in the form f(t)≈fN(t)=F(t)F where the entries of F would be calculated from (7). By substituting these matrix forms in (4), we get
(29)∑k=0m‍γkF(t)(D+αk)TY=F(t)F.
Hence, the matrix equation corresponding to the basic equation (4) can be written as AY=F or in the form of augmented matrix [A;F], where
(30)A=[aij]={∑k=0m‍γk(D+αk)T}.
Therefore, the basic equation (4) is changed into a system of (2N+1) linear algebraic equations with unknown coefficients y0,y1,…,yN,y^1,…,y^N which can be written in the form
(31)AY=F,
or in augmented matrix form
(32)[A;F]=[a00a01a02…a0(2N);f0a10a11a12…a1(2N);f1a20a21a22…a2(2N);f2⋮⋮⋮⋱⋮;⋮a(2N)0a(2N)1a(2N)2…a(2N)(2N);f^N].
According to the structure of Fourier operational matrices of transmissions, the matrix coefficients A of the above system is tridiagonal (see for instance (14)), and for solving system of linear equations AY=F, we can use some iterative krylov subspace methods and determine the vector Y and hence the approximated solution yN(t)=F(t)Y is obtained.


## 3.3. 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). To do this, let the solution of (5) be approximated by the first (2N+1)-terms truncated Fourier series in the form y(t)≈yN(t)=F(t)Y, where the Fourier coefficient vector Y is unknown (should be determined) and the Fourier vector F(t) is given in (9). According to (12), the kth derivative of yN(t) can be expressed in the matrix form yN(k)(t)=F(k)(t)Y=F(t)(MT)kY.By substituting the above approximations into (5), we get
(33)F(t)(MT)mY=∑j=0J‍∑k=0m-1‍pjk(t)F(λjkt+μjk)(MT)kY+g(t).
To find the solution yN(t), we first collocate (33) at the nodes ti,i=0,1,…,2N-m as the roots of P^2N-m+1(t) (where P^2N-m+1(t) denotes the shifted Legendre polynomial of (2N-m+1) degree in the interval [-π,π]) yields
(34)F(ti)(MT)mY=∑j=0J‍∑k=0m-1‍pjk(ti)F(λjkti+μjk)(MT)kY+g(ti)nnnnnnnnnnnnnnnnnnnnni=0,1,…,2N-m,
where
(35)F(λjkt+μjk)=[12cos(λjkt+μjk)cos(2(λjkt+μjk))⋯sin(N(λjkt+μjk))].
The matrix system (34) can be rewritten as follows:
(36){D(MT)m-∑j=0J‍∑k=0m-1‍PjkE(λjk,μjk)(MT)k}Y=G,
where
(37)D=[F(t0)F(t1)⋮F(t2N-m)]=[12cos(t0)cos(2t0)…sin(Nt0)12cos(t1)cos(2t1)…sin(Nt1)⋮⋮⋮⋱⋮12cos(t2N-m)cos(2t2N-m)…sin(Nt2N-m)],E(λjk,μjk)=[F(λjkt0+μjk)F(λjkt1+μjk)⋮F(λjkt2N-m+μjk)]=[F(λjkt0)D+μjkF(λjkt1)D+μjk⋮F(λjkt2N-m)D+μjk]=[F(λjkt0)F(λjkt1)⋮F(λjkt2N-m)]D+μjk=[12cos(λjkt0)cos(2λjkt0)…sin(Nλjkt0)12cos(λjkt1)cos(2λjkt1)…sin(Nλjkt1)⋮⋮⋮⋱⋮12cos(λjkt2N-m)cos(2λjkt2N-m)…sin(Nλjkt2N-m)]D+μjk,Pjk=[pjk(t0)00…00pjk(t1)0…000pjk(t2)…0⋮⋮⋮⋱⋮000…pjk(t2N-m)],G=[g(t0)g(t1)g(t2)⋮g(t2N-m)].
Hence, the matrix equation (36) corresponding to the basic equation (5) can be written as WY=G or in the form of augmented matrix [W;G], where
(38)W=[wij]={D(MT)m-∑j=0J‍∑k=0m-1‍PjkE(λjk,μjk)(MT)k}.
Therefore, the basic equation (5) is changed into a system of (2N-m+1) linear algebraic equations with unknown coefficients y0,y1,…,yN,y^1,…,y^N, which can be written in the form
(39)WY=G,
or in augmented matrix form
(40)[W;G]=[w00w01w02…w0(2N);g(t0)w10w11w12…w1(2N);g(t1)w20w21w22…w2(2N);g(t2)⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(N-m)(2N);g(t2N-m)].Now we should impose the initial conditionsy(k)(-π)=θk to the associated algebraic system of equations in a proper manner. These initial conditions can be written as
(41)F(-π)(MT)kY=[θk],k=0,1,…,m-1.
Therefore, the matrix form of the initial conditions can be written as
(42)UkY=[θk]or[Uk;θk],k=0,1,…,m-1,
where
(43)Uk=F(-π)(MT)k=[uk0,uk1,…,uk(2N)],k=0,1,…,m-1.
Finally, from (40) and (42), the generalized Pantograph equation (5) subject to the initial conditions y(k)(-π)=θk reduces to the following system of algebraic equations:
(44)W~Y=G~,
where
(45)[W~;G~]=[w00w01w02…w0(2N);g(t0)w10w11w12…w1(2N);g(t1)w20w21w22…w2(2N);g(t2)⋮⋮⋮⋱⋮;⋮w(2N-m)0w(2N-m)1w(2N-m)2…w(2N-m)N;g(t2N-m)u00u01u02…u0(2N);θ0u10u11u12…u1(2N);θ1u20u21u22…u2(2N);θ2⋮⋮⋮⋱⋮;⋮u(m-1)0u(m-1)1u(m-1)2…u(m-1)(2N);θm-1].
The matrix coefficients W~ of the above system may be sparse, and for solving system of linear equations W~Y=F~, we can use some iterative krylov subspace methods and determine the vector Y and hence the approximated solution yN(t)=F(t)Y 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.


## 4. 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) in the form L[y(t)]=g(t), where
(46)L[y(t)]=y(m)(t)-∑j=0J‍∑k=0m-1‍pjk(t)y(k)(λjkt+μjk),t∈[-π,π],
with the initial conditions y(k)(-π)=θk,k=0,1,…,m-1. On the other hand, our collocation scheme finds the approximate solution yN(t) of the problem L[yN(t)]=g(t)+RN(t), where the residual function RN(t) vanishes at the collocation nodes and
(47)L[yN(t)]=yN(m)(t)-∑j=0J‍∑k=0m-1‍pjk(t)yN(k)(λjkt+μjk),t∈[-π,π],
with the initial conditions yN(k)(-π)=θk,k=0,1,…,m-1. Now, we define the error function as eN(t)=y(t)-yN(t). Trivially one can write L[eN(t)]=L[y(t)]-L[yN(t)]=-RN(t) with the initial conditions e(k)(-π)=0,k=0,1,…,m-1. In other words, we construct the error problem as follows:
(48)eN(m)(t)-∑j=0J‍∑k=0m-1‍pjk(t)eN(k)(λjkt+μjk)=-RN(t),t∈[-π,π],
with the initial conditions e(k)(-π)=0,k=0,1,…,m-1. Similar to the previous collocation method, we now solve the above error problem and approximate eN(t) by the aid of truncated Fourier series (with the 2M^+1 terms) in the following form:
(49)eN(t)≈eN,M^(t)=y0*2+∑n=1M^‍(yn*cos(nt)+y^n*sin(nt)),t∈[-π,π],M^>N,2M^+1>m.
If the exact solution of (5) is unknown, then the absolute errors |eN(ti)|=|y(ti)-yN(ti)| are not found. However, the absolute errors can be computed approximately with the aid of the estimated absolute error function |eN,M^(t)|.


## 5. 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 solutiony(t), the approximate solution yN(t), the absolute error function eN(t)=|y(t)-yN(t)|, and estimate absolute error function eN,M^(t) (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–18] could not obtain the exact solutions.Example 2.
As the first example we consider the following ODE with the exact solutiony(t)=sin(t)(50)y′(t)+y(t)=sin(t)+cos(t),t∈(-π,π),
subject to the initial condition y(-π)=0.
According to (8), one can approximate the solution of the above ODE with assumption N=1 as follows:
(51)y(t)≈y1(t)=y02+y1cos(t)+y^1sin(t)=[12cos(t)sin(t)][y0y1y^1]=F(t)Y.
Our aim is to determine the vector components Y with the aid of Fourier operational matrix of differentiation. By using (10) (i.e., F′(t)=F(t)MT), we have
(52)y′(t)=F(t)′Y=F(t)MTY=[12cos(t)sin(t)][0000010-10]Y.
Moreover, we can approximate the known function f(t)=sin(t)+cos(t) by the aid of the truncated Fourier series in the form
(53)sin(t)+cos(t)≈[12cos(t)sin(t)][011].
By substituting the matrix forms of (51)–(53) into (50) we have
(54)([0000010-10]+[100010001])Y=[1000110-11]Y=[011].
Since y(-π)=0, then F(-π)Y=[1/2-10]Y=0. Thus, the above matrix equation transforms into the following matrix equation:
(55)[10001112-10]Y=[010].
The solution of the above system is Y=[001]T and hence
(56)yN(t)=F(t)Y=[12cos(t)sin(t)][001]=sin(t),
which is the exact solution.Example 3.
As the second example we consider the following linear difference equation with the exact solutiony(t)=cos(t)(57)y(t-π)+y(t+π2)=-sin(t)-cos(t),t∈(-π,π).
According to (8), one can approximate the solution of the above difference equation with assumption N=1 as follows:
(58)y(t)≈y1(t)=y02+y1cos(t)+y^1sin(t)=[12cos(t)sin(t)][y0y1y^1]=F(t)Y.
Our aim is to determine the vector components Y with the aid of Fourier operational matrix of transmission. By using (16) (i.e., F(t+αk)=F(t)(D+αk)T,  F(t-αk)=F(t)(D-αk)T), we have
(59)y(t-π)=F(t-π)Y=F(t)(D-π)TY=[12cos(t)sin(t)][1000-1000-1]Y,(60)y(t+π2)=F(t+π2)Y=F(t)(D-π/2)TY=[12cos(t)sin(t)][1000010-10]Y.
Moreover, we can approximate the known function f(t)=-sin(t)-cos(t) by the aid of the truncated Fourier series in the form
(61)-sin(t)-cos(t)≈[12cos(t)sin(t)][0-1-1].
By substituting the matrix forms of (59)–(61) into (57), we have
(62)([1000-1000-1]+[1000010-10])Y=[2000-110-1-1]Y=[0-1-1].
The solution of the above system is Y=[010]T and hence
(63)yN(t)=F(t)Y=[12cos(t)sin(t)][010]=cos(t),
which is the exact solution.Example 4.
As the third example, we consider the following Pantograph equation with the exact solutiony(t)=exp(t)sin(2t)(64)y′(t)=exp(t2)y(t2)NNNN+(sin(2t)+2cos(2t)-sin(t))exp(t),nnnnnnnnnnnnnnnnnnnny(0)=0,t∈[0,1].
Since the computational interval of this problem is [0,1], we need the roots of the shifted Legendre polynomial P^2N-1(t)=P2N-1(2t-1) for any arbitrary values of N. Now, we solve this Pantograph equation by means of the proposed method in Section 3 for different values of N such as N=2 and 3. These approximate solutions are as follows:
(65)y2(t)=-6.979679330+10.44580744cos(t)+3.727316997sin(t)-3.466128106cos(2t)-.848624025sin(2t),y3(t)=-5.726862705+9.652054236cos(t)+1.209226542sin(t)-4.336815209cos(2t)+.652069044sin(2t)+.4116236736cos(3t)-.1710465663sin(3t).In Table1, we provide the absolute values of errors for N=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.Table 1
Comparison of the actual and estimate absolute errors forN=2,3 and M^=3,4,5 of Example 4.
t
i
The actual absolute errors
Estimate absolute errors
|
e
2
(
t
i
)
|
=
|
y
(
t
i
)
-
y
2
(
t
i
)
|
|
e
2,3
(
t
i
)
|
|
e
2,4
(
t
i
)
|
0
4.0000
e
-
09
0
3.4000
e
-
10
0.2
2.1080
e
-
04
2.0730
e
-
04
2.1080
e
-
04
0.4
7.7779
e
-
04
7.8162
e
-
04
7.7778
e
-
04
0.6
1.0039
e
-
03
1.0076
e
-
03
1.0039
e
-
03
0.8
1.5349
e
-
04
1.5161
e
-
04
1.5346
e
-
04
1
3.1644
e
-
04
3.1748
e
-
04
3.1646
e
-
04
|
e
3
(
t
i
)
|
=
|
y
(
t
i
)
-
y
3
(
t
i
)
|
|
e
3,4
(
t
i
)
|
|
e
3,5
(
t
i
)
|
0
4.4000
e
-
09
7.0000
e
-
10
8.5000
e
-
10
0.2
3.5096
e
-
06
3.5003
e
-
06
3.5037
e
-
06
0.4
3.8180
e
-
06
3.8078
e
-
06
3.8252
e
-
06
0.6
3.7010
e
-
06
3.6681
e
-
06
3.6935
e
-
06
0.8
1.8690
e
-
06
1.8527
e
-
06
1.8792
e
-
06
1
1.0250
e
-
06
1.0301
e
-
06
1.0370
e
-
06Example 5.
As the final example we consider the linear Pantograph differential equation of first order(66)y′(t)=-y(0.8t)-y(t),y(0)=1.
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] are not efficient as other Taylor, Hermite, and presented Fourier methods. Therefore, for illustrating the efficiency of the proposed method, we depict the history of the residual of the obtained solution of our scheme together with the associated residuals of Taylor and Hermite methods in Figure 1. From this figure, one can see that the presented method is more applicable and efficient for solving such Pantograph equations.Table 2
Comparison between solutions obtained by the present method and the other methods in Example5.
t
i
Walsh  series method [1]
Laguerre series method (N=20) [2]
Taylor series method (N=8) [3]
Hermit series method (N=8) [4]
Presented method (N=6)
0
1.000000
0.999971
1.000000
1.000000
1
0.2
0.665621
0.664703
0.6664691
0.664691
0.6646909997
0.4
0.432426
0.433555
0.433561
0.433561
0.4335607773
0.6
0.275140
0.276471
0.276482
0.276482
0.2764823286
0.8
0.170320
0.171482
0.171484
0.171484
0.1714841124
1
0.100856
0.102679
0.102744
0.102670
0.1026701264Figure 1
Comparison of the residuals of Taylor [3], Hermite [4], and presented methods of Example 5 for N=8.


## 6. Conclusions and Future Works

The aim of this paper is to introduce two new operational matrices of Fourier series for solving a huge class of high-order linear differential equations. It should be noted that operational matrices of differentiation have more accuracy with respect to the integration ones usually. Also Fourier series approximation is very effective and useful for ODEs with periodic solutions. These two advantages together with the superior results of the proposed scheme (in the numerical illustrations) encourage the authors to extend this method for solving high-order linear and nonlinear Fredholm integro-differential equations [16], systems of Volterra integral equations [30], high-order linear partial differential equations [20, 21], and optimal control problems [31].



---

*Source: 198926-2013-04-10.xml*