# Numerical Solutions to Fractional Perturbed Volterra Equations

Authors: B. Bandrowski; A. Karczewska; P. Rozmej
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

In the paper, a class of perturbed Volterra equations of convolution type with three kernel functions is considered. The kernel functionsgα=tα−1/Γ(α), t>0, α∈[1,2], correspond to the class of equations interpolating heat and wave equations. The results obtained generalize our previous results from 2010.


---

## Body

## 1. Introduction

We study perturbed Volterra equations of the form(1.1)u(x,t)=u(x,0)+∫0t[gα(t-s)+(gα*k)(t-s)]Δu(x,s)ds+∫0tb(t-s)u(x,s)ds,
where x∈ℝd,t>0,gα(t)=tα-1/Γ(α),Γ is the gamma function, gα*k denotes the convolution, α∈[1,2], b,k∈Lloc1(ℝ+;ℝ), 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 considered Volterra equations can be used. Unfortunately, the resolvent approach proposed by Karczewska and Lizama may be used to (1.1) for some particular kernel functions b,k and gα, α∈(0,1/2) for t>1, only. Hence, in our case, the method proposed in [2] cannot be applied for (1.1).Motivation for the study of the fractional integro-differential equations comes from several problems appearing in physics, biology, and/or engineering. There are many phenomena well modeled by deterministic or stochastic fractional equations, see, for example, [3–7]. Several results on stochastic Volterra equations, fractional as well, particularly on the existence of strong solutions to those equations have been obtained by one of us [8–14].Equation (1.1) is an interesting example of using the so-called fractional calculus in the theory of “classical” equations. Let us emphasize that (1.1) is a generalization of the equations which interpolate the heat and wave equations [15, 16]. Two convolutions appearing in (1.1) with the kernel functions b and k, respectively, represent some perturbation acting on the Volterra equation of convolution type.Fractional calculus is a generalization of ordinary differentiation and itegration to arbitrary order [4, 17–19]. There is an increasing interest in applications of fractional calculus in many fields of mathematics [20], mechanics [5, 21, 22], physics [23, 24], and even in biology [6, 25]. A thorough and comprehensive survey of analytical and numerical methods used in solving many problems with applications of fractional calculus is contained in a recent monograph by Baleanu et al. [17].Spectral methods belong to frequently used tools to obtain approximate solutions to complicated problems like fluid dynamic equations, weather predictions, and many others (see e.g., the monograph of Canuto et al. [26]). Recently, these methods have been used as a tool for calculation of fractional derivatives and integrals [27] and to solve Volterra equations with fractional time [28]. In general, spectral methods consist in representation of the solution to the equation under consideration in a finite subspace whereas the exact solution belongs to space of infinite dimension. The method presented in the present paper belongs to that class.The paper is organized as follows. In Section2, a general idea of Galerkin method to integral equations is presented and approximation by the use of finite dimensional Hilbert space is explained. Section 3 presents a system of linear equations obtained from (1.1) by a discrete formulation enabling for numerical solutions. The detailed form of matrices appearing in that approximation is presented for one dimensional case in Section 3.1 and for two spatial dimension case in Section 3.2. The set of basis functions is represented in Section 3.3 and numerical methods used to solve large-scale sparse linear systems are discussed in Section 3.4, as well. Examples of numerical solutions to (1.1) are exhibited and discussed in detail in Section 4, whereas error estimations for the precision of approximate results are given in Section 5.


## 2. Galerkin Method

Let{ϕi:i=1,2,…,∞} represent a set of orthonormal functions on the interval [0,t], spanning a Hilbert space H.Definition 2.1.
Letf,g∈H. The number
(2.1)〈f(t),g(t)〉:=∫0tf(τ)g(τ)Θ(τ)dτ,
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(2.2)∀i,j〈ϕi(t),ϕj(t)〉=δij,
where δij is the Kronecker delta.We are looking for an approximate solution to (1.1) as an element of the subspace Hnϕ, spanned on nϕ first basic functions {ϕj:j=1,2,…,nϕ}(2.3)unϕ(x,t)=∑j=1nϕ‍cj(x)ϕj(t).
For simplicity of notations, let us consider (1.1) in one spatial dimension only. Inserting (2.3) into (1.1), one obtains
(2.4)unϕ(x,t)=u(x,0)+∫0t[a(t-s)+(a*k)(t-s)]Δunϕ(x,s)ds+∫0tb(t-s)unϕ(x,s)ds+ϵnϕ(x,t),
where function ϵnϕ represents the approximation error function. From (2.3) and (2.4), one gets
(2.5)ϵnϕ(x,t)=∑j=1nϕ‍cj(x)ϕj(t)-∫0t[a(t-s)+(a*k)(t-s)]∑j=1nϕ‍d2dx2cj(x)ϕj(s)ds-∫0tb(t-s)∑j=1nϕ‍cj(x)ϕj(s)ds-u(x,0).Definition 2.2.
TheGalerkin approximation of (1.1) is the function unϕ∈Hnϕ, such that ϵnϕ⊥Hnϕ, that is,
(2.6)∀j=1,2,…,n〈ϵnϕ(x,t),ϕj(t)〉=0.It follows from Definitions2.2 and 2.1 and (2.5) that
(2.7)0=∫0t[∑j=1nϕ‍cj(x)ϕj(τ)]ϕi(τ)Θ(τ)dτ-∫0tu(x,0)ϕi(τ)Θ(τ)dτ-∫0t[∫0τ[a(τ-s)+(a*k)(τ-s)]∑j=1nϕ‍d2dx2cj(x)ϕj(s)ds]ϕi(τ)Θ(τ)dτ-∫0t[∫0τb(τ-s)∑j=1nϕ‍cj(x)ϕj(s)ds]ϕi(τ)Θ(τ)dτfori=1,2,…,nϕ.
Therefore
(2.8)∫0tu(x,0)ϕi(τ)Θ(τ)dτ=∫0t[∑j=1nϕ‍cj(x)ϕj(τ)]ϕi(τ)Θ(τ)dτ-∫0t[∫0τ[a(τ-s)+(a*k)(τ-s)]∑j=1nϕ‍d2dx2cj(x)ϕj(s)ds]ϕi(τ)Θ(τ)dτ-∫0t[∫0τb(τ-s)∑j=1nϕ‍cj(x)ϕj(s)ds]ϕi(τ)Θ(τ)dτ,i=1,2,…,nϕ.
Using (2.2), (2.8) can be written in an abbreviated form
(2.9)gi(x)=ci(x)-∑j=1nϕ‍aijd2dx2cj(x)-∑j=1nϕ‍bijcj(x),
where
(2.10)gi(x)=u(x,0)∫0tϕi(τ)Θ(τ)dτ,(2.11)aij=∫0t[∫0τ[a(τ-s)+(a*k)(τ-s)]ϕj(s)ds]ϕi(τ)Θ(τ)dτ,(2.12)bij=∫0t[∫0τb(τ-s)ϕj(s)ds]ϕi(τ)Θ(τ)dτ.
In general aij≠aji.The solution of the set ofnϕ coupled differential equations (2.9) for coefficients cj(x),j=1,2,…,nϕ provides Galerkin approximation (2.3) to (1.1).


## 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:ℝ→ℝ is given by
(3.1)f′′(x)≈f(x-h)-2f(x)+f(x+h)h2+O(h3).
Then the set of equations  (2.9) takes the following form:
(3.2)gi(xl)=ci(xl)+1h2∑j=1nϕ‍aij[-cj(xl-1)+2cj(xl)-cj(xl+1)]-∑j=1nϕ‍bijcj(xl)=ci(xl)+1h2∑j=1nϕ‍[-aijcj(xl-1)+(2aij-h2bij)cj(xl)-aijcj(xl+1)],
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
(3.3)gi(xl,ym)=ci(xl,ym)+1h2∑j=1nϕ‍aij[-cj(xl-1,ym)-cj(xl,ym-1)+4cj(xl,ym)-cj(xl+1,ym)-cj(xl,ym+1)]-∑j=1nϕ‍bijcj(xl,ym)=ci(xl,ym)+1h2∑j=1nϕ‍[(4aij-h2bij)-aijcj(xl-1,ym)-aijcj(xl,ym-1)+(4aij-h2bij)cj(xl,ym)-aijcj(xl+1,ym)-aijcj(xl,ym+1)].
Both sets of linear equations (3.2) and (3.3) can be written in a matrix form
(3.4)g=𝒜c,
where vectors g, c and matrix 𝒜 have block forms
(3.5)g=(G1G2⋮Gnϕ),c=(C1C2⋮Cnϕ),𝒜=([A11]⋯[A1nϕ][A21]⋯[A2nϕ]⋮⋱⋮[Anϕ1]⋯[Anϕnϕ]).
Detailed structure of blocks occurring in (3.5) is given below.


### 3.1. One-Dimensional Case

BlocksGi and Gi are nϕ-dimensional column vectors. For the sake of space, we present their transpositions
(3.6)GiT=(gi(x1),gi(x2),…,gi(xnh)),CiT=(ci(x1),ci(x2),…,ci(xnh)).Blocks[Aij] have the form
(3.7)[Aij]=(μijηij0⋯0θijηijμijηij⋯000ηijμij⋯00⋮⋮⋮⋱⋮⋮000⋯μijηijθij00⋯ηijμij)nh×nh,
where μij=δij+(2/h2)aij-bij, ηij=-(1/h2)aij, whereas
(3.8)θij={ηijforperiodicboundaryconditions,0forclosedboundaryconditions.Vectorsg and c are nϕnh-dimensional, whereas the matrix 𝒜 has dimension nϕnh×nϕnh. The matrix 𝒜 is a sparse one. The number of nonzero elements of the matrix 𝒜 is at most nϕ2(3nh-2) (with closed boundary conditions) or 3nϕ2nh (with periodic boundary conditions).


### 3.2. Two-Dimensional Case

In two-dimensional case,nϕnh2-dimensional vectors GiT and CiT read as
(3.9)GiT=(gi(x1,y1),gj(x1,y2),…,gi(x1,ynh),gi(x2,y1),gi(x2,y2),…,gi(xnh,ynh)),CiT=(ci(x1,y1),ci(x1,y2),…,ci(x1,ynh),ci(x2,y1),ci(x2,y2),…,ci(xnh,ynh)).
Blocks [Aij] have the form of embedded blocks
(3.10)[Aij]=((αij)(βij)(0)⋯(0)(0)(γij)(βij)(αij)(βij)⋯(0)(0)(0)(0)(βij)(αij)⋯(0)(0)(0)⋮⋮⋮⋱⋮⋮⋮(0)(0)(0)⋯(αij)(βij)(0)(0)(0)(0)⋯(βij)(αij)(βij)(γij)(0)(0)⋯(0)(βij)(αij))nh×nh,
where each term (·) is an embedded block of the size nh×nh. In particular
(3.11)(γij)={(βij)forperiodicboundaryconditions,(0)forclosedboundaryconditions,
block (0) is a matrix nh×nh with all null elements, block (αij) is again a sparse matrix of the form
(3.12)(αij)=(μijηij0⋯0θijηijμijηij⋯000ηijμij⋯00⋮⋮⋮⋱⋮⋮000⋯μijηijθij00⋯ηijμij)nh×nh,
where μij=δij+(4/h2)aij-bij, ηij=-(1/h2)aij,
(3.13)θij={ηijforperiodicboundaryconditions,0forclosedboundaryconditions
and block (βij) is diagonal
(3.14)(βij)=(-1h2aij00⋯0000-1h2aij0⋯00000-1h2aij⋯000⋮⋮⋮⋱⋮⋮⋮000⋯-1h2aij00000⋯0-1h2aij0000⋯00-1h2aij)nh×nh.
The matrix 𝒜 is the sparse matrix of nϕnh2×nϕnh2 elements. However, only at most nϕ2(5nh-4)nh (with closed boundary conditions) or 5nϕ2nh2 (with periodic boundary conditions) elements are nonzero.


### 3.3. 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 Pl which are solutions of the Legendre differential equation
(3.15)ddx[(1-x2)ddxPl(x)]+l(l+1)Pl(x)=0
for l=0,1,2,…. The Legendre polynomials are orthogonal on the interval [-1,1] with the weight function Θ≡1(3.16)∫-11Pl(x)Pm(x)dx=22l+1δlm.
Taking the basis function in the form
(3.17)P~l(x)=2j-1tPl-1(2xt-1),l=1,2,3,…,
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.


### 3.4. Methods for Solving Large Linear Systems

The matrices𝒜, 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𝒜 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.


## 3.1. One-Dimensional Case

BlocksGi and Gi are nϕ-dimensional column vectors. For the sake of space, we present their transpositions
(3.6)GiT=(gi(x1),gi(x2),…,gi(xnh)),CiT=(ci(x1),ci(x2),…,ci(xnh)).Blocks[Aij] have the form
(3.7)[Aij]=(μijηij0⋯0θijηijμijηij⋯000ηijμij⋯00⋮⋮⋮⋱⋮⋮000⋯μijηijθij00⋯ηijμij)nh×nh,
where μij=δij+(2/h2)aij-bij, ηij=-(1/h2)aij, whereas
(3.8)θij={ηijforperiodicboundaryconditions,0forclosedboundaryconditions.Vectorsg and c are nϕnh-dimensional, whereas the matrix 𝒜 has dimension nϕnh×nϕnh. The matrix 𝒜 is a sparse one. The number of nonzero elements of the matrix 𝒜 is at most nϕ2(3nh-2) (with closed boundary conditions) or 3nϕ2nh (with periodic boundary conditions).


## 3.2. Two-Dimensional Case

In two-dimensional case,nϕnh2-dimensional vectors GiT and CiT read as
(3.9)GiT=(gi(x1,y1),gj(x1,y2),…,gi(x1,ynh),gi(x2,y1),gi(x2,y2),…,gi(xnh,ynh)),CiT=(ci(x1,y1),ci(x1,y2),…,ci(x1,ynh),ci(x2,y1),ci(x2,y2),…,ci(xnh,ynh)).
Blocks [Aij] have the form of embedded blocks
(3.10)[Aij]=((αij)(βij)(0)⋯(0)(0)(γij)(βij)(αij)(βij)⋯(0)(0)(0)(0)(βij)(αij)⋯(0)(0)(0)⋮⋮⋮⋱⋮⋮⋮(0)(0)(0)⋯(αij)(βij)(0)(0)(0)(0)⋯(βij)(αij)(βij)(γij)(0)(0)⋯(0)(βij)(αij))nh×nh,
where each term (·) is an embedded block of the size nh×nh. In particular
(3.11)(γij)={(βij)forperiodicboundaryconditions,(0)forclosedboundaryconditions,
block (0) is a matrix nh×nh with all null elements, block (αij) is again a sparse matrix of the form
(3.12)(αij)=(μijηij0⋯0θijηijμijηij⋯000ηijμij⋯00⋮⋮⋮⋱⋮⋮000⋯μijηijθij00⋯ηijμij)nh×nh,
where μij=δij+(4/h2)aij-bij, ηij=-(1/h2)aij,
(3.13)θij={ηijforperiodicboundaryconditions,0forclosedboundaryconditions
and block (βij) is diagonal
(3.14)(βij)=(-1h2aij00⋯0000-1h2aij0⋯00000-1h2aij⋯000⋮⋮⋮⋱⋮⋮⋮000⋯-1h2aij00000⋯0-1h2aij0000⋯00-1h2aij)nh×nh.
The matrix 𝒜 is the sparse matrix of nϕnh2×nϕnh2 elements. However, only at most nϕ2(5nh-4)nh (with closed boundary conditions) or 5nϕ2nh2 (with periodic boundary conditions) elements are nonzero.


## 3.3. 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 Pl which are solutions of the Legendre differential equation
(3.15)ddx[(1-x2)ddxPl(x)]+l(l+1)Pl(x)=0
for l=0,1,2,…. The Legendre polynomials are orthogonal on the interval [-1,1] with the weight function Θ≡1(3.16)∫-11Pl(x)Pm(x)dx=22l+1δlm.
Taking the basis function in the form
(3.17)P~l(x)=2j-1tPl-1(2xt-1),l=1,2,3,…,
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.


## 3.4. Methods for Solving Large Linear Systems

The matrices𝒜, 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𝒜 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.


## 4. Examples of Numerical Solutions

In this section, several examples of approximate numerical solutions to (1.1) are presented. The function
(4.1)u(x,0)=11+exp((|x|-r1)/r2),
where |x|=∑i=1d‍xi2 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 r1=3 and r2=0.3 in (4.1) were chosen for a clear graphical presentation of the results.


### 4.1. One-Dimensional Case

For presentation of approximate numerical results we chose an intervalx∈[-10,10] with nh=201 equidistant grid points. The Hilbert space Hnϕ 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.Figure1 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.Figure 1
Numerical solution to (1.1) for α=1.5, t∈[0,5], k(t)=c·cos(t), b≡0.The effect of the perturbation of the formb=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.Figure 2
Numerical solution to (1.1) for α=1.5, t∈[0,5], k≡0, b(t)=c·e-t.


### 4.2. Two-Dimensional Case

In this subsection, we show some results obtained for two-dimensional case. The calculations have been performed on the grid ofnh×nh points, where nh=101, x,y∈[-5,5], with nϕ=20 basic functions spanning the Hilbert space Hnϕ.Figures3, 4, and 5 illustrate several examples of the numerical solutions to (1.1) for α=1.5 and different perturbation functions k and b.Figure 3
Numerical solution to (1.1) for α=1.5, t∈[0,3.5], k≡0, b≡0. Closed boundary conditions.Figure 4
Numerical solution to (1.1) for α=1.5, t∈[0,3.5], k(t)=cos(t), b≡0. Periodic boundary conditions.Figure 5
Numerical solution to (1.1) for α=1.5, t∈[0,3.5], k(t)=cos(t), b(t)=e-t. Periodic boundary conditions.The casek≡b≡0, corresponds to the unperturbed equation. Its solution, as seen from Figure 3, evolves in the wave manner with a significant influence of diffusion due to 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 one-dimensional case, the presence of perturbation in the form ofk(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 withb≠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.


## 4.1. One-Dimensional Case

For presentation of approximate numerical results we chose an intervalx∈[-10,10] with nh=201 equidistant grid points. The Hilbert space Hnϕ 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.Figure1 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.Figure 1
Numerical solution to (1.1) for α=1.5, t∈[0,5], k(t)=c·cos(t), b≡0.The effect of the perturbation of the formb=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.Figure 2
Numerical solution to (1.1) for α=1.5, t∈[0,5], k≡0, b(t)=c·e-t.


## 4.2. Two-Dimensional Case

In this subsection, we show some results obtained for two-dimensional case. The calculations have been performed on the grid ofnh×nh points, where nh=101, x,y∈[-5,5], with nϕ=20 basic functions spanning the Hilbert space Hnϕ.Figures3, 4, and 5 illustrate several examples of the numerical solutions to (1.1) for α=1.5 and different perturbation functions k and b.Figure 3
Numerical solution to (1.1) for α=1.5, t∈[0,3.5], k≡0, b≡0. Closed boundary conditions.Figure 4
Numerical solution to (1.1) for α=1.5, t∈[0,3.5], k(t)=cos(t), b≡0. Periodic boundary conditions.Figure 5
Numerical solution to (1.1) for α=1.5, t∈[0,3.5], k(t)=cos(t), b(t)=e-t. Periodic boundary conditions.The casek≡b≡0, corresponds to the unperturbed equation. Its solution, as seen from Figure 3, evolves in the wave manner with a significant influence of diffusion due to 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 one-dimensional case, the presence of perturbation in the form ofk(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 withb≠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.


## 5. Precision of Numerical Results

A comparison between the analytic and numerical solutions to (1.1) is possible only for one-dimensional 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
(5.1)Δunϕ,nh(t)=max|unϕ,nhanal(xi,t)-unϕ,nhnum(xi,t)|i=1nh,
where maximum is taken over all grid points xi. For d=1 and α=1 and 2, nϕ>20, nh>100,  t≤6 the error estimate Δunϕ,nh(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]).
Ifunϕ(x,t) (see (2.3)) is the Legendre-Gauss-Lobatto interpolation of u(x,t), u∈Hr(I), I=[0,t], t>0, then
(5.2)∥u-unϕ∥L∞(I)≤Cnϕ3/4-r∥u∥Hr(I),r∈ℕ,
where C is a positive constant (see [26]). Therefore, we can get the following error bounds:
(5.3)∥∫0tgα(t-s)(u(x,s)-unϕ(x,s))ds∥L∞(I)≤Cnϕ3/4-r∥u∥Hr(I),α>0,r≥1.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 (Section3). In this case, the numerical error can be estimated by O(h3), 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 to10-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Δϕunϕ,nh(t) denote the maximum difference between two solutions obtained for the same t and the same grid (defined by nh) but with different numbers of basic functions nϕ and nϕ-2.
(5.4)Δϕunϕ,nh(t):=maxi|unϕ,nh(xi,t)-unϕ-2,nh(xi,t)|.Then letΔhunϕ,nh(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) nh and nh-10(5.5)Δhunϕ,nh(t):=maxi|unϕ,nh(xi,t)-u~nϕ,nh-10(xi,t)|,
where u~nϕ,nh-10(xi,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(5.6)Δϕunϕ,nh(t):=maxi,j|unϕ,nh(xi,yj,t)-unϕ-2,nh(xi,yj,t)|,(5.7)Δhunϕ,nh(t):=maxi,j|unϕ,nh(xi,yj,t)-u~nϕ,nh-10(xi,yj,t)|.
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.Figure 6
Error estimate (5.4) for one-dimensional case, for t=1.8, nh=121, α=1,5, k(t)=cos(t), and b(t)=e-t. Periodic boundary conditions.Figure 7
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.Figure 8
Error estimate (5.6) for two-dimensional case, for t=1.8, nh=101, α=1,5, k(t)=cos(t), and b(t)=e-t. Closed boundary conditions.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 ofnϕ and nh.



---

*Source: 529602-2012-12-04.xml*