FAST METHODS FOR THE SOLUTION OF SINGULAR INTEGRO-DIFFERENTIAL AND DIFFERENTIAL EQUATIONS

Uniform methods based on the use of the Galerkin method and different Chebyshev expansion sets are developed for the numerical solution of linear integrodifferential equations of the first order. These methods take a total solution time O(N21n N) using N expansion functions, and a_]so provide error extimates which are cheap to compute. These methods solve both singular and regular integro-differential equations. The methods are also used in solving differential equations.

We assume that (x,y) and Pk(X) are either regular or have singularities provided that the singularities are of known and standard form like for example, weak or logarithmic singularities.For a given Chebyshev expansion set {ho(x)} c H, the solution f(x) defines approximations Using Galerkin technique given in Mikhlin [I], we modify the linear system of equations where L L(N) _a(N) g (N)   (1.3) (N) is the (N + 1 x N + 1) leading minor of the matrix L with elements 1 L..13 .dx T i(x) Lhj(x)//l-x p i, j 0,i, N -i (1.4) and g(N) is the leading (N + l)-vector of the vector g with elements i gi I dx ri(x) gx)/l-x 2 i 0,i ,N -i (1.5) (N) T. is the i th Chebyshev polynomial, as usual.Now to determine the vector a we replace the first equation of the system, (1.3), by the boundary condition equation N (N)f.3 a' where fj hj(a) (1.6) Thus the linear system of equations (1.3) is replaced by L,(N) a_(N) g, (N)   (1.7) where for i >_ i the i-chequation is the (i i)-equation in (1.3) and the vector (N) can be determined by solving (1.7).According to [2], provided the set {h.} Further, fN f" In this paper, we consider two different Chebyshev expansion sets: {h. (x) T. (x)} and {h 0(x) I, h I(x) x, hi(x) (i x2) Tj_2 leading to three different methods (I), (II), (III). (i.i0} (x), j -> 2} (i.ii) Methods based on different techniques have been described before for solving integro-differential equations of the first order; Linz [2], Ei-Gendi [3], Abd-elal [4] in all these papers integro-differential equations of the type (i.i) with Kl(X,y) 0 are reduced to integral equations and a quadrature rule is used to establish numerical procedures.All of these methods are limited to integro-diff- erential equations with no f' under the integral sign, also they do not treat boundary conditions in a very uniform way.Ei-Gendi's method [3] used Chebyshev expansion (1.2, i. I0) in approximating the solution of the equation and produce the solution in time 0(N3).The methods we describe in this paper not only over- come these limitations, but also (the last two methods) produce the solution at a cost of total solution time 0(N21n N) and give reliable error estimates which are cheap to compute.Method (I) is a straightforward method in which a Chebyshev expansion set (i.i0) is used to approximate the solution f(x), and then we solve the linear system of equations (1.7) to get the coefficient vector a(N); hence, we consider it a standard method.Method (II) uses a modified Chebyshev expansion set (I.Ii) to approximate f(x) and so we consider it a modified method.Method (III) uses Chebyshev expansion set (I.i0) to approximate not only f(x) by expansion (1.2), but also f (x) d'(N)l T. (x).
The three methods effectively handle singularities in any or all of (x,y), k 0, i, g(x), the solution f(x), and its derivative f (x), provided that the singularities are of known form and have a known Chebyshev expansion (see [b]).
These requirements limit the applicability of the method to those cases where the singularitie-s which appear are of "standard" form-for example, weak singularities or logarithmic singularities.The methods can also treat some other types of singularities modifying the integro-differential equation; for example, a simple pole can be changed to a logarithmic singularity using integration by parts.We give in section 2 the analysis which leads to the structure of the matrix L (N) for the three methods, while a comparison between the convergence rate attained by the methods is given in section 3. Section 4 shows, by example, that in the three cases rapid convergence is obtained.

THE MATRIX L(N).
We wish to investigate the construction of the matrix L (N) for the three where A (k) and B (k) are (N + I x N + i) matrices with elements A.
Bi j(k) dx T i(x)//l-'x '2 dy (x,y) --dy k {hj(y)} (2.3)   -I -I The integrals appearing in Equations(2.2),(2.3) and (1.5) must be approximated numerically.We do this by relating Aij(k), Bij(k), and gi' i, j O,1,...,N to coefficients n the expansions of Pr(X), Kr(x,y)/l'y7 and g(x) Chebyshev respectively.These later coefficients are evaluated numerically using the fast alorithm given by Delves, Abd-Elal, and Hendry [6] in which Fourier transform technique is used Thls algorithm leads to small quadrature errors whether K (x,y) P (x) and r r g(x) are singular or regular functions; also, it takes 0(N21n N) operations for evaluating the coefficients of the expansion for K (x,y) l-y 2 and 0(N in N) for r evaluating the coefficients of the expansions for P (x) and g(x).Indeed to eval- (r) (r) uate Aj BI] and gl of Equations (2.2), (2.3), and (1.5) numerically, let us assume that the functions P (x) and Kr(x,y)/1-y2 have Chebyshev expansions r 2' (r) T i(x) Tj(y) r 0,i (2.5) where the expansion coefficients 1 pj(r) 2 I dx rj(x) Pr(X)//-r 0,I -i satisfy the inequality Ipj(r) -< C r r 0,i and the expansion coefficients which we can replace by the weaker bounds .-Yr(r) (2.9) (2.10) Also, gi of Equation (1.5) has the bound Igil -< G 1 i _> 0 i when j 0 j when j >-1 and Cr, Dr' G are constants, denotes a sum with first term halved.
Knowing the numerical values for the coefficients p.(k), K.. (k) using the fast algorithm [6], we can easily calculate the elements A.
(k) and B (k) i lJ lJ j 0,1 ,N of the matrices A (k) and B (k) for methods I, II, and III as follows" Method (I).In this method we choose the expansion set (i.i0) and by substituting Ai-'3 Ai j(0) 0 for all i,j except A00 i --for i >_i 2 () (i) 0 for all i,j except for j > 1 i When Pl(X) i, then Aij lj and j of different parity where by different parity, we mean one even and one odd.
(i) Aoo lj [s] is the integer part of s, and means halving the term with 2r + 2( j + i j + i 2 2 1) O.

14)
Notice that we take 0(N21n N) operations to get the matrices A (0), B (0), but from (2.12) and (2.14) it is clear that we take 0(N3) operations to obtain the matrices A (I) and B (I) and in general, tliis makes method I take 0(N3) operations to set up the matrix L *(N" unless explicit forms for A (i) (for example, case Pl(X) i) and Bo (i) are achieved; then, it takes 0(N21n N) operations to get the matrix L *(N)" Method II.In this method we choose the expansion set (i.ii).For i 0,i ,N A..
Method (III)..In this method we use expansion set (i.i0), and hence approximate f(x) by (1.2) and f'(x) by (1.13).Using the relation connecting a i and d i [7]: where means that the term Tj l(X)/ 2(j i) 0 for j 0,i.Now consider the unknown vector _.c [c 0, c 1 CN+ 1] where c 0 ao, cj + i dj, j 0 then the integral equation (i.i) is now reduced to where L (N) L(N) c(N) g(N) (2.217 is the N + i x N + 2) matrix defined by (2.1); hence we need one equation more, and this comes from the boundary condition f(a) , which we write in the form N+I Z c'(N) e. (2.30) Corollary (III.i. (0) < A0(i-j) j i > j -0 -i < AO( j i) j j > i Corollary (111.2).

ASYMPTOTICALLY DIAGONAL MATRICES.
Definition: Given an infiniet matrix L, let the matrix F have elements ILijl Fij (ILiil ILjj I) 1/2 (3.1) The matrix L is said to be "Asymptotically lower diagonal (A.L.D.) of type B(p,r;c)" if constants p,r e 0, c > 0 exist such that Fo. _< c -P (i j)-r i > j (3.2) and is said to by "Asymptotically upper diagonal" (A.U.D.) of the same type if F.. -< c -P (ji) -r j > i (3.3)In an obvious notation we shall then refer to systems of type B as: Type B(PL, PU' rL' ru; CL' Cu) (3.4) For systems (A.D.) of type B, Freeman and Delves [8] provide estimates of the convergence rate.In order to compare the convergence rate attained by the three methods of this paper we need to study the matrices given by each methe@.
We construct the elements of the matrix L *(N) where f0 i, fl a, fj (i a 2) Tj_2(a), 2 _< j < N(Ifjl < i, j 0,i N) Method (I): As is clear from equation (2.12) the matrix A (I) is not (A.D.); indeed, for j > the elements Aij(1) increase with j.So L *(N) is not (A.U.D.) and hence the analysis of Freeman and Delves [8] is not applicable to method I, so as we will see later we do not suggest a value for the truncation error.Also we can not use the iterative method given by Delves [5] to solve the linear system (1.7) and hence any standard method for solving linear equations (Gauss elimination method) can be used.

ERROR ESTIMATES.
The error in any of the three methods contains three distinct components: (a) The truncation error due to cutting off the expansion (1.2) at the th N term.
(b) The discretization error stemming from the quadrature errors A B in the matrix A-B, and g in the vector g.
(c) There are also in principle errors arising from the numerical solution of the linear equations. (N) i(N) are the computed coefficients (a) Truncation error estimates: (S + S 2) From Theorems2 and 3 and for Methods (II) and (III), S 2 dominates S and so S + S 2 S 2 Hence for Method (II) For Method (I), since we are unable to apply Theorem I, we do not suggest a value for S + S 2.
To solve the linear system of equations (1.7) or (2.24), for Method II or II< we use an iterative scheme given in [5] with 0(N2) operations.The error due to this iterative solution is small and so we neglect it, but for Method I and accord- ing to Delves [5], we cannot use this iterative method and hence the error could be relatively large due to error cancellation.We will now see that this error cancellation has no serious effect in a numerical example.

NUMERICAL EXAMPLE.
We give in this section the numerical results for a singular integro-differ- ential equation of the first order.i I f' (x) + f(x) g(x) + I K0(x'Y) f(Y) dy + ] K l(x,y) f'(y) dy -I -i K0(x,y) 2 in(x-y) (x-2y + 2xy 2 2y 3) e-2e(x + 1) ln(x + 1) + 2xe with boundary condition f(a) , a E [-i,i].(x) dx where x.] cos (jw/N), j 0,1,...,N and e N fN-fexact (i) As is clear from the analysis, although Method (I) is a standard method, in fact Methods (II) and (III) are preferable because they provide easy error estimates.
(2) All the three methods work very well when applied to the numerical example, and this suggests that Method (I) is probably a stable method.
(3) The three methods can be applied to ordinary differential equations of the first order which have the form (i.I) but with Hence the same analysis holds with the matrix B 0.
(4) The three methods represent a uniform way of treating boundary conditions, so we recommend methods (II) and (III) be extended to include integro-differential equations of the second order.We are now working on this.
(5) A standard Galerkin calculation has an operations count of 0(N 3) for both setting up and solving the linear equations defining the coefficient vector" however, Methods (II) and (III) of this paper need"

Call for Papers
Thinking about nonlinearity in engineering areas, up to the 70s, was focused on intentionally built nonlinear parts in order to improve the operational characteristics of a device or system.Keying, saturation, hysteretic phenomena, and dead zones were added to existing devices increasing their behavior diversity and precision.In this context, an intrinsic nonlinearity was treated just as a linear approximation, around equilibrium points.Inspired on the rediscovering of the richness of nonlinear and chaotic phenomena, engineers started using analytical tools from "Qualitative Theory of Differential Equations," allowing more precise analysis and synthesis, in order to produce new vital products and services.Bifurcation theory, dynamical systems and chaos started to be part of the mandatory set of tools for design engineers.
This proposed special edition of the Mathematical Problems in Engineering aims to provide a picture of the importance of the bifurcation theory, relating it with nonlinear and chaotic dynamics for natural and engineered systems.Ideas of how this dynamics can be captured through precisely tailored real and numerical experiments and understanding by the combination of specific tools that associate dynamical system theory and geometric tools in a very clever, sophisticated, and at the same time simple and unique analytical environment are the subject of this issue, allowing new methods to design high-precision devices and equipment.
Authors should follow the Mathematical Problems in Engineering manuscript format described at http://www .hindawi.com/journals/mpe/.Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking System at http:// mts.hindawi.com/according to the following timetable:

1980
ATEoIATICS SUBJECT CLASSIFICATION CODES.45J05.i. INTRODUCTION.consider the Galerkin solution for those integro-differential equations of the first order having the form Lfg are elements of a Hilbert space H I and L: H -H is linear.

2 x
It has the exact solution f(x) e The computed errors o N are defined to be ON j Nj0 eN2 (xj)/N i e N

First
Round of Reviews March 1, 2009