Matrix Based Operatorial Approach to Differential and Integral Problems

are also included in this general form. Moreover, nonlinear problems where the r.h.s. f (x) is replaced by f (x, y(x), y′(x), ..., y(m−1)(x)) can be solved using Newton’s method in the functional space Cm [a, b] by solving a sequence of linear problems (1)+(2). MATLAB uses different methods to solve initial condition problems (ode family) or boundary value problems (bvp4c or bvp5c) based on Runge-Kutta, Adams-Bashforth-Moulton, BDF algorithms, etc. One of the most effective methods for solving (1)+(2) is to shift the problem to the interval [−1, 1] and then to use the Chebyshev spectral methods, i.e. to approximate the solution y by a finite sum of the Chebyshev series


Introduction
Many problems of the real life lead us to linear differential equations of the form with the general conditions where x (1) ij y (j−1) (x are also included in this general form.Moreover, nonlinear problems where the r.h.s.f (x) is replaced by f (x, y(x), y ′ (x), ..., y (m−1) (x)) can be solved using Newton's method in the functional space C m [a, b] by solving a sequence of linear problems (1)+(2).MATLAB uses different methods to solve initial condition problems (ode family) or boundary value problems (bvp4c or bvp5c) based on Runge-Kutta, Adams-Bashforth-Moulton, BDF algorithms, etc.One of the most effective methods for solving (1)+( 2) is to shift the problem to the interval [−1, 1] and then to use the Chebyshev spectral methods, i.e. to approximate the solution y by a finite sum of the Chebyshev series give the 64 approximated eigenvalues successively in the vector L and they are represented in Fig. 1.We see that DMS and Chebfun calculate accurately only a small number of eigenvalues, while Chebpack gives exactly all 64 eigenvalues.The proposed package Chebpack, which is described in this chapter, is based on the representation (3) of the unknown functions and uses the tau method for linear operators (such as differentiation, integration, product with the independent variable,...) and the pseudospectral method for nonlinear operators -nonlinear part of the equations.
The tau method was invented by Lanczos (1938Lanczos ( , 1956) ) and later developed in an operatorial approach by Ortiz and Samara (Ortiz & Samara, 1981).In the past years considerable work has been done both in the theoretical analysis and numerical applications.
Chebpack is freely accessible at https://sites.google.com/site/dvtrif/and at (Trif, 2011).All the running times in this chapter are the elapsed times for a 1.73GHz laptop and for MATLAB 2010b under Windows XP.

Chebpack, basic module
The package contains, at the basic module -level0, the tools which will be used in the next modules.Let us start with the Chebyshev series expansion of a given function y (Boyd, 2000): THEOREM 1.If y is Lipschitz continuous on [-1, 1], it has a unique representation as an absolutely and uniformly convergent series y(x)= c 0 2 T 0 (x)+c 1 T 1 (x)+c 2 T 2 (x)+...,

39
Matrix Based Operatorial Approach to Differential and Integral Problems www.intechopen.comand the coefficients are given by the formula dx, k = 0, 1, 2, .... Let y N−1 be the truncation of the above Chebyshev series and dom =[ −1, 1] be the working interval.Any interval [a, b] can be shifted and scaled to [−1, 1] by using the shifted Chebyshev polynomials First of all, we need a set of N collocation points x 1 , ..., x N ∈ dom in order to find a good transformation between the above spectral approximation (5) of y and its physical representation y(x 1 ), y(x 2 ), ..., y(x N ).Let be the unique polynomial obtained by interpolating y(x) through the points x 1 , ..., x N , see (Trefethen, 2000) for more details on the coefficients a k versus c k .The next theorem of (Boyd, 2000) estimates the error: THEOREM 2. Let y have at least N derivatives on dom.Then for some ξ on the interval spanned by x and the interpolation points.The point ξ depends on the function y, upon N, upon x and upon the location of the interpolation points.
Consequently, the optimal interpolation points are the roots of the Chebyshev polynomial T N (x),(Chebyshev points of the first kind) For these points x, the polynomials {p N−1 } are generally nearly as good approximations to y as the polynomials {y N−1 } and if y is analytic on dom, then both y − y N−1 and y − p N−1 decrease geometrically as N → ∞.This is the spectral convergence property, i.e. the convergence of y − y N−1 and y − p N−1 towards zero is faster than any power of 1 N as N → ∞.Numerical integration and Lagrangian interpolation are very closely related.The standard formulas for a continuous function where w k are the quadrature weights The Gauss quadrature formulas are based on the optimal Legendre points x k , k = 1, ..., N and these formulas are exact for polynomials f up to degree 2N − 1.The idea of Clenshaw-Curtis quadrature is to use Chebyshev points x instead of the above optimal nodes.By using the Chebyshev points of the first kind (7) one obtains the "classical" Clenshaw-Curtis formula while by using the zeros of the first derivative of a Chebyshev polynomial plus the endpoints ±1, i.e. the Chebyshev extrema in [−1, 1] (the so called Chebyshev points of the second kind) one obtains the "practical" Clenshaw-Curtis formula.Both formulas have all the good properties of the Gaussian quadrature, see (Trefethen, 2008) for more details.Consequently, we may use Chebyshev points of the first kind or of the second kind both for quadrature formulas and for physical representation of a function y on [−1, 1].Any interval [a, b] may be scaled to [−1, 1] and we obtain the corresponding formulas.Moreover, by using the mapping x = cos θ and T k (x)=cos kθ we see that the following two series A Chebyshev series is in fact a Fourier cosine series so that the FFT and iFFT may be used to transform the spectral representation of y into the physical one and conversely, the physical representation into the spectral representation.The quadrature weights w could also be calculated by a fast algorithm given in (Waldvogel, 2006).
The first code of level0, inspired from chebpts.m of chebfun (Trefethen et al., 2011) is The input parameters are N -the dimension of the vectors x and w, dom -the computational domain [a, b] and kind which can be 1 or 2 in order to calculate x as the Chebyshev points of the first or of the second kind.Some short tests show the performances of this code.Let us approximate are inspired by chebpolyval.mand chebpoly.mfrom chebfun (Trefethen et al., 2011).These codes perform the correspondence between the spectral representation c of a function f and its physical representation v = f (x) on Chebyshev points of the first or second kind.It is important to remark that linear operators are better represented in the spectral space, while the nonlinear operators are easily handled in the physical space.
In t2x and x2t, c and v are matrices of the same dimension, each column represents the coefficients or the values for some other function, the number of rows is the above dimension N, while kind specifies the type of Chebyshev points used in v.For example, the code n=16;dom=[0,1];kind=2;x=pd(n,dom,kind); vs=sin(x);vc=cos(x);ve=exp(x);c=x2t ([vs,vc,ve],kind); gives in the columns of c the coefficients of the Chebyshev series (6) of sin(x), cos(x) and exp(x) calculated by using the values of these functions on Chebyshev points of the second kind on [0, 1].We remark here that, taking into account the term c 0 T 0 /2, the coefficient c 0 is doubled.
Another code from level0, inspired from bary.m of chebfun (Trefethen et al., 2011) (Berrut & Trefethen, 2004), the barycentric formula is For Chebyshev points one can give explicit formula for barycentric weights w.For the Chebyshev points of the first kind we have MATLAB -A Ubiquitous Tool for the Practical Engineer www.intechopen.comand for the Chebyshev points of the second kind we have We remark that for a general interval dom =[a, b] and if the sign changes for all x k and w k the weights must be multiplicated by ±2 N−1 (b − a) 1−N .This factor cancels out in the barycentric formula so that it is no need to include it.Let us calculate now the differentiation matrix D such that if f is the column of the Chebyshev coefficients of a function f , then Df is the column of the Chebyshev coefficients of the derivative function from where Consequently, D is a sparse upper triangular matrix with Of course, the differentiation could be iterated, i.e. the coefficients of f (p) are D p f.The corresponding code from level0 is where n is the dimension of the matrix D. For dom =[ a, b] the above matrix D is multiplied by 2/(b − a).
Similarly, the code calculates the sparse integration matrix J such that the coefficients of x f (t)dt are Jf.Here the first coefficient of the result Jf may be changed in order to obtain the coefficients for a specific primitive of f .For example, the coefficients of the primitive which vanishes at a = dom(1) are obtained by using J 0 f.The basic formulas for dom =[−1, 1] are

43
Matrix Based Operatorial Approach to Differential and Integral Problems www.intechopen.com As an important example, let us calculate the coeficients of a specific primitive F(x) of the function f (x).We must then solve the initial-value problem If c are the Chebyshev coefficients of F and f are the coefficients of f , the equation is discretized in spectral space as Dc = f.In order to implement the initial condition, we remark that can be written as Tc = α where This means that we can replace the last row of D by T and the last entry of f by α, thus obtaining a new matrix D and a new vector f .Finally, c = D −1 f are the coefficients of the specific primitive.
The following code from level0 T=cpv(n,xc,dom) (chebyshev polynomial values) implements such conditions.Here x c is an arbitrary vector in dom =[a, b] and cpv calculates the values of the Chebyshev polynomials The code is based on the recurrence formulas of Chebyshev polynomials on [−1, 1] The test code Chebpack\examples\ex_level0\quad_ex3 performs these calculations for the special case y ′ = cos x, y(0)=0, with the solution y = sin(x).The coefficients c are obtained by using the differentiation matrix, cc are the coefficients of the exact solution, ccc are obtained by using the integration matrix J and cccc are obtained by using the integration matrix J 0 .
We also remark that if T=cpv(n,x,dom), These transforms are performed by FFT in the codes x2t and t2x, but for a small dimension N we may use this direct matrix multiplication.

44
MATLAB -A Ubiquitous Tool for the Practical Engineer www.intechopen.com As another example, let us calculate the values at the grid points x of the specific primitive which vanishes at a = dom(1) Starting with the values f = f (x) we have the Chebyshev coefficients T −1 f , then J 0 T −1 f are the Chebyshev coefficients of the specific primitive on [a, b] and finally, Another code from level0 X=mult(n,dom) calculates the sparse multiplication matrix X such that if f is the column vector of the Chebyshev coefficients of a function f (x), then Xf is the column vector of the coefficients of the function xf(x).The code is based on the formulas Then, in general, the coefficients of x p f (x) are given by X p f and the coefficients of a(x) f (x) are given by a(X)f for analytical functions a(x), where a(X) is the matricial version of the function a.Moreover, if f (x) x p has no singularity at the origin, then its coefficients are X −p f.Of course, X is a tri-diagonal matrix, X 2 is a penta-diagonal matrix and so on but, generally, the matrix version funm(full(X))of the scalar function a(x) or X −p = [inv(X)] p are not sparse matrices.For a general interval dom =[a, b], X is replaced by b−a 2 X + b+a 2 I N where I N is the sparse unit matrix speye(N).Another method to calculate a(X) is to pass from the values a(x) at the Chebyshev grid x to the Chebyshev coefficients a using x2t and to approximate Here m must be chosen sufficiently large, but m ≤ N so that the known function a(x) is correctly represented by a 0 , a 1 , ..., a m−1 .
In order to calculate the coefficients of the product we may use the formula

45
Matrix Based Operatorial Approach to Differential and Integral Problems www.intechopen.com The needed coefficients are given by Af where the matrix A ≈ a(X) is given by the code A=fact(a,m) from level0.
The test code Chebpack\examples\ex_level0\fact_ex1.m calculates cosh(X) using fact and funm from Matlab.We remark that X is a sparse matrix, so that funm must be applied to full(X).

Chebpack, linear module
At the first level -level1, the package contains subroutines to solve -initial and boundary value problems for high order linear differential equations -initial value problems for first order linear differential systems -linear integral equations -eigenvalues and eigenfunctions for differential problems.
The main method used is the so called tau-method, see (Mason & Handscomb, 2003) or (Canuto et al., 1988) for more theoretical details and the implementation.It is based on -discretization using the differentiation matrix D -discretization using the integration matrix J -splitting the interval dom =[d 1 , d 2 , ..., d p ].

Discretization using the differentiation matrix D
The corresponding code is Chebpack\level1\ibvp_ode.m [x,solnum]=ibvp_ode(n,dom,kind) where n, dom, kind have the same significance as above, x is the Chebyshev grid and solnum is the numerical solution in the physical space calculated at the grid x.The general form for initial/boundary value problems for high order linear differential equations is (1) and its discrete form is where the unknown y is represented in the spectral space by its Chebyshev coefficients c, while b are the Chebyshev coefficients of the r.h.s.f (x).
We remark that the coefficients P k (X) of the equation can be defined in myDE -directly, for example P k (x)=−x gives P k (X)=−X -using f unm, for example P k (x)=sin x gives P k (X)= f unm( full(X),@sin), i.e. using the Taylor series of sin X -using f act, for example P k (x)=cos x gives P k (X)= f act(x2t(cos(x), kind), m), i.e. using the Chebyshev series of cos X -ifP k (x) is a constant, say P k , then The boundary conditions of the general type (2) are implemented using cpv.For example, for y(x 1c ) − y ′ (x 2c )+2y ′′ (x 3c )=y c we calculate T=cpv (N,[x1c,x2c,x3c],dom).One of the last rows of A is replaced by T(1,:)-T(2,:) * D+2 * T(3,:) * D^2 and the corresponding entry of the vector b is replaced by y c .

Discretization using the integration matrix J
The corresponding code is Chebpack\level1\ibvp_ode_int.m [x,solnum]=ibvp_ode_int(n,dom,kind) We remark that the discretization using the differentiation matrix D does not work well for large N.For example, this type of discretization for the problem y(−1)=−2, y(1)=0 with N = 2048 and an error about 5.46 × 10 −11 leads us to a sparse system Ac = b but with a sparsity factor about 25% that increases the computational time to 6.4 sec, see the example ibvp_ode_ex.m in the folder Chebpack\examples\ex_level1.A better idea is to integrate two times the above equation using the much more sparse integration matrix J.This integration make the coefficients c 0 and c 1 arbitrary and we may fix their values by using the boundary conditions, this time at the first two rows of A and b.
Precisely, the first and the second integration of the equation ( 13 ibvp_ode_int_ex.m in the same folder as above gives the same accuracy for the same N = 2048, but needing only 0.12 sec.The new matrix A has now a sparsity factor of about 0.2439% for the dimension 2048.This higher sparsity qualifies the integration method to be used for splitting the interval as well as for differential systems, where the dimension N of matrices is multiplied by the number of subintervals or by the number of differential equations in the system.We give the formulas for first order equations and a general formula.For the first order we have in myDE Generally, if we denote the differentiation operator on functions P by dP, the identity operator by I and the formal k power of the operator I − Jd by [] (k) , we obtain, after m integrations, Ac = b where For example, for m = 3 we have It is important to remark that the absolute value of the Chebyshev coefficients gives us some information about the necessary dimension N of the discretized problem in order to capture the correct behavior of the solution.For example, let us consider the problem ibvp_ode_int_ex2.m in the folder Chebpack\examples\ex_level1 gives the results in Fig. 2. We see that up to a dimension about N = 400, the algorithm cannot resolve y accurately, due to its highly oscillatory behavior.After that, the Chebyshev series begins to converge rapidly.For N = 1024 the elapsed time is about 0.05 sec.

Discretization using the integration matrix J and splitting the interval
The corresponding code is Chebpack\level1\ibvp_ode_split.m  , med = d i+1 +d i 2 .Using the matrix J for the discretization we obtain on each subinterval i = 1, 2, 3 the discretized form A (i) c (i) = b (i) where the matrix A (i) and the vector b (i) are given by ( 14) as above, while c (i) are the Chebyshev coefficients of the solution y (i) on that subinterval i.Now, using the Kronecker product and the reshape command of Matlab, we build a large (but very sparse) system Ac = b ⎛ ⎝ The boundary conditions are now the given boundary conditions say at d 1 and d 4 supplemented by smoothness conditions at d 2 , d 3 15)

Linear first-order systems
Let us consider a first order linear differential system for x ∈ [a, b] If we denote by c =( c (1) , ..., c (m) ) T the Chebyshev coefficients of y 1 (x), ..., y m (x) and by f (1) , ..., f (m) the coefficients of f 1 , ..., f m then the discretized version of the system is Ac = b where f (1)  . . .
The initial conditions are implemented like in ( 15): the last row of each block in the above matrix is replaced by T or zeros for the corresponding columns and the last entry of each block in the r.h.s. is replaced by y ka such that [O, ..., T, ..., O]c = y ka for each k.The corresponding code from level1 is where y0 is the column vector of the initial values [y 1a , ..., y ma ] T .Of course, we may use the integration matrix J instead of D for discretization, obtaining again a system Ac = b where Jf (1)  . . .
with the implementation of the initial conditions on the first row of each block, see ibvp_sys_ex3_int.m, or we may consider systems of higher order, see ibvp_sys_ex2x2.m from the folder Chebpack\examples\ex_level1.

Linear integral equations
Let us consider a Fredholm integral equation The Fredholm integral operator A(y) becomes after discretization with shifted Chebyshev polynomials Consequently, if c =( c 0 , ..., c N−1 ) T are the coefficients of y, then Kc are the coefficients of A(y), given by the matrix K =(k jk ) j,k=0,...,N−1 .
The matrix K can be calculated starting from the physical values In matrix form, this means where T=cpv(n,x,dom) and then we apply x2t, see also (Driscoll, 2010).
The Fredholm integral equation ( 16) becomes (I N − K)c = f where f are the Chebyshev coefficients of f and we obtain the solution by solving this linear system.The corresponding model code from the folder Chebpack\level1 is Similarly, for a Volterra integral equation we obtain, using ( 11), for the Volterra integral operator where y are the values of y(x) at the grid points x.Consequently, the physical representation of the Volterra integral operator is the matrix V phys = TJ 0 T −1 .* K(x i , x j ), see again (Driscoll, 2010) while its spectral representation is Of course, these codes work well only if the kernel K(x, t) is sufficiently smooth such that it can be spectrally represented by an acceptable number of Chebyshev polynomials.

Eigenvalues/eigenfunctions for Sturm-Liouville problems
Let us consider now the second order spectral problem P 2 (x)y ′′ + P 1 (x)y ′ + P 0 (x)y = λR(x)y with homogeneous boundary value conditions as above.Using tau method, we get the following N -dimensional matrix eigenproblem (P 2 (X)D 2 + P 1 (X)D + P 0 (X))c = λR(X)c 52 MATLAB -A Ubiquitous Tool for the Practical Engineer www.intechopen.com of the form Ac = λBc.The conditions give Tc = 0 and combining these equations we derive the generalized eigenproblem where we retain only the first N rows of the matrices to obtain A and B.Here λ * is chosen by the user as a large and known value.We remark that for λ = λ * we get Tc = 0 as above but now the matrix B behaves well.Consequently, the eigenproblem has instead of two infinite eigenvalues two known λ * eigenvalues that can be eliminated.The same procedure is applied for higher order problems.
The corresponding model code from the folder Chebpack\level1 is The folder Chebpack\examples\ex_level1 contains some other examples in the files eig_ode2_ex * , eig_ode3_ex and eig_ode4_ex.An older package LiScEig is also freely accessible at (Trif, 2005).

Chebpack, nonlinear module
At the second level -level2, we have subroutines to solve -initial and boundary value problems for nonlinear differential equations -nonlinear integral equations.
Here the codes of the first level are used at each step of the Newton method applied in the functional space.Another method could be the successive iteration method.

Successive iteration method
Let us consider, as an example, the nonlinear Emden boundary value problem If the starting spectral approximation of y is, for example, Y 0 = 1 then the discretization of the problem is AY = F.Here A = XD 2 + 2D and F = F(x, Y 0 ) are the coefficients of the r.h.s.modified by using the boundary value conditions.We apply successively If Y n converges, then it converges to a solution of the bvp.The Matlab codes are ibvp_ode_succapprox.m from the folder Chebpack\level2 or ibvp_ode_ex1.mfrom the folder Chebpack\examples\ex_level2.
Of course, the discretization could be performed using the integration matrix J instead of D.
Let us consider, for example, the Lotka-Volterra system

53
Matrix Based Operatorial Approach to Differential and Integral Problems www.intechopen.comWe transform this system to integral form and the discretized form in spectral representation is where F is the column vector of the coefficients of y old 1 (t)y old 2 (t) obtained by using t2x.m and x2t.m.In this discretized form we implement the initial conditions as usually and we must solve this system which has a sparsity factor of about 6% for n = 32 and about 3% for n = 64.For a long interval dom given as dom =[d 0 , d 1 , ..., d m ], we apply the successive approximations method at each subinterval.The initial approximation for the following subinterval is given by the final values of the solution for the current subinterval.To test the numerical solution, we remark that the Lotka-Volterra system has as invariant Λ = By 1 + By 2 − C log y 1 − A log y 2 .The code is used by the command [x,solnum]=ibvp_sys_succapprox(32,[0:10:100],2,[0.5,0.5]); and the result is given in Fig. 4, with the value of the invariant Λ = 0.3039720770839.starting with a suitable y old .For each x = x s in the grid and at each iteration, the integral is evaluated as where {t k , k = 1, ..., n} is the Chebyshev grid on [0, b] and {w k , k = 1, ..., n} are the corresponding weights.Consequently, we obtain Taking into account the nonlinearities, all the calculations are performed into the physical space.Next, y old ←− y new until y new − y old < ε.The corresponding code is ibvp_int_succapprox.m from the folder Chebpack\level2.

Newton method
Let us consider again a nonlinear differential problem of the form where L is a linear differential operator, such as Ly(x)=xy ′′ (x)+2y ′ (x) and f is the nonlinear part.We have also the boundary or initial conditions BC/IC.If we denote by P(y)(x)=Ly(x) − f (x, y(x)), BC/IC the problem is of the form P(y)(x)=0 where P is the operator of the problem and y is the unknown.The Newton method works as follows.Starting with an initial approximation y 0 (x) verifying the initial or boundary conditions, we must solve at each step n the linear problem for y n+1 − y n , with the corresponding homogeneous IC/BC.Next, y n+1 = y n + (y n+1 − y n ) and we perform these iterations until y n+1 − y n < ε.For our problem, the linear step is The corresponding code is ibvp_ode_newton.mfrom the folder Chebpack \level2.I t starts with y 0 (x) verifying or not the IC/BC and solves at each step the above linear problem for y n+1 with the nonhomogeneous IC/BC.
A nonlinear system (of order 2 for example) At the linear step it becomes We remark that in the linear step we use the integration which becomes in the discrete form i.e AZ = B.In the l.h.s the matrix J (A(X) − JacF(X, Y n )) is calculated using the code fact.m.This code uses the physical values of A(x) − JacF(x, Y n ) and converts them into spectral coefficients.In the r.h.s the code also starts with the physical values and converts them into their spectral coefficients.The initial conditions are implemented now in the rows 1, n + 1, .... Of course, this code ibvp_sys_newton.m from the folder Chebpack\level2 can be used in a long-term integration algorithm that starts with the initial values y a1 , y b1 , calculates the solution on [a, b], extracts the final values y b1 , y b2 which become initial values for the same code on a new interval [b, c] and so on.
In the case of a nonlinear integral equation (of Fredholm or Volterra) type we start with an initial approximation y 0 (in physical space) and at each Newton step we obtain the linear equation for the correction z We solve this equation in spectral form as in the previous section and the corrected value of y is y 0 + z.We repeat this step until convergence, i.e. until z < ε.There are many examples in the folder Chebpack\examples\level2.

Chebpack, experimental module
Finally, the package contains an experimental level -level3, in progress, for -partial differential equations of evolution type -fractional differential equations (i.e.differential equations of non-integer order).

56
MATLAB -A Ubiquitous Tool for the Practical Engineer www.intechopen.compde_nonlin_ex4.m in the folder Chebpack\examples\ex_level3 for the Korteweg-de Vries problem This problem of the form u t = Lu + N(u) is numerically solved in spectral space c ′ (t)= Lc(t)+N(c(t)) by using the so called implicit exponential time differencing Euler method Here, the first fixed point problem is solved using successive iterations starting with c (k) , where c (k) are the Chebyshev coefficients of the numerical solution at the time level k.

Fractional differential equations
The fractional derivative D q f (x) with 0 < q < 1, 0 < x ≤ b, in the Riemann-Liouville version, is defined by (Podlubny, 1999) and we have Let us consider a function f : [0, b] → IR, with the spectral representation Using the spectral derivative matrix D, we have and using the linearity of the fractional derivative of order q ∈ (0, 1), we obtain By calculating the physical values of the above integrals in the columns k of a matrix I, each row corresponding to a sample of x from the Chebyshev grid, the formula for the fractional derivative is where T is the matrix given by cpv.m.This means that the Caputo fractional differentiation matrix D in physical form is given by i.e.D phys times the vector of physical values of f gives the vector of physical values of D q * f .For the spectral form, i.e.D spec times the vector T −1 f (x) of the coefficients of f gives the vector of the coefficients of D q * f .If q > 1 then let ex be [q] and q will be replaced by q − [q].In this case, the differentiation matrix will be ex+1) .
In order to avoid the singularity of the fractional derivative D q f if f (0) = 0, we suppose that f(0)=0.The problem of computing the integrals for each x of the grid may be solved iteratively, see (Piessens & Branders, 1976).Indeed, we by direct calculation, using the recurrence formula for Chebyshev polynomials, as well as for the derivatives of the Chebyshev polynomials, for k = 3, 4, ..
This recurrence is as stable as the recurrence which calculates the Chebyshev polynomials.The calculation of the fractional derivative matrix D spec is coded in deriv_frac.m from the folder Chebpack\level3.Next, if needed, Using an idea of Sugiura & Hasegawa (Sugiura & Hasegawa, 2009), let J(s; f ) be J(s; f )= s 0 f ′ (t)(s − t) −q dt 60 MATLAB -A Ubiquitous Tool for the Practical Engineer www.intechopen.comand, approximating f (t) by a combination p n (t) of Chebyshev polynomials on (0, 1) we have the approximations |J(s; f ) − J(s; p n )| ∼ O(nρ −n ), for some ρ > 1.Of course, this method is not suitable for the functions with a singularity in [0, 1] or singularities of lower-order derivatives, like x 0.1 for example.In this case, n must be excessively large.For the initial value problems for fractional differential equations, let us consider the problem D q * y(x)=x 2 + 2 Γ(3 − q) x 2−q − y(x), y(0)=0, q = 0.5, x ∈ [0, 1] .
The physical discretization is where, in order to implement the initial condition, the first row of A is replaced by [1, 0, ..., 0] and b(1) is replaced by 0. The solution is now y(x)=A −1 B. The example is coded in deriv_frac_ex1.mfrom the folder Chebpack\examples\ex_level3.

Conclusion
The new package Chebpack (Trif, 2011) is a large and powerful extension of LiScEig (Trif, 2005).It is based on the Chebyshev tau spectral method and it is applied to linear and nonlinear differential or integral problems, including eigenvalue problems.An experimental module contains applications to linear or nonlinear partial differential problems and to fractional differential problems.Each module is illustrated by many examples.A future version will handle also piecewise functions as well as functions with singularities.The following comparisons with MATLAB codes bvp4c, bvp5c as well as with DMS or Chebfun prove the efficiency of Chebpack.The elapsed time was evaluated after three code executions.The first test problem (test1.m in the folder Chebpack\tutorial)is x 1/2 , u(0)=u(1)=0, x ∈ [0, 1] with the exact solution u ex (x)=x 3/2 − x 5/2 .The elapsed time and the errors are presented in ij ∈ [a, b], ∀i, j, k = 1, ..., m.These multipoint conditions include (for m = 2, for example) -initial value conditions, y(a)=β 1 , y ′ (a)=β 2 , -boundary value conditions α 11 y(a)+α 12 y ′ (a)=β 1 , α 21 y(b)+α 22 y ′ (b)=β 2 , -periodic conditions y(a) − y(b)=0, y ′ (a) − y ′ (b)=0.Eigenvalue problems for linear differential operators m ∑ k=1 P k (x) d k y dx k + (P 0 (x) − λw(x)) y = 0, x ∈ [a, b] m ∑ j=1 α

Fig. 1 .
Fig. 1.Eigenvalues for the problem (4) [x,w] = pd(N,dom,kind) (pd means "physical domain").It calculates the grid x =[x 1 , ..., x N ] T (column vector) and the quadrature weights w =[w 1 , ..., w N ] (row vector) for the quadrature formula b a f ) and c =( c 0 , ..., c N−1 ) T , we have f (x)=Tc, for x ∈ [a, b].The code cpv could be used to transform between the spectral representation f of the function f and the physical representation cos(πx) − π xsin(πx) .The discrete form is εI N + JX − J 2 c = J 2 f where c are the Chebyshev coefficients of the solution y and f are the Chebyshev coefficients of the r.h.s.The new code 47 Matrix Based Operatorial Approach to Differential and Integral Problems www.intechopen.com

Fig. 2 .
Fig. 2. The Chebyshev coefficients and the numerical solution for ex2 Sometimes, the solution of a differential problem has a different behavior in different subintervals.For example, for small ε the solution of the problem (13) has a shock near the origin and we need a very large N in order to capture its correct behavior.In these cases it is better to split the working interval dom[a, b] into disjoint subintervals [d 1 , d 2 ] ∪ [d 2 , d 3 ] ∪ ... ∪ [d p−1 , d p ]=[ a, b] adapted to the behavior of the solution.The great advantage is to use a small N for each subinterval.The partial solutions on each subinterval are connected by some level of smoothness.Precisely, let us consider for example a second order differential problem on [a, b] and let us split the interval as [a, b]=[ d 1 , d 2 ] ∪ [d 2 , d 3 ] ∪ [d 3 , d 4 ].This splitting is given by dom = [d 1 , d 2 , d 3 , d 4 ] on input.If we calculate the basic ingredients xs, Xs, Ds, Js for the standard interval [−1, 1], then for each subinterval [d i , d i+1 ], i = 1, 2, 3 the corresponding ingredients become x = len * xs + med, X = len * Xs + med * I N , D = Ds len , J = len * Js where len = d i+1 −d i 2 Fig. 3.The sparsity structure of the matrix A with boundary conditions implemented for 4 subintervals

51
Matrix Based Operatorial Approach to Differential and Integral Problems www.intechopen.com phys T. The Volterra integral equation becomes in physical representation (I N − K)y = f where f are now the values of f at the grid points x.The corresponding model code from Chebpack is [x,solnum]=volt_eq(n,dom,kind) from the folder Chebpack\level1.The folder Chebpack\examples\ex_level1 contains some examples in the files fred_eq_ex * and volt_eq_ex * .

Fig. 4 .
Fig. 4. The solution of the Lotka-Volterra problem In the case of nonlinear integral equations, for example y(x)= f (x)+ b 0 K(x, t, y(t))dt, we perform successive iterations (if this method works) y new (x)= f (x)+ b 0 K(x, t, y old (t))dt

Table 1 .
Here Chebpack uses the differentiation matrix.The second test problem (test2.m in the folder Chebpack\tutorial) is (13) with ε = 10 −4 .elapsedtimeand the errors are presented in Table2.Here Chebpack uses the integration matrix.