# On the Computation of Blow-up Solutions for Semilinear ODEs and Parabolic PDEs

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

---

## Abstract

We introduce an adaptive numerical method for computing blow-up solutions for
ODEs and well-known reaction-diffusion equations. The method is based on the
implicit midpoint method and the implicit Euler method. We demonstrate that
the method produces superior results to the adaptive PECE-implicit method
and the MATLAB solver of comparable order.


---

## Body

## 1. Introduction

Reaction-diffusion equations model a wide range of problems in physics, biology, and chemistry. They explain how the concentration of one or more substances distributed in space changes under the influence of two processes: chemical reactions and diffusion. These substances can be basic particles in physics, bacteria, molecules, cells, and so forth. The substances reside in a regionΩ⊂ℝd,d≥1.The reaction-diffusion equation is a semilinear parabolic partial differential 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∈∂ΩEquation (1.1) can be viewed as a heat conduction problem, where u(x,t) is the temperature of a substance in a bounded domain Ω⊂ℝd and f(t,x,u) represents a heat source. Δu(t,x) is referred to as the diffusion term and f(t,x,u) as the reaction term. In this case, convection does not take place, so f does not depend on ∇u.For sufficiently large initial functionu0(x) the solution of (1.1) will blowup in finite time. Blow-up occurs when the solution of the partial differential equation ceases to exist in finite time; that is, there is Tb<∞ (blow-up time) so that
(1.2)lim⁡t→Tb-‖u(t,⋅)‖=+∞.Bebernes and Eberly [1] state that a necessary condition for blow-up in finite time is if (1.3)∫u0∞f-1(s)ds<∞.
Kaplan [2] showed that for convex source terms f=f(u) satisfying (1.3) diffusion cannot prevent blow-up if the initial state is large enough. In most papers, blow-up properties are discussed in the case where the nonlinear term in (1.1) is of the form f=f(u) where Ω⊂ℝd is bounded and t∈(0,T), where T is finite. In most of the work the case where d=1 has been studied. However, more recently, Brunner et al. [3] studied the numerical solution of blow-up problems within the context of unbounded domains.Stuart and Floater [4] 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. Bandle and Brunner [5] present a survey of the theory and the numerical analysis of blow-up solutions for quasilinear reaction-diffusion equations. Budd et al. [6] proposed the use of moving mesh partial differential equation (MMPDE) methods for solving (1.1). In this method the function u(x,t) is discretized to give the solution values ui(t) defined on a moving mesh xi(t),i=0,…,N. A more general study of the MMPDE is presented in [7] and the references therein. More recently Ma et al. [8, 9] have used the moving mesh methods to numerically study blow-up in nonlocal reaction diffusion equations and partial integrodifferential equations in general.In this paper we will use the method of lines (MOLs) to solve (1.1). In this method the PDE is discretized in space, which leads to a system of ODEs with initial conditions. The numerical solution can then be obtained by solving the ODE initial value problems (see [10]). We introduce an adaptive method based on the implicit midpoint method and the implicit Euler method to solve the resulting system of ODEs. More work has been done on PsDEs with autonomous nonlinear reaction term. In this work we also give numerical results for PDEs with nonautonomous nonlinear term.


## 2. Description of Methods Used

In this work we use variable step methods to compute the blow-up time. We use one-point collocation methods and compare the results with MATLAB solvers ode45 and ode15s. 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 adjust the stepsize in accordance with |S1-S2|.


### 2.1. One-Point Collocation Methods

One-point collocation methods are a family of methods of the form(2.1)yn+1=yn+hnf(tn+c1hn,(1-c1)yn+c1yn+1),
where c1∈[0,1]. The specified cases are (i)
c1=0 corresponds to explicit Euler method;(ii)
c1=1/2 corresponds to implicit midpoint method;(iii)
c1=1 corresponds to implicit Euler method.Note that all the one-point collocation methods are of order 1; however, the implicit midpoint method (c1=1/2) achieves order 2 local superconvergence (see [11]).


### 2.2. Adaptive PECE-Implicit Euler Method

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.2)S1=yn+1=yn+hnf(tn+1,yp),To getS2, we use Newton's method as the solver to deal with the implicit nature of the implicit Euler method.


### 2.3. Adaptive Implicit Midpoint-Implicit Euler Method

We use the implicit Euler method with Newton's method as the solver to getS1. To get S2 we use midpoint Euler method given by(2.3)yn+1=yn+hnf(tn+hn2,yn+yn+12).First we getyn+1 using the implicit Euler method and substitute in (2.3) to get S2, that is,(2.4)yp=yn+hnf(tn+1,yn+1),
then,
(2.5)S2=yn+1=yn+hnf(tn+hn2,yn+yp2).


### 2.4. 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 with the adaptive implicit midpoint-implicit Euler method we are introducing.


## 2.1. One-Point Collocation Methods

One-point collocation methods are a family of methods of the form(2.1)yn+1=yn+hnf(tn+c1hn,(1-c1)yn+c1yn+1),
where c1∈[0,1]. The specified cases are (i)
c1=0 corresponds to explicit Euler method;(ii)
c1=1/2 corresponds to implicit midpoint method;(iii)
c1=1 corresponds to implicit Euler method.Note that all the one-point collocation methods are of order 1; however, the implicit midpoint method (c1=1/2) achieves order 2 local superconvergence (see [11]).


## 2.2. Adaptive PECE-Implicit Euler Method

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.2)S1=yn+1=yn+hnf(tn+1,yp),To getS2, we use Newton's method as the solver to deal with the implicit nature of the implicit Euler method.


## 2.3. Adaptive Implicit Midpoint-Implicit Euler Method

We use the implicit Euler method with Newton's method as the solver to getS1. To get S2 we use midpoint Euler method given by(2.3)yn+1=yn+hnf(tn+hn2,yn+yn+12).First we getyn+1 using the implicit Euler method and substitute in (2.3) to get S2, that is,(2.4)yp=yn+hnf(tn+1,yn+1),
then,
(2.5)S2=yn+1=yn+hnf(tn+hn2,yn+yp2).


## 2.4. 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 with the adaptive implicit midpoint-implicit Euler method we are introducing.


## 3. Blow-up for ODEs

We will first consider this simple case(3.1)y′(t)=λy(t)+byp(t),t>0,y(0)=y0>0,
with λ<0 and b>0.


### 3.1. Analytic Solution

Theorem 3.1.
Given the system (3.1), its solution will blowup in finite time for
(3.2)y0>(-bλ)1/(1-p)
and at
(3.3)Tb=1λ(1-p)ln⁡bλy01-p+b
(see Brunner [11]).Proof.
Note that (3.1) is a Bernoulli equation, whose solution is
(3.4)y(t)=[-bλ+(y0(1-p)+bλ)eλ(1-p)t]1/(1-p).
Then, y→∞ when
(3.5)-bλ+(y0(1-p)+bλ)eλ(1-p)t=0.
Thus
(3.6)Tb=1λ(1-p)ln⁡bλy01-p+b.
Tb is finite if
(3.7)λy01-p+b>0.
Thus
(3.8)y0>(-bλ)1/(1-p),
as required.We compute the blow-up time forλ=-1,b=1,y0=2, and p=1.1,1.5,2,3. The results are shown in Table 1, and a graph showing the analytic solution of (3.1) is shown in Figure 1.Table 1
Blow-up solutions of (3.1).
pBlow-up time (Tb)1.127.035551.52.45589420.69314730.143841Figure 1
Analytic solution of (3.1).


### 3.2. Numerical Computation

We solve (3.1) with λ=-1, b=1, y0=2, and p=1.1,1.5,2,3. Tables 2 and 3 show the blow-up times and the errors of each method, respectively, and Figure 2 shows the graphs of the solution of (3.1). The tolerance used for the computations is 1e-6.Table 2
Blow-up time for (3.1).
pPECE-implicit EulerMidpoint-implicit Eulerode45ode15s1.126.8973327.0353027.0355527.035281.52.4421152.4558862.4558952.45585720.68897820.69314510.69314740.693129530.14290360.14384020.14384110.1438358Table 3
Error for (3.1).
pPECE-implicit EulerMidpoint-implicit Eulerode45ode15s1.10.8617800.00025000.0002701.50.0137790.00000800.00003820.0041690.00000200.00001830.0009370.00000100.000005Numerical solution of (3.1).
(a)
By PECE-implicit Euler(b)
By midpoint-implicit Euler(c)
By ode45(d)
By ode15s


### 3.3. Discussion

The blow-up results for the different methods are very close to the analytic value as shown in Table2. From Table 3, we see that ode45, which is of higher order than the other three methods, gives the best results than the other three methods. The adaptive implicit midpoint-implicit Euler gives a better result than the other two methods of comparable order, that is, adaptive PECE-implicit Euler method and ode15s. The adaptive PECE-implicit Euler method gives a quite large error and requires a very small tolerance to get a result which is close to the exact value. From the results, we observe that as the value of p in (3.1) is increased the blow-up time occurs earlier and is much later for values of p much closer to 1.We seek to determine whether the performance of the methods is the same in the case where we have a reaction-diffusion equation.


## 3.1. Analytic Solution

Theorem 3.1.
Given the system (3.1), its solution will blowup in finite time for
(3.2)y0>(-bλ)1/(1-p)
and at
(3.3)Tb=1λ(1-p)ln⁡bλy01-p+b
(see Brunner [11]).Proof.
Note that (3.1) is a Bernoulli equation, whose solution is
(3.4)y(t)=[-bλ+(y0(1-p)+bλ)eλ(1-p)t]1/(1-p).
Then, y→∞ when
(3.5)-bλ+(y0(1-p)+bλ)eλ(1-p)t=0.
Thus
(3.6)Tb=1λ(1-p)ln⁡bλy01-p+b.
Tb is finite if
(3.7)λy01-p+b>0.
Thus
(3.8)y0>(-bλ)1/(1-p),
as required.We compute the blow-up time forλ=-1,b=1,y0=2, and p=1.1,1.5,2,3. The results are shown in Table 1, and a graph showing the analytic solution of (3.1) is shown in Figure 1.Table 1
Blow-up solutions of (3.1).
pBlow-up time (Tb)1.127.035551.52.45589420.69314730.143841Figure 1
Analytic solution of (3.1).


## 3.2. Numerical Computation

We solve (3.1) with λ=-1, b=1, y0=2, and p=1.1,1.5,2,3. Tables 2 and 3 show the blow-up times and the errors of each method, respectively, and Figure 2 shows the graphs of the solution of (3.1). The tolerance used for the computations is 1e-6.Table 2
Blow-up time for (3.1).
pPECE-implicit EulerMidpoint-implicit Eulerode45ode15s1.126.8973327.0353027.0355527.035281.52.4421152.4558862.4558952.45585720.68897820.69314510.69314740.693129530.14290360.14384020.14384110.1438358Table 3
Error for (3.1).
pPECE-implicit EulerMidpoint-implicit Eulerode45ode15s1.10.8617800.00025000.0002701.50.0137790.00000800.00003820.0041690.00000200.00001830.0009370.00000100.000005Numerical solution of (3.1).
(a)
By PECE-implicit Euler(b)
By midpoint-implicit Euler(c)
By ode45(d)
By ode15s


## 3.3. Discussion

The blow-up results for the different methods are very close to the analytic value as shown in Table2. From Table 3, we see that ode45, which is of higher order than the other three methods, gives the best results than the other three methods. The adaptive implicit midpoint-implicit Euler gives a better result than the other two methods of comparable order, that is, adaptive PECE-implicit Euler method and ode15s. The adaptive PECE-implicit Euler method gives a quite large error and requires a very small tolerance to get a result which is close to the exact value. From the results, we observe that as the value of p in (3.1) is increased the blow-up time occurs earlier and is much later for values of p much closer to 1.We seek to determine whether the performance of the methods is the same in the case where we have a reaction-diffusion equation.


## 4. Reaction-Diffusion Equation: One Space Dimension

In this section we compute the blow-up time for a one space dimension reaction-diffusion equation with an autonomous reaction term and a nonautonomous reaction term.


### 4.1. Autonomous Reaction Term

We solve the system (1.1) with the autonomous reaction term f=up(t,x), where p>1 and the domain Ω is just the real line, that is, d=1 and Ω=(0,1). The system becomes(4.1)ut(t,x)-uxx(t,x)=up(t,x),t>0,x∈(0,1),(p>1),u(0,x)=u0(x)≥0,x∈(0,1),u(t,x)=0,t>0,x∈∂Ω.We use the method of lines (MOLs) to discretize (4.1) in space. For the spatial discretization, we choose a uniform mesh Dh:={xm:0=x0<x1<⋯<xM+1=1} (with xm=mh) on Ω and replace uxx(t,xm)(1≤m≤M) by the standard central difference approximation. We use the function Asin(πx) as the initial function with different values of A>0. We get the following system of ODEs for Um(t)≈u(t,xm)(1≤m≤M):(4.2)Um′(t)=Um+1(t)-2Um(t)+Um-1(t)h2+Ump(t),(1≤m≤M),U0(t)=UM+1(t)=0,Um(0)=Asin⁡(πxm).We solve the system (4.2) with p=2 and h=(1-0)/M. Tables 4, 5, and 6 show the blow-up results obtained with M=50,100, and 200, respectively, and with different values of A in the initial function Asin(πx). Figures 3, 4, 5, and 6 show the graphs of the solution of (4.1) for A=10 and A=12.Table 4
Blow-up time for (4.1) with M=50.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s8.0No blow-upNo blow-upNo blow-upNo blow-up10.0No blow-upNo blow-upNo blow-upNo blow-up10.5No blow-upNo blow-upNo blow-upNo blow-up11.0No blow-upNo blow-upNo blow-upNo blow-up11.1No blow-upNo blow-upNo blow-upNo blow-up11.20.61840.62040.62250.619811.30.46110.46290.46330.462511.40.40190.40360.40380.403311.50.36490.36660.36680.366312.00.27290.27450.27460.274312.50.22800.22950.22960.229413.00.19890.20050.20040.200213.50.17770.17910.17910.179014.00.16130.16260.16270.162614.50.14800.14930.14940.149315.00.13700.13830.13840.138315.50.12770.12900.12900.129016.00.11970.12090.12090.1209Table 5
Blow-up time for (4.1) with M=100.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s8.0No blow-upNo blow-upNo blow-upNo blow-up10.0No blow-upNo blow-upNo blow-upNo blow-up10.5No blow-upNo blow-upNo blow-upNo blow-up11.0No blow-upNo blow-upNo blow-upNo blow-up11.1No blow-upNo blow-upNo blow-upNo blow-up11.2No blow-upNo blow-upNo blow-upNo blow-up11.3No blow-upNo blow-upNo blow-upNo blow-up11.40.53220.53300.53400.532511.50.43230.43310.43350.432812.00.29230.29240.29250.292212.50.23820.23890.23900.238813.00.20620.20630.20630.206213.50.18260.18330.18330.183214.00.16510.16580.16580.165614.50.15110.15180.15180.151615.00.13960.14030.14030.140215.50.12990.13060.13060.130516.00.12160.12220.12230.1222Table 6
Blow-up time for (4.1) with M=200.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s8.0No blow-upNo blow-upNo blow-upNo blow-up10.0No blow-upNo blow-upNo blow-upNo blow-up10.5No blow-upNo blow-upNo blow-upNo blow-up11.0No blow-upNo blow-upNo blow-upNo blow-up11.1No blow-upNo blow-upNo blow-upNo blow-up11.2No blow-upNo blow-upNo blow-upNo blow-up11.3No blow-upNo blow-upNo blow-upNo blow-up11.4No blow-upNo blow-upNo-blow-upNo blow-up11.50.50740.50770.50850.506712.00.30340.30370.30380.303412.50.24410.24440.24440.244213.00.20940.20970.20970.209513.50.18530.18560.18560.185514.00.16720.16750.16750.167314.50.15280.15310.15310.153015.00.14100.14130.14130.141215.50.13110.13140.13140.131316.00.12270.12290.12290.1229Numerical solution of (4.1) obtained using PECE-implicit Euler.
(a)
withA=10(b)
withA=12Numerical solution of (4.1) obtained using midpoint-implicit Euler.
(a)
WithA=10(b)
WithA=12Numerical solution of (4.1) obtained using ode45.
(a)
WithA=10(b)
WithA=12Numerical solution of (4.1) obtained using ode15s.
(a)
WithA=10(b)
WithA=12


### 4.2. Nonautonomous Reaction Term

We now solve the system (1.1) with the non-autonomous reaction term f=tkxrup(t,x), where p>1 and the domain Ω is just the real line, that is, d=1 and Ω=(0,1). The system becomes(4.3)ut(t,x)-uxx(t,x)=tkxrup(t,x),t>0,x∈(0,1),(p>1),u(0,x)=u0(x)≥0,x∈(0,1),u(t,x)=0,t>0,x∈∂Ω.As in (4.1) we use the method of lines to discretize (4.3) to obtain the following system:
(4.4)Um′(t)=Um+1(t)-2Um(t)+Um-1(t)h2+tkxmrUmp(t),(1≤m≤M),U0(t)=UM+1(t)=0,Um(0)=Asin⁡⁡(πxm).We solve (4.4) for k=0,1 and r=1,2, Tables 7 and 8 show the blow-up times obtained from the different values of k and r.Table 7
Blow-up time for (4.1) with M=100 and k=0.
Blow-up time(Tb)rAPECE-implicit EulerMidpoint-implicit Eulerode45ode15s121.0No blow-upNo blow-upNo blow-upNo blow-up21.5No blow-upNo blow-upNo blow-upNo blow-up22.00.46700.46790.46810.467522.50.33050.33120.33130.331023.00.27800.27870.27870.278523.50.24540.24610.24610.246024.00.22210.22280.22280.2227236.0No blow-upNo blow-upNo blow-upNo blow-up36.5No blow-upNo blow-upNo blow-upNo blow-up37.00.36730.36810.36810.367937.50.30380.30450.30450.304338.00.26830.26900.26900.268838.50.24380.24440.24450.244339.00.22520.22580.22590.2257Table 8
Blow-up time for (4.1) with M=100 and k=1.
Blow-up time(Tb)rAPECE-implicit EulerMidpoint-implicit Eulerode45ode15s1215No blow-upNo blow-upNo blow-upNo blow-up2160.80410.83670.83670.84102170.68720.69810.69810.69562180.63100.63780.63780.63842190.59360.59870.59870.59902200.56550.56960.56970.56842370No blow-upNo blow-upNo blow-upNo blow-up3710.80540.86710.86710.87543720.69490.71260.71260.71423730.64270.65330.65330.65413740.60800.61570.61570.61623750.58190.58810.58810.5884


### 4.3. Discussion

We observe that as we increase the amplitude of the initial function,A, the blow-up time tends to occur earlier. We also observe that for smaller values of A, there is no blow-up.Considering the reaction-diffusion equation, and comparing the autonomous against the non-autonomous reaction term cases, we observe that the introduction of the non-autonomous term ensures that a much larger amplitude,A, in the initial function, is required for blow-up to occur. We also note that increasing k for fixed r or increasing r for fixed k increases the minimum amplitude for blow-up to occur.On the performance of the methods, we note a similar trend to what we observed in the ODE case. The adaptive implicit midpoint-implicit Euler method gives results that are significantly superior to the adaptive PECE-implicit Euler method and ode15s. In fact, its performance is comparable to ode45.


## 4.1. Autonomous Reaction Term

We solve the system (1.1) with the autonomous reaction term f=up(t,x), where p>1 and the domain Ω is just the real line, that is, d=1 and Ω=(0,1). The system becomes(4.1)ut(t,x)-uxx(t,x)=up(t,x),t>0,x∈(0,1),(p>1),u(0,x)=u0(x)≥0,x∈(0,1),u(t,x)=0,t>0,x∈∂Ω.We use the method of lines (MOLs) to discretize (4.1) in space. For the spatial discretization, we choose a uniform mesh Dh:={xm:0=x0<x1<⋯<xM+1=1} (with xm=mh) on Ω and replace uxx(t,xm)(1≤m≤M) by the standard central difference approximation. We use the function Asin(πx) as the initial function with different values of A>0. We get the following system of ODEs for Um(t)≈u(t,xm)(1≤m≤M):(4.2)Um′(t)=Um+1(t)-2Um(t)+Um-1(t)h2+Ump(t),(1≤m≤M),U0(t)=UM+1(t)=0,Um(0)=Asin⁡(πxm).We solve the system (4.2) with p=2 and h=(1-0)/M. Tables 4, 5, and 6 show the blow-up results obtained with M=50,100, and 200, respectively, and with different values of A in the initial function Asin(πx). Figures 3, 4, 5, and 6 show the graphs of the solution of (4.1) for A=10 and A=12.Table 4
Blow-up time for (4.1) with M=50.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s8.0No blow-upNo blow-upNo blow-upNo blow-up10.0No blow-upNo blow-upNo blow-upNo blow-up10.5No blow-upNo blow-upNo blow-upNo blow-up11.0No blow-upNo blow-upNo blow-upNo blow-up11.1No blow-upNo blow-upNo blow-upNo blow-up11.20.61840.62040.62250.619811.30.46110.46290.46330.462511.40.40190.40360.40380.403311.50.36490.36660.36680.366312.00.27290.27450.27460.274312.50.22800.22950.22960.229413.00.19890.20050.20040.200213.50.17770.17910.17910.179014.00.16130.16260.16270.162614.50.14800.14930.14940.149315.00.13700.13830.13840.138315.50.12770.12900.12900.129016.00.11970.12090.12090.1209Table 5
Blow-up time for (4.1) with M=100.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s8.0No blow-upNo blow-upNo blow-upNo blow-up10.0No blow-upNo blow-upNo blow-upNo blow-up10.5No blow-upNo blow-upNo blow-upNo blow-up11.0No blow-upNo blow-upNo blow-upNo blow-up11.1No blow-upNo blow-upNo blow-upNo blow-up11.2No blow-upNo blow-upNo blow-upNo blow-up11.3No blow-upNo blow-upNo blow-upNo blow-up11.40.53220.53300.53400.532511.50.43230.43310.43350.432812.00.29230.29240.29250.292212.50.23820.23890.23900.238813.00.20620.20630.20630.206213.50.18260.18330.18330.183214.00.16510.16580.16580.165614.50.15110.15180.15180.151615.00.13960.14030.14030.140215.50.12990.13060.13060.130516.00.12160.12220.12230.1222Table 6
Blow-up time for (4.1) with M=200.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s8.0No blow-upNo blow-upNo blow-upNo blow-up10.0No blow-upNo blow-upNo blow-upNo blow-up10.5No blow-upNo blow-upNo blow-upNo blow-up11.0No blow-upNo blow-upNo blow-upNo blow-up11.1No blow-upNo blow-upNo blow-upNo blow-up11.2No blow-upNo blow-upNo blow-upNo blow-up11.3No blow-upNo blow-upNo blow-upNo blow-up11.4No blow-upNo blow-upNo-blow-upNo blow-up11.50.50740.50770.50850.506712.00.30340.30370.30380.303412.50.24410.24440.24440.244213.00.20940.20970.20970.209513.50.18530.18560.18560.185514.00.16720.16750.16750.167314.50.15280.15310.15310.153015.00.14100.14130.14130.141215.50.13110.13140.13140.131316.00.12270.12290.12290.1229Numerical solution of (4.1) obtained using PECE-implicit Euler.
(a)
withA=10(b)
withA=12Numerical solution of (4.1) obtained using midpoint-implicit Euler.
(a)
WithA=10(b)
WithA=12Numerical solution of (4.1) obtained using ode45.
(a)
WithA=10(b)
WithA=12Numerical solution of (4.1) obtained using ode15s.
(a)
WithA=10(b)
WithA=12


## 4.2. Nonautonomous Reaction Term

We now solve the system (1.1) with the non-autonomous reaction term f=tkxrup(t,x), where p>1 and the domain Ω is just the real line, that is, d=1 and Ω=(0,1). The system becomes(4.3)ut(t,x)-uxx(t,x)=tkxrup(t,x),t>0,x∈(0,1),(p>1),u(0,x)=u0(x)≥0,x∈(0,1),u(t,x)=0,t>0,x∈∂Ω.As in (4.1) we use the method of lines to discretize (4.3) to obtain the following system:
(4.4)Um′(t)=Um+1(t)-2Um(t)+Um-1(t)h2+tkxmrUmp(t),(1≤m≤M),U0(t)=UM+1(t)=0,Um(0)=Asin⁡⁡(πxm).We solve (4.4) for k=0,1 and r=1,2, Tables 7 and 8 show the blow-up times obtained from the different values of k and r.Table 7
Blow-up time for (4.1) with M=100 and k=0.
Blow-up time(Tb)rAPECE-implicit EulerMidpoint-implicit Eulerode45ode15s121.0No blow-upNo blow-upNo blow-upNo blow-up21.5No blow-upNo blow-upNo blow-upNo blow-up22.00.46700.46790.46810.467522.50.33050.33120.33130.331023.00.27800.27870.27870.278523.50.24540.24610.24610.246024.00.22210.22280.22280.2227236.0No blow-upNo blow-upNo blow-upNo blow-up36.5No blow-upNo blow-upNo blow-upNo blow-up37.00.36730.36810.36810.367937.50.30380.30450.30450.304338.00.26830.26900.26900.268838.50.24380.24440.24450.244339.00.22520.22580.22590.2257Table 8
Blow-up time for (4.1) with M=100 and k=1.
Blow-up time(Tb)rAPECE-implicit EulerMidpoint-implicit Eulerode45ode15s1215No blow-upNo blow-upNo blow-upNo blow-up2160.80410.83670.83670.84102170.68720.69810.69810.69562180.63100.63780.63780.63842190.59360.59870.59870.59902200.56550.56960.56970.56842370No blow-upNo blow-upNo blow-upNo blow-up3710.80540.86710.86710.87543720.69490.71260.71260.71423730.64270.65330.65330.65413740.60800.61570.61570.61623750.58190.58810.58810.5884


## 4.3. Discussion

We observe that as we increase the amplitude of the initial function,A, the blow-up time tends to occur earlier. We also observe that for smaller values of A, there is no blow-up.Considering the reaction-diffusion equation, and comparing the autonomous against the non-autonomous reaction term cases, we observe that the introduction of the non-autonomous term ensures that a much larger amplitude,A, in the initial function, is required for blow-up to occur. We also note that increasing k for fixed r or increasing r for fixed k increases the minimum amplitude for blow-up to occur.On the performance of the methods, we note a similar trend to what we observed in the ODE case. The adaptive implicit midpoint-implicit Euler method gives results that are significantly superior to the adaptive PECE-implicit Euler method and ode15s. In fact, its performance is comparable to ode45.


## 5. Reaction-Diffusion Equation: Two Space Dimensions

We solve the system (1.1) with the reaction term f=up(t,x,y), where p>1 with Ω=ℝ2. The system becomes
(5.1)ut(t,x,y)-uxx(t,x,y)-uyy(t,x,y)=up(t,x,y),t>0,x,y∈(0,1),(p>1),u(0,x,y)=u0(x,y)≥0,x,y∈(0,1),u(t,x,y)=0,t>0,x,y∈∂Ω.We use the method of lines (MOLs) to discretize (5.1) in space. For the spatial discretization, we choose uniform meshes for x and y, Dh:={xm:0=x0<x1<⋯<xM+1=1} and Ih:={yn:0=y0<y1<⋯<yN+1=1}, respectively, (with xm=mh and yn=nh) on Ω. We replace uxx(t,xm,yn) and uyy(t,xm,yn) (1≤m≤M,1≤n≤N) by the standard central difference approximation. We use the function Asin(πx)sin(πy) as the initial function with different values of A>0. We get the following system of ODEs for Um,n(t)≈u(t,xm,yn)(1≤m≤M,1≤n≤N):
(5.2)Um,n′(t)=Um+1,n(t)+Um-1,n(t)-4Um,n+Um,n+1(t)+Um,n-1(t)h2+Um,np(t),U0,n(t)=UM+1,n(t)=Um,0=Um,N+1=0,Um,n(0)=Asin⁡(πxm)sin⁡(πyn),1≤m≤M,1≤n≤N.We solve the system (5.2) with p=2 and h=(1-0)/M. Figure 7 shows the solution of (5.2), and Table 9 shows the blow-up times obtained using the four methods.Table 9
Blow-up time for (5.1) with M=10 and N=10.
Blow-up time(Tb)APECE-implicit EulerMidpoint-implicit Eulerode45ode15s21no blow-upno blow-upno blow-upno blow-up22no blow-upno blow-upno blow-upno blow-up230.47179710.47210220.47210280.4720943240.20341350.20369980.20370040.2036986250.15649470.15675470.15675530.1567541260.13065370.13089700.13089760.1308967270.11333130.11356190.11356250.1135617280.10059050.10081090.10081150.1008108290.09069000.09090170.09090230.0909017300.08270900.08291330.08291390.0829133Numerical solution of (5.1).
(a)
By PECE-implicit Euler(b)
By midpoint-implicit Euler(c)
By ode45(d)
By ode15s


### 5.1. Discussion

We observe similar results for the two space dimensions reaction-diffusion equation to the one space dimension case; that is, the adaptive implicit midpoint-implicit Euler method gives significantly better results than the adaptive PECE-implicit Euler method and ode15s. Its results are comparable with the higher order, computationally costly ode45.


## 5.1. Discussion

We observe similar results for the two space dimensions reaction-diffusion equation to the one space dimension case; that is, the adaptive implicit midpoint-implicit Euler method gives significantly better results than the adaptive PECE-implicit Euler method and ode15s. Its results are comparable with the higher order, computationally costly ode45.


## 6. Future Work

We would like to extend the blow-up computations to the Volterra integrodifferential equations, that is, cases where the reaction term is nonlocal.Given the relative simplicity of this new method, and cheaper computational expense, we conclude that it is much better than the higher-order RK-based solver, for implementation on problems of the nature studied in this paper. We would further like to test the new method on other problems in engineering and applied science, with the objective of proposing it for wider implementation.



---

*Source: 162034-2012-02-09.xml*