On the Computation of Blow-Up Solutions for Nonlinear Volterra Integrodifferential Equations

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.


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 T b < ∞ if the solution of the differential equation becomes unbounded at t T b .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 Δu t, x is referred to as the diffusion term and f t, x, u as the reaction term.It is well known from blow-up theories that for sufficiently large initial function u 0 x 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, f t, x, u t 0 k t − s u p s, x ds so that 1.1 becomes In this paper, we consider the Volterra integrodifferential equation

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 midpointimplicit 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 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 S1 and S2.We approximate the local truncation error by computing the difference between the two solutions, that is, |S1 − S2|.If the error is smaller than the specified acceptable error, we move to the next step and accept S2 as the solution for the current step.Otherwise, we set a new stepsize and redo the current step.

Adaptive PECE-Implicit Euler Method (PECE-IE)
This method is based on the implicit Euler method.We compute S1 using a predictorcorrector method in which y p is obtained using explicit Euler's method so that we have To compute S2 we use Newton's method as the solver to deal with the implicit nature of the implicit Euler method.

Adaptive Implicit Midpoint-Implicit Euler Method (IMIE)
We use the implicit Euler method with Newton's method as the solver to get S1.To compute S2 we use implicit midpoint method given by First we compute y i 1 using the implicit Euler method and substitute in 2.4 to get S2 that is, 2.6

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.

Numerical Computation
We consider y t λy t t 0 k t − s y p s ds, t > 0, y 0 y 0 > 0 3.1 with λ < 0 and k t − s 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 .

Using Repeated Trapezoidal Rule
We let t 0 0, 3.4

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 y t i λy t i t i 0 k t i − s y p s ds, t i > 0, y 0 y 0 > 0, 3.5 t 0 0, t i 1 t i h i , i 0, 1, 2, 3, . . ., n.

3.6
The approximation of 3.5 using the block-by-block method yield the following system: 3.7

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 .

Mathematical Problems in Engineering
We rewrite 3. where the convolution quadrature weights w j are the coefficients of the generating power series ∞ j 0 w j ξ j F δ ξ h .

3.12
Here, F s is the Laplace transform of the convolution kernel k t 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:

Numerical Examples
In this section we compute examples of 3.1 .y p s ds, t > 0, y 0 y 0 > 0, 4.1 with λ −1, y 0 y 0 3 and p 2, 2.5, 3. Tables 1 and 2 show the blow-up results by the repeated trapezoidal rule and the blockby-block method, respectively, and Figure 1 shows the solution of 4.1 computed by the repeated trapezoidal rule.3 and 4 show the blow-up results by repeated trapezoidal rule and the block-byblock method, respectively, and Figure 2 shows the solution of 4.2 computed by trapezoidal rule.

Discussion
From the results obtained, we observe that blow-up occurs earlier as the exponent p 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

Figure 1 :
Figure 1: Numerical solution of 4.1 by repeated trapezoidal rule.

Figure 2 :
Figure 2: Numerical solution of 4.2 by repeated trapezoidal rule.
The approximation of every integral in 3.3 by repeated trapezoidal rule will yield the following system with y i y t i and k i,j k t i , t j k t i − s y p s ds, i 1, 2, 3, . .., n.3.3

Table 1 :
Blow-up time for 4.1 by repeated trapezoidal rule.

Table 2 :
Blow-up time for 4.1 by block-by-block method.

Table 3 :
Blow-up time for 4.2 by repeated trapezoidal rule.

Table 4 :
Blow-up time for 4.2 by block-by-block method.