Abstract
We make use of an adaptive numerical method to compute blow-up solutions for nonlinear ordinary Volterra integrodifferential equations (VIDEs). The method is based on the implicit midpoint method and the implicit Euler method and is named the implicit midpoint-implicit Euler (IMIE) method and was used to compute blow-up solutions in semilinear ODEs and parabolic PDEs in our earlier work. We demonstrate that the method produces superior results to the adaptive PECE-implicit Euler (PECE-IE) method and the MATLAB solver of comparable order just as it did in our previous contribution. We use quadrature rules to approximate the integral in the VIDE and demonstrate that the choice of quadrature rule has a significant effect on the blow-up time computed. In cases where the problem contains a convolution kernel with a singularity we use convolution quadrature.
1. Introduction
Many solutions of differential equations modeling physical problems blow-up in finite time. the phenomenon of blow-up is said to have occurred at if the solution of the differential equation becomes unbounded at . The blow-up time often represents an important change in the properties of the model and it is important that it should be accurately reproduced by a numerical computation.
Much work in the study of blow-up has been on reaction-diffusion equations. The reaction-diffusion equation is a semilinear parabolic equation of the form where is referred to as the diffusion term and as the reaction term. It is well known from blow-up theories that for sufficiently large initial function the solution of (1.1) will blow-up in finite time.
Most methods for solving evolution equations lose accuracy as the solution becomes large (approaches blow-up) and hence adaptive methods are used. Bandle and Brunner [1] present a survey of the theory and the numerical analysis of blow-up solutions for quasilinear reaction-diffusion equations of the form (1.1). Budd et al. [2] proposed the use of moving mesh partial differential equation methods (MMPDE) for solving (1.1). A more recent study of the MMPDE is presented in [3] and the references therein.
Berger and Kohn [4] used a mesh refinement strategy which is able compute the solution accurately over the entire physical interval even as the solution grows in magnitude. Stuart and Floater [5] showed that fixed step methods, both explicit and implicit fail to reproduce blow-up time for a scalar ODE. They also examined variable step methods. They used time stepping strategies which are based on a rescaling of the time variable in the underlying differential equation. They also apply these ideas to a PDE.
In some applications, for example, in population dynamics, the reaction term is nonlocal [6]. For example, so that (1.1) becomes
In this paper, we consider the Volterra integrodifferential equation (VIDE), with and is the kernel.
Recently Ma [7] has studied the finite-time blow-up theory and blow-up conditions for nonlinear ordinary VIDEs of the form (1.3) and nonlinear partial VIDEs. Assuming(a) is an integrable positive function such that , where ;(b) is nonnegative, nondecreasing and continuous for , for , and ,
he proves that the solution of the VIDE (1.3) blows up if and only if for any , where .
The earliest author to look at blow-up solutions of parabolic equations with nonlinear memory is Bellout [8]. The theory of blow-up for nonlocal reaction-diffusion problems is given by authors such as [9, 10]. Ma et al. [11] implement the moving mesh methods to solve reaction-diffusion equations with nonlocal nonlinear terms such as (1.2) which blows-up in finite time.
2. Description of Methods Used
We use quadrature rules to approximate the integral in (1.3) and solve the resulting system of ODEs of the form:
We use the adaptive PECE-implicit Euler (PECE-IE) and adaptive implicit midpoint-implicit Euler (IMIE) methods introduced in Dlamini and Khumalo [12] to solve the resulting system of ODEs and compare the results with matlab ODE solvers ode45 and ode15s. As the blow-up time is approached changes occur on an increasingly smaller time scales and so to be more accurate in obtaining the blow-up time, variable stepsize methods are used. We select a mesh with .
In our procedures we specify the acceptable error per step and if it is not met the procedure adjusts the step size so that each step introduces an error that is not more than the acceptable error. At each step we solve the problem using two different algorithms, giving two different solutions, say and . We approximate the local truncation error by computing the difference between the two solutions, that is, . If the error is smaller than the specified acceptable error, we move to the next step and accept as the solution for the current step. Otherwise, we set a new stepsize and redo the current step.
2.1. Adaptive PECE-Implicit Euler Method (PECE-IE)
This method is based on the implicit Euler method. We compute using a predictor-corrector method in which is obtained using explicit Euler’s method so that we have
To compute we use Newton’s method as the solver to deal with the implicit nature of the implicit Euler method.
2.2. Adaptive Implicit Midpoint-Implicit Euler Method (IMIE)
We use the implicit Euler method with Newton’s method as the solver to get . To compute we use implicit midpoint method given by
First we compute using the implicit Euler method and substitute in (2.4) to get that is,
then,
2.3. Matlab Solvers (ode45 and ode15s)
ode45 is a one-step Matlab solver that is based on an explicit Runge-Kutta (4, 5) scheme. It varies the size of the step of the independent variable in order to meet the accuracy specified. On the other hand, ode15s is a multistep Matlab variable order solver based on implicit methods. We use these solvers to compare the PECE-IE and IMIE methods.
3. Numerical Computation
We consider with and is the kernel.
We solve (3.1) by first approximating the integral using a quadrature rule and then solve the resulting system of ODEs using the methods mentioned above. We use the repeated trapezoidal rule and block-by-block-method (see [13, 14] and the references therein). The effect of the quadrature rule in the computation of the blow-up time will be investigated. In cases where the convolution kernel has a singularity we use convolution quadrature discussed in Lubich [15–17].
3.1. Using Repeated Trapezoidal Rule
We let
The approximation of every integral in (3.3) by repeated trapezoidal rule will yield the following system (with and )
3.2. Using Block-by-Block Method
In this method we use Simpson’s rule and a quadratic interpolation to approximate the integral in (3.1). We rewrite (3.1) in the form The approximation of (3.5) using the block-by-block method yield the following system:
3.3. Using Convolution Quadrature
In problems where the convolution kernel has a singularity, closed Newton’s cotes formulas fail in the computation of the solution because of the singularity. To overcome this problem we use convolution quadrature. An introduction to convolution quadrature is found in Lubich [15, 16] and he also gives a review of all aspects of convolution quadrature in [17].
We rewrite (3.1) in the form
Here the integral in (3.8) is approximated by convolution quadrature with a step size , where the convolution quadrature weights are the coefficients of the generating power series
Here, is the Laplace transform of the convolution kernel and is the quotient of the generating polynomials of a linear multistep method. We will use a convolution quadrature formula based on the order 2 BDF method, that is, we use , (see Lubich [17]). The approximation of (3.8) by the convolution quadrature method yields the following system:
4. Numerical Examples
In this section we compute examples of (3.1).
Example 4.1.
with , and .
Tables 1 and 2 show the blow-up results by the repeated trapezoidal rule and the block-by-block method, respectively, and Figure 1 shows the solution of (4.1) computed by the repeated trapezoidal rule.
(a) By PECE-IE
(b) By IMIE
(c) By ode45
(d) By ode15s
Example 4.2.
with and .
Tables 3 and 4 show the blow-up results by repeated trapezoidal rule and the block-by-block method, respectively, and Figure 2 shows the solution of (4.2) computed by trapezoidal rule.
(a) By PECE-IE
(b) By IMIE
(c) By ode45
(d) By ode15s
Example 4.3.
with and .
We use same parameters as in Ma [7] and compare the results.
Equation (4.3) has a convolution kernel with a singularity and hence we use convolution quadrature to solve it. We have
where is given by the generating power series (4.4), with and . .
Table 5 shows the blow-up results and Figure 3 shows the solution of (4.3).
(a) By PECE-IE
(b) By IMIE
(c) By ode45
(d) By ode15s
4.1. Discussion
From the results obtained, we observe that blow-up occurs earlier as the exponent increases as shown in Figures 1, 2, and 3. Comparing the the blow-up results computed when we use the repeated trapezoidal rule and the block-by-block method, we note a significant difference in the blow-up computed, which means the choice of quadrature rule has an effect on the blow-up time. In Example 4.3, where we use convolution quadrature we get results comparable to Ma’s results, (see [7]).
On the performance of the methods, we note that the IMIE method gives results which are much closer to those for ode45 which is of higher order than the other three. Thus based on our results we conclude that the adaptive IMIE method is superior than the adaptive PECE-IE method and ode15s.