A numerical method for solving systems of higher order linear functional differential equations

1 Introduction

In physics, chemistry, biology and engineering, a lot of problems are modelled by differential equations, delay differential equations [10-13] and their systems [1-9,14-17].

Also, solving these equation may be analytical difficult and therefore numerical techniques are needed. Until now many analytical and numerical solution techniques such as method of steps, Euler’s method, Runga Kutta method, shooting method, spline method [18-20], variational iteration method [6], Adam’s method [8], Adomian decomposition method [21], homotopy perturbation method [14] etc. are used for differential equations and systems of these. For this purpose, in this study we will focus on the numerical solutions of systems of functional differential equations. Some matrix and collocation methods [20, 22-27] which was applied successfully for ordinary differential equations, partial differential equations, integral equations and difference equations previously.

In this paper we will consider the systems of equations in [0,1] of the form,

for i = 1, 2,..., k under the conditions

Here a matrix-collocation method will be applied to the problem to get the approximate solutions in the truncated series form

where L_n(x) are the Müntz-Legendre polynomials defined by the formula

In the problem yj(0)(x)=yj(x) are the unknown functions, ηi,jn,r(x) and g_i(x) functions defined in the interval 0 ≤ x ≤ 1. On the other hand αi,jn,r,βi,jn,r,φn,i,ψn,i,λi ,are real constants, a_j,n, n = 0,1, 2,..., N are unknown Müntz-Legendre coefficients.

2 Matrix relations required for solution method

At the beginning let us consider the equation (1) and try to construct the matrix form of each term in the equation. The approximate solutions y_j(x) the truncated series of Müntz-Legendre polynomials and its derivatives can be written in the matrix form as,

where

On the other hand L(x) matrix can be represented as,

where

T(x) = [1x . . . x^N]

and the matrix F is defined by

Using the relations (5) and (6), y_j(x) can be written in the form

Then substituting (αi,jn,rx+βi,jn,r)  instead of x in Equation (7) yields for j = 1 , 2 , . . . , k

Also the relation between T(αi,jn,rx+βi,jn,r)  and T(x) matrices can be defined by

where by using the binomial expansion,

Besides the relation T(x) and its derivative T^(k)(x) can be written as

T^(k)(x) = T(x) (B^T)^(k)

where

Thus using these relations we obtain the matrix form of

Furthermore, the matrix form of [y^(k)(α^n,r_i,jX + β^n,r_i,j)] is given below:

where

3 Method of solution

In this Section, the method is constructed by using matrix relations in Section 2, collocation points and matrix operations. Firstly, let us write the matrix equation corresponding to the system (1) as ,

Here μ_n,r(x) and g(x) matrices are defined as

and

Using the collocation points [10–13, 18, 19, 21, 26, 28] defined by the following formula

in Equation (10) yields,

So the fundamental matrix [15, 22, 23] equation can be represented as

where

From the relation (9) and the collocation points given by (11) we have

where

Hence, fundamental matrix equation is

The dimensiones of the matrices μ,T,B¯α,β,B¯,F¯ in Equation (14) are k(N +1) χ k(N +1) and the dimensions of A and G are k(N +1) χ 1. Thus the fundamental matrix equation corresponding to the Equation (1) can be written briefly in the form

where

Equation (15) is a system of equations consisting k(N + 1) unknown Müntz-Legendre coefficients and k(N + 1) linear algebraic equations. On the other hand matrix form that corresponds to the conditions (2) can be written in the form of

where for i = 1, 2, . . . , k

So that the matrix form of the conditions can be represented briefly in the form

where

Finally if we replace any rows of [ W; G] by the rows of [U λ] we have the new augmented matrix [24, 25] equation:

that is

If

, we could say

Therefore the unknown Müntz-Legendre coefficients matrix can be And by solving this linear equation system. Finally substituting the coefficients a_j,0, a_j,1,..., a_j,N, j = 1, 2,..., k in Equation (3)gives us the approximate solutions:

4 Error estimation and improved approximate solutions

In this section an error analysis based on the residual function is given for the method defined in the previous section. When the approximate solutions obtained by our method substituted in Equation (1), the equation will be satisfied approximately, that is;

where R_j,N(x) is the residual function [23–26]. Let us write the error function e_j,N = y_j(x) - y_j,N(x) and to obtain the error problem let us subtract the equations in (21) from the equations (1) side by side. So we could reach the error problem

Applying the same procedure for the conditions gives us the new homogenous conditions of the form,

And solving the error problem (22)–(23) by the same method defined in section 3, an approximation e_j,N,M can be find to the error function e_j,N , j = 1, 2,..., k . When solving this error problem it is better to choose M and replacing the matrix

instead of the matrix G . Therefore the approximate solution y_j,N(x) can be improved as y_j,N,M(x) = y_j,N(x) + e_i,N,M(x). So a new error function can be defined as improved absolute error function by the relation |E_j,N,M (x)| = y_j (x) - y_j,N,M (x)|.

(5) Numerical examples

In this section some examples will be given to explain the method in details and to show the numerical results. All the computations and graphs are performed by a code written in MATLAB R2007b.

Example 1.First of all let us consider the system of equations of order two with two unknown functions in 0 ≤ x ≤ 1

under the conditions y1(0)=1,y1(1)(0)=1,y2(0)=1,y2(2)(0)=−1. The exact solutions of the problem are

y₁(x)=e^x and y₂(x)=e^-x. Here k=2,m=2,R=1,

g₁(x)=e^3x/10-1/10+e^-3x/10+1/10+2e^x-1/5-3e^-x+3/10.and g₂(x)=e^x/5-1/10-e^-x/10+3/10+e^x-3/10-4e^-x+1/5.

Now by taking N = 2 for the problem we search for the approximate solutions of the form

The set of the collocation points for N =2 is calculated as

The fundamental matrix equation corresponding to the problem is given by

and the matrices in this equation is defined by