Abstract and Applied Analysis

Volume 2012 (2012), Article ID 529602, 19 pages

doi:10.1155/2012/529602

Research Article

Numerical Solutions to Fractional Perturbed Volterra Equations

B. Bandrowski,¹ A. Karczewska,¹ and P. Rozmej²

¹Faculty of Mathematics, Computer Science and Econometrics, University of Zielona Góra, ul. Szafrana 4a, 65-516 Zielona Góra, Poland

²Faculty of Physics and Astronomy, University of Zielona Góra, ul. Szafrana 4a, 65-516 Zielona Góra, Poland

Received 8 June 2012; Revised 27 October 2012; Accepted 27 October 2012

Academic Editor: Dumitru Bǎleanu

Copyright © 2012 B. Bandrowski et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In the paper, a class of perturbed Volterra equations of convolution type with three kernel functions is considered. The kernel functions , , , correspond to the class of equations interpolating heat and wave equations. The results obtained generalize our previous results from 2010.

1. Introduction

We study perturbed Volterra equations of the form

where is the gamma function, denotes the convolution, , , 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 and , for , 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 and , 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 Section 2, 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 represent a set of orthonormal functions on the interval , spanning a Hilbert space .

Definition 2.1. Let . The number

where is a weight function, is called the scalar product of functions on the interval .

Let us recall that two functions are orthonormal when

where is the Kronecker delta.

We are looking for an approximate solution to (1.1) as an element of the subspace , spanned on first basic functions

For simplicity of notations, let us consider (1.1) in one spatial dimension only. Inserting (2.3) into (1.1), one obtains

where function represents the approximation error function. From (2.3) and (2.4), one gets

Definition 2.2. The Galerkin approximation of (1.1) is the function , such that , that is,

It follows from Definitions 2.2 and 2.1 and (2.5) that

Therefore

Using (2.2), (2.8) can be written in an abbreviated form

where

In general .

The solution of the set of coupled differential equations (2.9) for coefficients 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 , where . The grid approximation of a second derivative of a function is given by

Then the set of equations (2.9) takes the following form:

where and .

In two-dimensional case, with the grid , where for , the set of equations (2.9) takes the form

Both sets of linear equations (3.2) and (3.3) can be written in a matrix form

where vectors , and matrix have block forms

Detailed structure of blocks occurring in (3.5) is given below.

3.1. One-Dimensional Case

Blocks and are -dimensional column vectors. For the sake of space, we present their transpositions

Blocks have the form

where , , whereas

Vectors and are -dimensional, whereas the matrix has dimension . The matrix is a sparse one. The number of nonzero elements of the matrix is at most (with closed boundary conditions) or (with periodic boundary conditions).

3.2. Two-Dimensional Case

In two-dimensional case, -dimensional vectors and read as

Blocks have the form of embedded blocks

where each term is an embedded block of the size . In particular

block is a matrix with all null elements, block is again a sparse matrix of the form

where , ,

and block is diagonal

The matrix is the sparse matrix of elements. However, only at most (with closed boundary conditions) or (with periodic boundary conditions) elements are nonzero.

3.3. Basis Functions

The basis functions have to be orthogonal on the interval with respect to a weight function . We use the set of Legendre polynomials which are solutions of the Legendre differential equation

for . The Legendre polynomials are orthogonal on the interval with the weight function

Taking the basis function in the form

ensures that the functions fulfill the orthonormality relations (2.2) on the interval . Therfeore they can be used as a basis in the Galerkin method. In principle, any set of functions orthonormal on the interval 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

where and 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 and 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 interval with equidistant grid points. The Hilbert space was spanned on the basis of functions described in Section 3.3.

In our previous study [28], we have shown that when increases from to the solution of unperturbed evolution governed by (1.1) with changes from pure heat (diffusive) behaviour to pure wave motion. Below we present results for fractional cases (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 , where represents an intensity of the perturbation, whereas the function . 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, , (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 , , , .

The effect of the perturbation of the form , when , is presented in Figure 2. As in the former case 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, , produces the opposite effect, the amplitude of decreases, like in the case of dumping.

Figure 2: Numerical solution to (1.1) for , , , .

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 of points, where , , with basic functions spanning the Hilbert space .

Figures 3, 4, and 5 illustrate several examples of the numerical solutions to (1.1) for and different perturbation functions and .

Figure 3: Numerical solution to (1.1) for , , , . Closed boundary conditions.

Figure 4: Numerical solution to (1.1) for , , , . Periodic boundary conditions.

Figure 5: Numerical solution to (1.1) for , , , . Periodic boundary conditions.

The case , 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 , between the pure diffusion case () and the pure wave case (). 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 of 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 influences mainly the amplitude of the solution. Comparing Figure 5 to Figure 4, one notices that the amplitude increases with . Perturbation with the opposite sign 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 or . Despite the existence of the analytic solutions for this case for an arbitrary , given (for 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 . For and and 2, , , the error estimate was always less than .

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 . In a special case, when , 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 (see (2.3)) is the Legendre-Gauss-Lobatto interpolation of , , , , then

where 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 , where 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 , 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 denote the maximum difference between two solutions obtained for the same and the same grid (defined by ) but with different numbers of basic functions and .

Then let denote the maximum difference between two solutions obtained for the same and the same Hilbert space (the same ) but with different numbers of grid points (in one direction) and

where 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

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 , , , , and . Periodic boundary conditions.

Figure 7: Error estimate (5.5) for one-dimensional case, for , , , , and . Closed boundary conditions.

Figure 8: Error estimate (5.6) for two-dimensional case, for , , , , and . 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 of and .

Acknowledgment

A. Karczewska is thankful for support from “Stochastic Analysis Research Network,” Grant PIA-CONICYT-Scientific Research Ring no. 1112.

References

J. Prüss, Evolutionary Integral Equations and Applications, Birkhäuser, 1993.

A. Karczewska and C. Lizama, “Stochastic Volterra equations under perturbations,” preparation

H. Holden, B. Øksendal, J. Ubøe, and T. Zhang, Stochastic Partial Differential Equations, Birkhäuser Boston Inc., Boston, Mass, USA, 1996.

A. Kilbas, H. Srivastava, and J. Trujillo, Theory and Applications of Fractional Differential Equations, vol. 204 of North-Holland Mathematical Studies, Elsevier, 2006.

Leszczyński, An Introduction to Fractional Mechanics, The Publishing Office of Czestochowa University of Technology, Czestochowa, Poland, 2011.

R. L. Magin, Fractional Calculus in Bioengineering, Begel House Inc, Redding, Calif, USA, 2006.

F. Mainardi, “Fractional calculus: some basic problems in continuum and statistical mechanics,” in Fractals and Fractional Calculus in Continuum Mechanics (Udine, 1996), vol. 378, pp. 291–348, Springer, Vienna, Austria, 1997.

A. Karczewska, Convolution Type Stochastic Volterra Equations, vol. 10 of Lecture Notes in Nonlinear Analysis, Juliusz Schauder Center for Nonlinear Studies, 2007.

A. Karczewska and C. Lizama, “On stochastic fractional Volterra equations in Hilbert space,” Discrete and Continuous Dynamical Systems A, pp. 541–550, 2007.

A. Karczewska and C. Lizama, “Stochastic Volterra equations driven by cylindrical Wiener process,” Journal of Evolution Equations, vol. 7, no. 2, pp. 373–386, 2007.

A. Karczewska, “On difficulties appearing in the study of stochastic Volterra equations,” in Quantum probability and related topics, vol. 27, pp. 214–226, World Scientific, 2011.

A. Karczewska and C. Lizama, “Solutions to stochastic fractional relaxation equation,” Physica Scripta T, vol. 136, Article ID 014030, 2009.

A. Karczewska and C. Lizama, “Strong solutions to stochastic Volterra equations,” Journal of Mathematical Analysis and Applications, vol. 349, no. 2, pp. 301–310, 2009.

A. Karczewska and C. Lizama, “Solutions to stochastic fractional oscillation equations,” Applied Mathematics Letters, vol. 23, no. 11, pp. 1361–1366, 2010.

Y. Fujita, “Integrodifferential equation which interpolates the heat equation and the wave equation,” Osaka Journal of Mathematics, vol. 27, no. 2, pp. 309–321, 1990.

W. R. Schneider and W. Wyss, “Fractional diffusion and wave equations,” Journal of Mathematical Physics, vol. 30, no. 1, pp. 134–144, 1989.

D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo, Fractional Calculus Models and Numerical Methods, Complexity, Nonlinearity and Chaos, World Scientific, Boston, Mass, USA, 2012.

I. Podlubny, Fractional Differential Equations, vol. 198, Academic Press, San Diego, Calif, USA, 1999.

S. G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional Integrals and Derivatives, Gordon and Breach Science, Linghorne Pa, USA, 1993.

D. Baleanu and J. I. Trujillo, “A new method of finding the fractional Euler-Lagrange and Hamilton equations within Caputo fractional derivatives,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 5, pp. 1111–1115, 2010.

R. Gorenflo and F. Mainardi, “Fractional calculus: integral and differential equations of fractional order,” in Fractals and Fractional Calculus in Continuum Mechanics (Udine, 1996), vol. 378 of CISM Courses and Lectures, pp. 223–276, Springer, Vienna, Austria, 1997.

J. Sabatier and O. P. Agrawal, Tenreiro Machado, Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, Springer, Berlin, Germany, 2007.

R. Hilfer, Ed., Applications of Fractional Calculus in Physics, World Scientific, Singapore, 2000.

R. Magin, X. Feng, and D. Baleanu, “Solving the fractional order Bloch equation,” Concepts in Magnetic Resonance Part A, vol. 34, no. 1, pp. 16–23, 2009.

R. L. Magin, “Fractional calculus models of complex dynamics in biological tissues,” Computers & Mathematics with Applications, vol. 59, no. 5, pp. 1586–1593, 2010.

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods, Springer, Berlin, Germany, 2006.

C. Li, F. Zeng, and F. Liu, “Spectral approximations to the fractional integral and derivative,” Fractional Calculus and Applied Analysis, vol. 15, no. 3, pp. 383–406, 2012.

B. Bandrowski, A. Karczewska, and P. Rozmej, “Numerical solutions to integral equations equivalent to differential equations with fractional time,” International Journal of Applied Mathematics and Computer Science, vol. 20, no. 2, pp. 261–269, 2010.

R. Barrett, M. Berry, T. F. Chan, and et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 1994.

H. A. van der Vorst, “Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems,” Society for Industrial and Applied Mathematics, vol. 13, no. 2, pp. 631–644, 1992.

Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” Society for Industrial and Applied Mathematics, vol. 7, no. 3, pp. 856–869, 1986.

D. Kincaid and W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, University of Texas at Austin, 2002.