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

Authors: P. G. Dlamini; M. Khumalo
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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 atTb<∞ if the solution of the differential equation becomes unbounded at t=Tb. 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(1.1)ut(t,x)-Δu(t,x)=f(t,x,u),t>0,x∈Ω⊂Rd,u(0,x)=u0(x)≥0,x∈Ω,u(t,x)=0,t>0,x∈∂Ω,
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 u0(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)=∫0tk(t-s)up(s,x)dsso that (1.1) becomes(1.2)ut(t,x)-uxx(t,x)=∫0tk(t-s)up(s,x)ds,t>0,x∈(0,1).In this paper, we consider the Volterra integrodifferential equation (VIDE),(1.3)y′(t)=λy(t)+∫0tk(t-s)g(y(s))ds,t>0,y(0)=y0>0
with λ<0 and k(t-s) 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)
k(t) is an integrable positive function such that limt→∞K(t)=∞, where K(t)=∫0tk(t)ds;(b)
g(t) is nonnegative, nondecreasing and continuous for t>0, g≡0 for t≤0, and limy→∞(g(y)/y)=∞,he proves that the solution of the VIDE (1.3) blows up if and only if(1.4)∫α∞[sg(s)]1/β<∞,
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:(2.1)y′=f(t,y),t≥0.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(2.2){ti:0=t0<t1<t2<⋯},i=0,1,2,3,…
with ti-ti-1=hi.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, sayS1 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.


### 2.1. Adaptive PECE-Implicit Euler Method (PECE-IE)

This method is based on the implicit Euler method. We computeS1 using a predictor-corrector method in which yp is obtained using explicit Euler’s method so that we have(2.3)S1=yi+1=yi+hif(ti+1,yp).To computeS2 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 getS1. To compute S2 we use implicit midpoint method given by(2.4)yi+1=yi+hif(ti+hi2,yi+yi+12).First we computeyi+1 using the implicit Euler method and substitute in (2.4) to get S2 that is,
(2.5)yp=yi+hif(ti+1,yi+1),then,(2.6)S2=yi+1=yi+hif(ti+hi2,yi+yp2).


### 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.


## 2.1. Adaptive PECE-Implicit Euler Method (PECE-IE)

This method is based on the implicit Euler method. We computeS1 using a predictor-corrector method in which yp is obtained using explicit Euler’s method so that we have(2.3)S1=yi+1=yi+hif(ti+1,yp).To computeS2 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 getS1. To compute S2 we use implicit midpoint method given by(2.4)yi+1=yi+hif(ti+hi2,yi+yi+12).First we computeyi+1 using the implicit Euler method and substitute in (2.4) to get S2 that is,
(2.5)yp=yi+hif(ti+1,yi+1),then,(2.6)S2=yi+1=yi+hif(ti+hi2,yi+yp2).


## 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(3.1)y′(t)=λy(t)+∫0tk(t-s)yp(s)ds,t>0,y(0)=y0>0
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].


### 3.1. Using Repeated Trapezoidal Rule

We let(3.2)t0=0,ti+1=ti+hi,i=0,1,2,3,…,n,(3.3)y′(ti)=λy(ti)+∑j=1i∫tj-1tjk(ti-s)yp(s)ds,i=1,2,3,…,n.The approximation of every integral in (3.3) by repeated trapezoidal rule will yield the following system (with yi=y(ti) and ki,j=k(ti,tj))(3.4)y1′=λy1+h12[k1,0y0p+k1,1y1p],y2′=λy2+h22[k2,0y0p+k2,1y1p+k2,2y2p],y3′=λy3+h32[k3,0y0p+k3,1y1p+k3,2y2p+k3,3y3p],⋮yn′=λyn+hn2[kn,0y0p+kn,1y1p+⋯+kn,n-1yn-1p+kn,nynp].


### 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(3.5)y′(ti)=λy(ti)+∫0tik(ti-s)yp(s)ds,ti>0,y(0)=y0>0,(3.6)t0=0,ti+1=ti+hi,i=0,1,2,3,…,n.
The approximation of (3.5) using the block-by-block method yield the following system:
(3.7)y1′=λy1+5h112k1,0y0p+2h13k1,1y1p-h112k1,2y2p,y2′=λy2+h23[k2,0y0p+4k2,1y1p+k2,2y2p],⋮y2n+1′=λy2n+1+h2n+13[k2n+1,0y0p+4k2n+1,1y1p+⋯+k2n+1,2ny2np]+5h2n+112k2n+1,2ny2np+2h2n+13k2n+1,2n+1y2n+1p-h2n+112k2n+1,2n+2y2n+2p,y2n+2′=λy2n+2+h2n+23[k2n+2,0y0p+4k2n+2,1y1p+⋯+4k2n+2,2n+1y2n+1p+k2n+2,2n+2y2n+2p].


### 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(3.8)y′(ti)=λy(ti)+∫0tik(ti-s)yp(s)ds,ti>0,y(0)=y0>0,(3.9)t0=0,ti+1=ti+hi,i=0,1,2,3,…,n.Here the integral(3.10)∫0tik(ti-s)yp(s)ds
in (3.8) is approximated by convolution quadrature with a step size h>0,(3.11)∑j=0iwjyp(tj),i=1,2,3,…,n,
where the convolution quadrature weights wj are the coefficients of the generating power series(3.12)∑j=0∞wjξj=F(δ(ξ)h).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:(3.13)y1′=λy1+w0y0p+w1y1p,y2′=λy2+w0y0p+w1y1p+w2y2p,y3′=λy3+w0y0p+w1y1p+w2y2p+w3y3p,⋮yn′=λyn+w0y0p+w1y1p+⋯+wn-1yn-1p+wnynp.


## 3.1. Using Repeated Trapezoidal Rule

We let(3.2)t0=0,ti+1=ti+hi,i=0,1,2,3,…,n,(3.3)y′(ti)=λy(ti)+∑j=1i∫tj-1tjk(ti-s)yp(s)ds,i=1,2,3,…,n.The approximation of every integral in (3.3) by repeated trapezoidal rule will yield the following system (with yi=y(ti) and ki,j=k(ti,tj))(3.4)y1′=λy1+h12[k1,0y0p+k1,1y1p],y2′=λy2+h22[k2,0y0p+k2,1y1p+k2,2y2p],y3′=λy3+h32[k3,0y0p+k3,1y1p+k3,2y2p+k3,3y3p],⋮yn′=λyn+hn2[kn,0y0p+kn,1y1p+⋯+kn,n-1yn-1p+kn,nynp].


## 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(3.5)y′(ti)=λy(ti)+∫0tik(ti-s)yp(s)ds,ti>0,y(0)=y0>0,(3.6)t0=0,ti+1=ti+hi,i=0,1,2,3,…,n.
The approximation of (3.5) using the block-by-block method yield the following system:
(3.7)y1′=λy1+5h112k1,0y0p+2h13k1,1y1p-h112k1,2y2p,y2′=λy2+h23[k2,0y0p+4k2,1y1p+k2,2y2p],⋮y2n+1′=λy2n+1+h2n+13[k2n+1,0y0p+4k2n+1,1y1p+⋯+k2n+1,2ny2np]+5h2n+112k2n+1,2ny2np+2h2n+13k2n+1,2n+1y2n+1p-h2n+112k2n+1,2n+2y2n+2p,y2n+2′=λy2n+2+h2n+23[k2n+2,0y0p+4k2n+2,1y1p+⋯+4k2n+2,2n+1y2n+1p+k2n+2,2n+2y2n+2p].


## 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(3.8)y′(ti)=λy(ti)+∫0tik(ti-s)yp(s)ds,ti>0,y(0)=y0>0,(3.9)t0=0,ti+1=ti+hi,i=0,1,2,3,…,n.Here the integral(3.10)∫0tik(ti-s)yp(s)ds
in (3.8) is approximated by convolution quadrature with a step size h>0,(3.11)∑j=0iwjyp(tj),i=1,2,3,…,n,
where the convolution quadrature weights wj are the coefficients of the generating power series(3.12)∑j=0∞wjξj=F(δ(ξ)h).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:(3.13)y1′=λy1+w0y0p+w1y1p,y2′=λy2+w0y0p+w1y1p+w2y2p,y3′=λy3+w0y0p+w1y1p+w2y2p+w3y3p,⋮yn′=λyn+w0y0p+w1y1p+⋯+wn-1yn-1p+wnynp.


## 4. Numerical Examples

In this section we compute examples of (3.1).Example 4.1.
(4.1)y′(t)=λy(t)+∫0typ(s)ds,t>0,y(0)=y0>0,
with λ=-1, y(0)=y0=3 and p=2,2.5,3.
Tables1 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.Table 1
Blow-up time for (4.1) by repeated trapezoidal rule.
pPECE-IEIMIEode45ode15s22.0582572.0649952.0650022.0649952.51.2538121.2591851.2591871.25917930.7641480.7682750.7682770.768271Table 2
Blow-up time for (4.1) by block-by-block method.
pPECE-IEIMIEode45ode15s21.8939821.8999681.8999701.8999672.51.2164451.2212891.2212911.22128430.6639210.6673930.6673950.667393Numerical solution of (4.1) by repeated trapezoidal rule.
(a)
By PECE-IE(b)
By IMIE(c)
By ode45(d)
By ode15sExample 4.2.
(4.2)y′(t)=λy(t)+∫0tes-typ(s)ds,t>0,y(0)=y0>0,
with λ=-1,y(0)=y0=3 and p=2,2.5,3.
Tables3 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.Table 3
Blow-up time for (4.2) by repeated trapezoidal rule.
pPECE-IEIMIEode45ode15s23.1418133.1575333.1575393.1575102.51.5839001.5927921.5928001.59272830.8785940.8842380.8842390.884229Table 4
Blow-up time for (4.2) by block-by-block method.
pPECE-IEIMIEode45ode15s22.7699462.7790282.7790332.7790192.51.5962931.6001801.6001821.60016930.9606310.9647310.9647320.964730Numerical solution of (4.2) by repeated trapezoidal rule.
(a)
By PECE-IE(b)
By IMIE(c)
By ode45(d)
By ode15sExample 4.3.
(4.3)y′(t)=λy+∫0t(t-s)-αyp(s)ds,t>0,y(0)=y0>0,
with λ=-1,α=2/3,y(0)=y0=1 and p=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
(4.4)∑j=0∞wjξj=F(δ(ξ)h)=Γ(13)[32h-2hξ+12hξ2]-1/3,
where wj is given by the generating power series (4.4), with F(s)=Γ(1/3)s-1/3 and δ(ξ)=(1-ξ)+(1/2)(1-ξ)2. Γ(1/3)=2.6789385347.
Table5 shows the blow-up results and Figure 3 shows the solution of (4.3).Table 5
Blow-up time for (4.3).
pPECE-IEIMIEode45ode15s21.5749121.5801571.5801591.5801532.50.8820660.8854170.8854180.88541330.6033340.6058260.6058270.605821Numerical 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 exponentp 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.


## 4.1. Discussion

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



---

*Source: 878497-2012-05-20.xml*