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𝑢𝑡(𝑡,𝑥)Δ𝑢(𝑡,𝑥)=𝑓(𝑡,𝑥,𝑢),𝑡>0,𝑥Ω𝑑,𝑢(0,𝑥)=𝑢0(𝑥)0,𝑥Ω,𝑢(𝑡,𝑥)=0,𝑡>0,𝑥𝜕Ω,(1.1) 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 𝑢0(𝑥) 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, 𝑓(𝑡,𝑥,𝑢)=𝑡0𝑘(𝑡𝑠)𝑢𝑝(𝑠,𝑥)𝑑𝑠so that (1.1) becomes𝑢𝑡(𝑡,𝑥)𝑢𝑥𝑥(𝑡,𝑥)=𝑡0𝑘(𝑡𝑠)𝑢𝑝(𝑠,𝑥)𝑑𝑠,𝑡>0,𝑥(0,1).(1.2)

In this paper, we consider the Volterra integrodifferential equation (VIDE),𝑦(𝑡)=𝜆𝑦(𝑡)+𝑡0𝑘(𝑡𝑠)𝑔(𝑦(𝑠))𝑑𝑠,𝑡>0,𝑦(0)=𝑦0>0(1.3) with 𝜆<0 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 lim𝑡𝐾(𝑡)=, where 𝐾(𝑡)=𝑡0𝑘(𝑡)𝑑𝑠;(b)𝑔(𝑡) is nonnegative, nondecreasing and continuous for 𝑡>0, 𝑔0 for 𝑡0, and lim𝑦(𝑔(𝑦)/𝑦)=,

he proves that the solution of the VIDE (1.3) blows up if and only if𝛼𝑠𝑔(𝑠)1/𝛽<,(1.4) for any 𝛼>0, where 𝛽>0.

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:𝑦=𝑓(𝑡,𝑦),𝑡0.(2.1)

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𝑡𝑖0=𝑡0<𝑡1<𝑡2<,𝑖=0,1,2,3,(2.2) with 𝑡𝑖𝑡𝑖1=𝑖.

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 𝑆1 and 𝑆2. We approximate the local truncation error by computing the difference between the two solutions, that is, |𝑆1𝑆2|. If the error is smaller than the specified acceptable error, we move to the next step and accept 𝑆2 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 𝑆1 using a predictor-corrector method in which 𝑦𝑝 is obtained using explicit Euler’s method so that we have𝑆1=𝑦𝑖+1=𝑦𝑖+𝑖𝑓𝑡𝑖+1,𝑦𝑝.(2.3)

To compute 𝑆2 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 𝑆1. To compute 𝑆2 we use implicit midpoint method given by𝑦𝑖+1=𝑦𝑖+𝑖𝑓𝑡𝑖+𝑖2,𝑦𝑖+𝑦𝑖+12.(2.4)

First we compute 𝑦𝑖+1 using the implicit Euler method and substitute in (2.4) to get 𝑆2 that is, 𝑦𝑝=𝑦𝑖+𝑖𝑓𝑡𝑖+1,𝑦𝑖+1,(2.5)

then, 𝑆2=𝑦𝑖+1=𝑦𝑖+𝑖𝑓𝑡𝑖+𝑖2,𝑦𝑖+𝑦𝑝2.(2.6)

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𝑦(𝑡)=𝜆𝑦(𝑡)+𝑡0𝑘(𝑡𝑠)𝑦𝑝(𝑠)𝑑𝑠,𝑡>0,𝑦(0)=𝑦0>0(3.1) with 𝜆<0 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 [1517].

3.1. Using Repeated Trapezoidal Rule

We let𝑡0𝑡=0,𝑖+1=𝑡𝑖+𝑖𝑦,𝑖=0,1,2,3,,𝑛,(3.2)𝑡𝑖𝑡=𝜆𝑦𝑖+𝑖𝑗=1𝑡𝑗𝑡𝑗1𝑘𝑡𝑖𝑦𝑠𝑝(𝑠)𝑑𝑠,𝑖=1,2,3,,𝑛.(3.3)

The approximation of every integral in (3.3) by repeated trapezoidal rule will yield the following system (with 𝑦𝑖=𝑦(𝑡𝑖) and 𝑘𝑖,𝑗=𝑘(𝑡𝑖,𝑡𝑗))𝑦1=𝜆𝑦1+12𝑘1,0𝑦𝑝0+𝑘1,1𝑦𝑝1,𝑦2=𝜆𝑦2+22𝑘2,0𝑦𝑝0+𝑘2,1𝑦𝑝1+𝑘2,2𝑦𝑝2,𝑦3=𝜆𝑦3+32𝑘3,0𝑦𝑝0+𝑘3,1𝑦𝑝1+𝑘3,2𝑦𝑝2+𝑘3,3𝑦𝑝3,𝑦𝑛=𝜆𝑦𝑛+𝑛2𝑘𝑛,0𝑦𝑝0+𝑘𝑛,1𝑦𝑝1++𝑘𝑛,𝑛1𝑦𝑝𝑛1+𝑘𝑛,𝑛𝑦𝑝𝑛.(3.4)

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𝑦𝑡𝑖𝑡=𝜆𝑦𝑖+𝑡𝑖0𝑘𝑡𝑖𝑦𝑠𝑝(𝑠)𝑑𝑠,𝑡𝑖>0,𝑦(0)=𝑦0𝑡>0,(3.5)0𝑡=0,𝑖+1=𝑡𝑖+𝑖,𝑖=0,1,2,3,,𝑛.(3.6) The approximation of (3.5) using the block-by-block method yield the following system: 𝑦1=𝜆𝑦1+51𝑘121,0𝑦𝑝0+213𝑘1,1𝑦𝑝11𝑘121,2𝑦𝑝2,𝑦2=𝜆𝑦2+23𝑘2,0𝑦𝑝0+4𝑘2,1𝑦𝑝1+𝑘2,2𝑦𝑝2,𝑦2𝑛+1=𝜆𝑦2𝑛+1+2𝑛+13𝑘2𝑛+1,0𝑦𝑝0+4𝑘2𝑛+1,1𝑦𝑝1++𝑘2𝑛+1,2𝑛𝑦𝑝2𝑛+52𝑛+1𝑘122𝑛+1,2𝑛𝑦𝑝2𝑛+22𝑛+13𝑘2𝑛+1,2𝑛+1𝑦𝑝2n+12𝑛+1𝑘122𝑛+1,2𝑛+2𝑦𝑝2𝑛+2,𝑦2𝑛+2=𝜆𝑦2𝑛+2+2𝑛+23𝑘2𝑛+2,0𝑦𝑝0+4𝑘2𝑛+2,1𝑦𝑝1++4𝑘2𝑛+2,2𝑛+1𝑦𝑝2𝑛+1+𝑘2𝑛+2,2𝑛+2𝑦𝑝2𝑛+2.(3.7)

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𝑦𝑡𝑖𝑡=𝜆𝑦𝑖+𝑡𝑖0𝑘𝑡𝑖𝑦𝑠𝑝(𝑠)𝑑𝑠,𝑡𝑖>0,𝑦(0)=𝑦0𝑡>0,(3.8)0𝑡=0,𝑖+1=𝑡𝑖+𝑖,𝑖=0,1,2,3,,𝑛.(3.9)

Here the integral𝑡𝑖0𝑘𝑡𝑖𝑦𝑠𝑝(𝑠)𝑑𝑠(3.10) in (3.8) is approximated by convolution quadrature with a step size >0,𝑖𝑗=0𝑤𝑗𝑦𝑝𝑡𝑗,𝑖=1,2,3,,𝑛,(3.11) where the convolution quadrature weights 𝑤𝑗 are the coefficients of the generating power series𝑗=0𝑤𝑗𝜉𝑗=𝐹𝛿(𝜉).(3.12)

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 𝛿(𝜉)=(1𝜉)+(1/2)(1𝜉)2, (see Lubich [17]). The approximation of (3.8) by the convolution quadrature method yields the following system:𝑦1=𝜆𝑦1+𝑤0𝑦𝑝0+𝑤1𝑦𝑝1,𝑦2=𝜆𝑦2+𝑤0𝑦𝑝0+𝑤1𝑦𝑝1+𝑤2𝑦𝑝2,𝑦3=𝜆𝑦3+𝑤0𝑦𝑝0+𝑤1𝑦𝑝1+𝑤2𝑦𝑝2+𝑤3𝑦𝑝3,𝑦𝑛=𝜆𝑦𝑛+𝑤0𝑦𝑝0+𝑤1𝑦𝑝1++𝑤𝑛1𝑦𝑝𝑛1+𝑤𝑛𝑦𝑝𝑛.(3.13)

4. Numerical Examples

In this section we compute examples of (3.1).

Example 4.1. 𝑦(𝑡)=𝜆𝑦(𝑡)+𝑡0𝑦𝑝(𝑠)𝑑𝑠,𝑡>0,𝑦(0)=𝑦0>0,(4.1) with 𝜆=1, 𝑦(0)=𝑦0=3 and 𝑝=2,2.5,3.
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.

Example 4.2. 𝑦(𝑡)=𝜆𝑦(𝑡)+𝑡0𝑒𝑠𝑡𝑦𝑝(𝑠)𝑑𝑠,𝑡>0,𝑦(0)=𝑦0>0,(4.2) with 𝜆=1,𝑦(0)=𝑦0=3 and 𝑝=2,2.5,3.
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.

Example 4.3. 𝑦(𝑡)=𝜆𝑦+𝑡0(𝑡𝑠)𝛼𝑦𝑝(𝑠)𝑑𝑠,𝑡>0,𝑦(0)=𝑦0>0,(4.3) with 𝜆=1,𝛼=2/3,𝑦(0)=𝑦0=1 and 𝑝=2,2.5,3.
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 𝑗=0𝑤𝑗𝜉𝑗=𝐹𝛿(𝜉)1=Γ33221𝜉+𝜉221/3,(4.4) where 𝑤𝑗 is given by the generating power series (4.4), with 𝐹(𝑠)=Γ(1/3)𝑠1/3 and 𝛿(𝜉)=(1𝜉)+(1/2)(1𝜉)2. Γ(1/3)=2.6789385347.
Table 5 shows the blow-up results and Figure 3 shows the solution of (4.3).

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.