Numerical solution of systems of differential equations using operational matrix method with Chebyshev polynomials

ABSTRACT In this study, we introduce an effective and successful numerical algorithm to get numerical solutions of the system of differential equations. The method includes operational matrix method and truncated Chebyshev series which represents an exact solution. The method reduces the given problem to a set of algebraic equations including Chebyshev coefficients. Some numerical examples are given to demonstrate the validity and applicability of the method. In Examples, we give some comparison between present method and other numerical methods. The obtained numerical results reveal that given method very good approximation than other methods. Moreover, the modelling of spreading of a non-fatal disease in a population is numerically solved. All examples run the mathematical programme Maple 13.


Introduction
Differential equation and systems are very useful materials both mathematical modelling and to find out some mathematical equations. For instance, the some Fredholm and Volterra integral equations can be transformed into a nonlinear differential equations with conditions. Recently, many problems in applied science are modelled mathematically by using a systems of ordinary differential equations, for example, pollution modelling and its numerical solutions [1], kinetic modelling of lactic acid [2], the prey and predator problem [3,4], modelling of the epidemiological model for computer viruses [5], modelling of mosquito dispersal [6], modelling a thermal explosion [7], dynamical models of happiness [8], stagnation point flow and Lorentz force [9], non-spherical particles sedimentation [10], boundary layer analysis of micropolar dusty fluid with TiO 2 nanoparticles in a porous medium [11], boundary layer flow of an Eyring-Powell non-Newtonian fluid over a linear stretching sheet [12], heat transfer [13], two-phase Couette flow analysis [14].
Due to these reasons, the solutions of these equations are great importance among scientists. To obtain the exact solutions of systems of differential equations is very crucial or impossible the classical methods especially nonlinear forms. Therefore, many scientists need confidential, quick, easy a numerical method. By the given reasons, many scientists have been motivated that scientists have been studied many numerical methods to solve systems of differential equations such as Runga kutta method [11], DTM and DQM [14], collocation method [12,13,15,16], homotopy perturbation method [17], Adomian decomposition method [18], pseudospectral method [19] and other [20][21][22].
In this study, a numerical algorithm is revealed to get numerical results the following systems of differential equations, for i = 1, 2, . . . , m, m j=1 P ij (x, y j , y (1) j , y (2) j , . . . , y using operational matrix method associate the truncated first kind shifted Chebshev polynomials with the (N + 1) terms as: where T * r (x) denotes the Chebyshev polynomials of the first kind, a j r are unknown Chebyshev coefficients and N is chosen any positive integer, f ij (x) and P ij are analytic functions.

Nomenclature
T n (x) The first kind Chebychev polynomials T * n (x) The first kind shifted Chebychev polynomials a j n The unknown coefficients y j N (x) The approximate solutions N 1 The absolute error
These polynomials have the following properties [26][27][28]: (i) [23] T * n+1 (x) has exactly n + 1 real zeroes on the interval [0, 1]. The ith zero x i is (ii) [23] We known that the relation between the powers x n and the shifted Chebyshev polynomials T * n (x) is where ′ denotes a sum whose first term is halved.
The given function y(x) [ L 2 [0, 1] can be approximated as a sum of shifted Chebyshev polynomials as: Now, we consider the truncated shifted Chebyshev polynomials with the first (N + 1) terms as Equation (3). Matrix representation of the approximate solution Equation (3) and its k derivatives are given where By the expression (5) and for n = 0, 1, … N, we find the matrix relation where Then, aid of (8), we get the following relations and To obtain the matrix X (k) (x) in terms of the matrix X(x), we can use the following relation: where Consequently, substituting the matrix forms (8) and (12) into (7) we have the matrix relation of the approximate solution

Method of solution
In this section, we introduce the numerical solution method for Equation (1) with initial conditions Equation (2). We suppose that it can be given the expansion of f i (x) by the shifted Chebyshev polynomials as: Using matrix representation of approximate solution and its derivatives, Equation (1) can be written as: The residual R i (x) for Equation (17) can be written as We use the Tau method in [17,[29][30][31][32][33], Equation (15) can be converted in m(N − (m ij − 1)) linear or nonlinear equations The initial conditions are given by Hence, we have the m(N + 1) sets of equations with m(N + 1) unknowns by Equations (17) and (18). We write the algorithm procedure in the Maple 13 and solve the m(N + 1) sets of equations with m(N + 1) unknowns, and so the approximate solution y j N (x) can be calculated.

Error estimation
We assume that y(x) is a sufficiently smooth function on [0, 1] and I N (x) is the interpolating polynomial to y at x i , where x i , i = 0, 1, . . . , n are the Chebyshev-Gauss grid points, then we have Therefore, we have [17,29,33] |y Theorem: Suppose that the known functions in Equation (1) are real (N + 1) times continuously differential functions on the [0, 1] and are the shifted Chebyshev polynomials expansion of the exact solution. Let be the approximate solution obtained by proposed method, then there exist real number a such that where A = a 0 a 1 · · · a N and A = a 0 a 1 · · · a N . Proof: Let y N (x) is real-valued polynomials of degree ≤ N and y N (x) is the best approximation of y(x). We can write Using (9) we obtain and we have Then from (19), (20) and (21) we found following error bound:
Example 4.2: Secondly, we take the following the nonlinear stiff problem [21][22] with initial conditions y 1 (0) = y 2 (0) = 1. The exact solutions of this problem are Solving this problem given present method for N = 6, 7, 8, we give the absolute errors for N = 6, 7, 8 in Table 1. In Table 2, we give the comparison of some numerical methods and present method. We plot the this numerical results in Figure 1 for y 1 (x) and Figure 2 for y 2 (x). Table 3 displays the maximum norm errors ||y 1 N − y 1 || and ||y 2 N − y 2 ||.

Example 4.3:
We consider the following linear differential equation with initial conditions y 1 (0) = 1 and y 2 (0) = 0. The solutions of this problem are The numerical results are given by Table 4 for y 1 (x) and y 2 (x). Moreover, Obtained numerical results are displayed in Figures 3 and 4.
Obtaining results are displayed in Figures 5-7. The comparison of numerical results of susceptible, infective and recovered populations is displayed in Figures 5-7 with N = 6, 8, 10, respectively. These results are coherent in [38].

Conclusion
In this article, Chebyshev operational method has been applied to numerically solve the systems of differential equations. Given problem has been transformed into an algebraic equations including unknown coefficients of Chebyshev series. The given algorithm has been written in Maple 13 in order to simply solving given examples. Several examples are given to demonstrate the effectiveness and accurate of the numerical method. The obtained results are compared with exact solutions and also the solutions which are obtained by some other numerical schemes in literature. Results say that present method gives acceptedly accurate results from tested problems. In Example 4.5, we numerically solve the modelling of spreading of a non-fatal disease in a population. From Example 4.5, each figure is shown that susceptible populations are decreasing while infective and recovered populations are increasing during t times.