Numerical Solutions to Fractional Perturbed Volterra Equations

and Applied Analysis 3 where Θ is a weight function, is called the scalar product of functions f, g on the interval 0, t . Let us recall that two functions are orthonormal when ∀i,j 〈 φi t , φj t 〉 δij , 2.2 where δij is the Kronecker delta. We are looking for an approximate solution to 1.1 as an element of the subspaceHnφ , spanned on nφ first basic functions {φj : j 1, 2, . . . , nφ} unφ x, t nφ ∑ j 1 cj x φj t . 2.3 For simplicity of notations, let us consider 1.1 in one spatial dimension only. Inserting 2.3 into 1.1 , one obtains unφ x, t u x, 0 ∫ t 0 a t − s a ∗ k t − s Δunφ x, s ds ∫ t 0 b t − s unφ x, s ds nφ x, t , 2.4 where function nφ represents the approximation error function. From 2.3 and 2.4 , one gets nφ x, t nφ ∑ j 1 cj x φj t − ∫ t 0 a t − s a ∗ k t − s nφ ∑ j 1 d2 dx2 cj x φj s ds − ∫ t 0 b t − s nφ ∑ j 1 cj x φj s ds − u x, 0 . 2.5 Definition 2.2. The Galerkin approximation of 1.1 is the function unφ ∈ Hnφ , such that nφ ⊥ Hnφ , that is, ∀j 1,2,...,n 〈 nφ x, t , φj t 〉 0. 2.6 4 Abstract and Applied Analysis It follows from Definitions 2.2 and 2.1 and 2.5 that 0 ∫ t 0 ⎡ ⎣ nφ ∑ j 1 cj x φj τ ⎤ ⎦φi τ Θ τ dτ − ∫ t 0 u x, 0 φi τ Θ τ dτ − ∫ t 0 ⎡ ⎣ ∫ τ 0 a τ − s a ∗ k τ − s nφ ∑ j 1 d2 dx2 cj x φj s ds ⎤ ⎦φi τ Θ τ dτ − ∫ t 0 ⎡ ⎣ ∫ τ 0 b τ − s nφ ∑ j 1 cj x φj s ds ⎤ ⎦φi τ Θ τ dτ for i 1, 2, . . . , nφ. 2.7 Therefore ∫ t 0 u x, 0 φi τ Θ τ dτ ∫ t 0 ⎡ ⎣ nφ ∑ j 1 cj x φj τ ⎤ ⎦φi τ Θ τ dτ − ∫ t 0 ⎡ ⎣ ∫ τ 0 a τ − s a ∗ k τ − s nφ ∑ j 1 d2 dx2 cj x φj s ds ⎤ ⎦φi τ Θ τ dτ − ∫ t 0 ⎡ ⎣ ∫ τ 0 b τ − s nφ ∑ j 1 cj x φj s ds ⎤ ⎦φi τ Θ τ dτ, i 1, 2, . . . , nφ. 2.8 Using 2.2 , 2.8 can be written in an abbreviated form gi x ci x − nφ ∑ j 1 aij d2 dx2 cj x − nφ ∑ j 1 bijcj x , 2.9 where gi x u x, 0 ∫ t 0 φi τ Θ τ dτ, 2.10 aij ∫ t 0 [∫ τ 0 a τ − s a ∗ k τ − s φj s ds ] φi τ Θ τ dτ, 2.11 bij ∫ t 0 [∫ τ 0 b τ − s φj s ds ] φi τ Θ τ dτ. 2.12 In general aij / aji. The solution of the set of nφ coupled differential equations 2.9 for coefficients cj x , j 1, 2, . . . , nφ provides Galerkin approximation 2.3 to 1.1 . Abstract and Applied Analysis 5 3. Discretization Equations can be solved using discretization in a space variable. In one-dimesional case, let us introduce a grid of points x1, x2, . . . , xnh , where xl − xl−1 h. The grid approximation of a second derivative of a function f : R → R is given by f ′′ x ≈ f x − h − 2f x f x h h2 O ( h3 ) . 3.1and Applied Analysis 5 3. Discretization Equations can be solved using discretization in a space variable. In one-dimesional case, let us introduce a grid of points x1, x2, . . . , xnh , where xl − xl−1 h. The grid approximation of a second derivative of a function f : R → R is given by f ′′ x ≈ f x − h − 2f x f x h h2 O ( h3 ) . 3.1 Then the set of equations 2.9 takes the following form: gi xl ci xl 1 h2 nφ ∑ j 1 aij [−cj xl−1 2cj xl − cj xl 1 ] − nφ ∑ j 1 bijcj xl ci xl 1 h2 nφ ∑ j 1 [ −aijcj xl−1 ( 2aij − hbij ) cj xl − aijcj xl 1 ] , 3.2 where i 1, 2, . . . , nφ and l 1, 2, . . . , nh. In two-dimensional case, with the grid x1, x2, . . . , xnh × y1, y2, . . . , ynh , where xl − xl−1 ym − ym−1 h for l,m 2, 3, . . . , nh, the set of equations 2.9 takes the form gi ( xl, ym ) ci ( xl, ym ) 1 h2 nφ ∑ j 1 aij [−cj(xl−1, ym) − cj(xl, ym−1) 4cj ( xl, ym ) − cj(xl 1, ym) − cj(xl, ym 1)] − nφ ∑ j 1 bijcj ( xl, ym ) ci ( xl, ym ) 1 h2 nφ ∑ j 1 [ − aijcj ( xl−1, ym ) − aijcj(xl, ym−1) ( 4aij − hbij ) cj ( xl, ym ) − aijcj(xl 1, ym) − aijcj(xl, ym 1) ] . 3.3 Both sets of linear equations 3.2 and 3.3 can be written in a matrix form


Introduction
We study perturbed Volterra equations of the form u x, t u x, 0 where x ∈ R d , t > 0, g α t t α−1 /Γ α , Γ is the gamma function, g α * k denotes the convolution, α ∈ 1, 2 , b, k ∈ L 1 loc R ; R , and Δ is the Laplace operator.The perturbation approach to Volterra equations of convolution type has been used by many authors, see, for example, 1 .Such approach may be applied to more general, not necessary convolution equations, too.Recently, perturbed Volterra equations, deterministic and stochastic as well, have been studied for instance by Karczewska and Lizama 2 .The authors consider the class of equations with three kernel functions which satisfy some scalar auxiliary equations.Such condition enables to construct the family of resolvent operators admitted by the Volterra equations.In consequence, the resolvent approach to the where Θ is a weight function, is called the scalar product of functions f, g on the interval 0, t .
Let us recall that two functions are orthonormal when where δ ij is the Kronecker delta.
We are looking for an approximate solution to 1.1 as an element of the subspace H n φ , spanned on n φ first basic functions {φ j : j 1, 2, . . ., n φ } u n φ x, t n φ j 1 c j x φ j t .

2.3
For simplicity of notations, let us consider 1.

2.8
Using 2.2 , 2.8 can be written in an abbreviated form where In general a ij / a ji .The solution of the set of n φ coupled differential equations 2.9 for coefficients c j x , j 1, 2, . . ., n φ provides Galerkin approximation 2.3 to 1.1 .

Discretization
Equations can be solved using discretization in a space variable.In one-dimesional case, let us introduce a grid of points x 1 , x 2 , . . ., x n h , where x l − x l−1 h.The grid approximation of a second derivative of a function f : R → R is given by Then the set of equations 2.9 takes the following form:

3.3
Both sets of linear equations 3.2 and 3.3 can be written in a matrix form g Ac, 3.4 where vectors g, c and matrix A have block forms Detailed structure of blocks occurring in 3.5 is given below.

One-Dimensional Case
Blocks G i and G i are n φ -dimensional column vectors.For the sake of space, we present their transpositions

3.6
Blocks A ij have the form where

Two-Dimensional Case
In two-dimensional case, n φ n 2 h -dimensional vectors G T i and C T i read as

3.9
Blocks A ij have the form of embedded blocks where each term • is an embedded block of the size n h × n h .In particular for periodic boundary conditions, 0 for closed boundary conditions, 3.11 block 0 is a matrix n h × n h with all null elements, block α ij is again a sparse matrix of the form where θ ij η ij for periodic boundary conditions, 0 for closed boundary conditions 3.13 and block β ij is diagonal

3.14
The matrix A is the sparse matrix of n φ n 2 h × n φ n 2 h elements.However, only at most n 2 φ 5n h − 4 n h with closed boundary conditions or 5n 2 φ n 2 h with periodic boundary conditions elements are nonzero.

Basis Functions
The basis functions {φ j : j 1, 2, 3, . ..} have to be orthogonal on the interval 0, t with respect to a weight function Θ.We use the set of Legendre polynomials P l which are solutions of the Legendre differential equation

3.16
Taking the basis function in the form ensures that the functions P l x fulfill the orthonormality relations 2.2 on the interval 0, t .Therfeore they can be used as a basis in the Galerkin method.In principle, any set of functions orthonormal on the interval 0, t can be used.For our purposes, however, the Lagrange polynomials appeared more efficient in practical applications than for instance Chebyshev polynomials.
Abstract and Applied Analysis 9

Methods for Solving Large Linear Systems
The matrices A, both in one-and two-dimensional case, are sparse matrices.In order to obtain a reasonable approximation of the solution to 1.1 , their sizes have to be large.Those facts suggest an application of iterative methods for solving the linear systems 3.4 .
In general, the matrix A is nonsymmetric.We have tested on our examples two iterative methods developed for solving large-scale linear systems of non-symmetric matrices.One of those metods is so called BiCGSTAB BiConjugate Gradient Stabilized method 29, 30 .The other one is the GMRES General Minimal Residual method 29, 31 .In both methods, a suitable preconditioning is necessary.
For cases discused in the paper, the GMRES method appeared to be more efficient.Usually, after a proper choice of auxiliary parameters of calculations, the GMRES was requiring less number of iterations and converging faster than the BiCGSTAB method.

Examples of Numerical Solutions
In this section, several examples of approximate numerical solutions to 1.1 are presented.The function where x 2 i and d is the space dimension has been taken as the initial condition.Such function, which is substantially different than zero only in a finite region, may represent a distribution of the temperature in a rod or plane heated locally or a distribution of gas or liquid particles which may diffuse in a nonhomogeneous medium.The values of constants r 1 3 and r 2 0.3 in 4.1 were chosen for a clear graphical presentation of the results.

One-Dimensional Case
For presentation of approximate numerical results we chose an interval x ∈ −10, 10 with n h 201 equidistant grid points.The Hilbert space H n φ was spanned on the basis of n φ 20 functions described in Section 3.3.
In our previous study 28 , we have shown that when α increases from α 1 to α 2 the solution of unperturbed evolution governed by 1.1 with k b 0 changes from pure heat diffusive behaviour to pure wave motion.Below we present results for fractional cases α 1.5 an intermediate case pointing out the effects of perturbations.For the sake of space, we show a few examples only, explaining the general influences of perturbations on the character of solutions.
Figure 1 illustrates the effect of perturbation in the form k t c•cos t , where c ∈ 0, 1 represents an intensity of the perturbation, whereas the function b ≡ 0. It is clearly seen that this perturbation, periodic in time variable, produces a wavy formations in space with amlitudes strongly depending on the intensity of perturbation.Results obtained with the opposite sign of the perturbation term, c → −c, not presented here for the sake of space , show that such perturbation decreases diffusive behaviour of the system and enforces a wavelike evolution.The effect of the perturbation of the form b c • e −t , when k ≡ 0, is presented in Figure 2. As in the former case c ∈ 0, 1 represents the intensity of perturbation.It is clear that the perturbation in such form generally increases the amplitude of the solution.The change of sign of the perturbation, c → −c, produces the opposite effect, the amplitude of u x, t decreases, like in the case of dumping.

Two-Dimensional Case
In this subsection, we show some results obtained for two-dimensional case.The calculations have been performed on the grid of n h × n h points, where n h 101, x, y ∈ −5, 5 , with n φ 20 basic functions spanning the Hilbert space H n φ .
Figures 3, 4 fractional value of α 1.5, between the pure diffusion case α 1 and the pure wave case α 2 .For more examples of unperturbed evolution of solutions to 1.1 in two-dimensional case, see 28 .
Results obtained for two-dimensional cases with nonzero perturbations generally exhibit the properties similar to those in one-dimensional case.Again, like in the onedimensional case, the presence of perturbation in the form of k t cos t results in an increase of a wave frequency see, e.g., Figure 4 .In other words the perturbation of such form produces additional wavy behaviour of the solution.The change of sign of the perturbation term changes the phase of that behaviour.
The presence of perturbation term with b / 0 influences mainly the amplitude of the solution.Comparing Figure 5 to Figure 4, one notices that the amplitude increases with b e −t .Perturbation with the opposite sign b t −e −t results in decrease of amplitude, like in the case of dumping.

Precision of Numerical Results
A comparison between the analytic and numerical solutions to 1.1 is possible only for onedimensional case when there is no perturbation and α 1 or α 2. Despite the existence of the analytic solutions for this case for an arbitrary α ∈ 1, 2 , given for d 1 case in terms of Mittag-Leffler functions 15, 16 , their computation is not practical.For non-perturbed case, we defined in 28 an error estimate as the maximum of the absolute value of the difference between the exact analytical solution and approximate numerical one where maximum is taken over all grid points x i .For d 1 and α 1 and 2, n φ > 20, n h > 100, t ≤ 6 the error estimate Δu n φ ,n h t was always less than 10 −5 .
When we consider presented method for obtaining numerical solution to fractional perturbed Volterra equation 1.1 , there are three levels of numerical errors.
The first level corresponds to the error of Galerkin approximation 2.3 which depends on basic functions number n φ .In a special case, when k ≡ 0, b ≡ 0 and Δ operator is replaced by identity operator, we can estimate the approximation error using the following result from 27 .
Remark 5.1 see 27, Remark 5.2 .If u n φ x, t see 2.3 is the Legendre-Gauss-Lobatto interpolation of u x, t , u ∈ H r I , I 0, t , t > 0, then where C is a positive constant see 26 .Therefore, we can get the following error bounds: Also at the first level, the integrals 2.11 and 2.12 are calculated.In our method, we use a Gauss-Legendre quadrature for numerical integration.The exact error of such quadrature can be found, for example, in Theorem 7.3.5 in 32 .
The second level is the Laplacian discretization Section 3 .In this case, the numerical error can be estimated by O h 3 , where h is the spatial grid step.
The last one is is the residual error of GMRES method for solving large linear systems.In our computations, the residual error threshold was set to 10 −9 , which was small enough to obtain reliable solution.
The joint error estimate from all three levels is not obvious.Moreover, in the considered perturbed case, the analytic solution to 1.1 is not known.Therefore, in order to estimate the accuracy of numerical solutions, we proceed in the following manner which is also applicable for two-and higher-dimensional cases.
When we are not able to confront the numerical solutions with analytic ones, we can investigate how does approximate solution change with increasing numbers of grid points and increasing number of basic functions.One can expect that increasing number of grid points and increasing number of basic functions should result in a better approximation of the true unknown solution.Taking appropriate sequences of those numbers and estimating the largest differences between consecutive solutions, one can show convergence of approximation errors.To do it let us define the following quantities.
Let Δ φ u n φ ,n h t denote the maximum difference between two solutions obtained for the same t and the same grid defined by n h but with different numbers of basic functions n φ and n φ−2 .

Δ φ u n φ ,n h t : max
5.4 Then let Δ h u n φ ,n h t denote the maximum difference between two solutions obtained for the same t and the same Hilbert space the same n φ but with different numbers of grid points in one direction n h and n h−10 Δ h u n φ ,n h t : max i u n φ ,n h x i , t − u n φ ,n h −10 x i , t , 5.5 The general conclusions of that investigation are the following.In all cases the error estimates decrease fast with increasing number of grid points or with increasing size of the basis.That decrease is in log plots seen as close to a straight line, that is, error estimates decrease almost exponentially.Then taking large enough grid and large enough set of basic functions, one can obtain, in principle, an error less then arbitrary small number.In practice increasing the sizes of basis and grid produces a sharp increase of numerical operations causing accumultion of rounding errors.That property can be compensated by increasing the precision of representation of real numbers using double or quadruple precision and so on.All those actions require higher and higher computer power to obtain results in a reasonable computing time.
Our estimates show, however, that a reasonable approximations can be obtained with relatively low values of n φ and n h .
− s g α * k t − s Δu x, s ds t 0 b t − s u x, s ds, 1.1
Definition 2.2.The Galerkin approximation of 1.1 is the function u n φ ∈ H n φ , such that n φ ⊥ H n φ , that is, 1 in one spatial dimension only.Inserting 2.3 into 1.1 , one obtains u n φ x, t u x, 0 t 0 a t − s a * k t − s Δu n φ x, s ds Error estimate 5.4 for one-dimensional case, for t 1.8, n h 121, α 1, 5, k t cos t , and b t e −t .Periodic boundary conditions.where u n φ ,n h −10 x i , t means the value in the node of bigger size obtained from values calculated for grid of smaller size by cubic-spline interpolation.In two-dimensional cases, the appropriate error estimates read Δ φ u n φ ,n h t : max n φ ,n h x i , y j , t − u n φ ,n h −10 x i , y j , t .5.7 Figures 6, 7, and 8 present some examples of the dependence of the above defined error estimates on the numbers of grid points and the size of the basis.Presented examples contain both one-and two-dimensional cases and two cases of boundary conditions. u Error estimate 5.5 for one-dimensional case, for t 1.8, n φ 20, α 1, 5, k t cos t , and b t e −t .Closed boundary conditions.Error estimate 5.6 for two-dimensional case, for t 1.8, n h 101, α 1, 5, k t cos t , and b t e −t .Closed boundary conditions.