Numerical simulation of universal finite time behavior for parabolic IVP via geometric renormalization group

A numerical procedure based on the renormalization group (RG) is presented. This procedure will compute the spatial profile and blow up time for self-similar behavior. This will be generated by a family of parabolic IVPs, which includes the semilinear heat equation. This procedure also handles different diffusion structures, finite time extinction problems and exponential absorption with trivial modifications to the power law version. Convergence of the procedure is proved for the semilinear heat equation, which is a marginal perturbation to the heat equation, with a typical class of initial data. Numerical experiments show the accuracy of the method for various related problems. The main feature is simplicity: it will be shown that an explicit numerical method with a fixed mesh size is capable of computing very sharp approximations to these behaviors. This will include a priori detection of the logarithmic correction in the case of the semilinear heat equation.


1.
Introduction. The intent of this article is to establish the convergence for a numerical method that captures the universal behavior exhibited by the semilinear heat equation. This includes the logarithmic correction in the spatial contraction. The behavior occurs as the solution of the PDE approaches its blow up time: the finite time at which the solution becomes unbounded. The family of problems considered here will be given in Section 2 with blow up being independent of the size of the initial data. Further, the logarithmic correction is actually detected a priori by the method, which is applicable to other scenarios and not just the semilinear heat equation.
The same numerical approach with minor adjustments will also capture the universal behavior exhibited in other parabolic finite time scenarios, such as the quasilinear extension of the Frank-Kamenetskii equation or problems involving finite time extinction. The universal behavior involved in all of these problems is self-similar in nature.
The method is based on the renormalization group (RG) since RG has applications to the study of non-stochastic PDE through self-similar intermediate asymptotic behavior (IAB) of initial-boundary value problems (IBVPs). See Barenblatt (1996) [4] for more on IAB and its relation with self-similarity. The approach here uses an idealized version of the problem, simply dubbed the IVP, whose solution eventually becomes self-similar, possibly with logarithmic corrections. The presence of logarithmic corrections is known as approximate self-similarity.
In particular, the geometry of the domain is chosen to minimize the effects of the spatial boundary on the solution. The boundary conditions under consideration will be homogeneous Dirichlet conditions, which are currently idealized. Dimensional analysis, or scale transformation group invariance, allows the consideration of either limit: t → ∞ or t → T − , so that either global in time or finite time self-similar behavior can be investigated. This is because both quantities, t and T − t, have dimensions of time. The concern here will be for finite time behavior: t → T − .
Following the approach of Wilson (1971) [29] for statistical mechanics, geometric RG will be used to construct a sequence of maps on the space of IVPs, which is covered in Section 3. Convergence of the sequence can be shown to be equivalent to self-similar behavior in the solution for the original IVP, which occurs as an asymptotic behavior near blow up.
Having established that the RG approach yields the correct finite time behavior for the semilinear heat equation, the RG maps are modified in Section 4 to handle a discretization of the original IVP. The discretized maps are dubbed nRG, or numerical RG. In addition, a similar convergence result is obtained for the discretized version.
Once nRG's convergence has been established for the semilinear heat equation, numerical results are then given in Section 5. These experiments will support the robustness of the approach and accuracy of the implementation in the face of roundoff error. If nRG is applied without knowledge of the logarithmic correction, it will be shown that the implementation will alert an observer that a correction is present. Once the correction's presence is known, its details, like being logarithmic or not and any exponents, can be computed by a small modification to the method, Furtado (2004) [19].
Numerically, little is involved beyond a method which is sufficiently accurate away from the finite time behavior. This is chosen here to be forward Euler with standard centered spatial differences. The IVPs considered will be Cauchy problems in one spatial dimension with homogeneous Dirichlet boundary conditions, and experiments will restrict the Cauchy problem to closed spatial intervals.
There is also a perturbative version of RG for ODEs, see Goldenfeld et al. (1996) [15] and (1994) [16], along with Ei et al. (2000) [18], whose examples include PDEs as well. The ODE ideas have been examined mathematically, see Ziane (2000) [30] and Deville et al. (2008) [17] . There are also perturbative PDE approaches from Moise and Ziane (2001) [25], and Abou Salem (2010) [2]. Berger & Kohn (1988) [6] simulated the semilinear heat equation using essentially a geometric RG approach, but the problem was treated as a scale invariant problem. Hence, it did not address the logarithmic correction in the spatial rescaling. In addition, the iterative method used for the simulations was presented formally and required a nontrivial connection of multiple simulations. The results here will contain both a proof for the original spirit of Berger and Kohn's approach as well as improvements upon its implementation and approximating abilities. Moving mesh approaches, see Budd et al. (2009) [13] attempt to approximate the original (non-idealized) BVP by accounting for the self-similarity found in the IAB.
2. Parabolic blow up with power law source. The family of parabolic IVPs will be scalar Cauchy problems, which arise in areas including chemical reactor and combustion theory, quantum and fluid mechanics and population dynamics, see references in Levine (1990) with m ≥ 0, and p ≥ 1 + m when m > 0, but p > 1 when m = 0. The space of initial data B D,λ will be discussed in the next paragraph. It is known that these IVPs have a unique solution, u(x, t) with x ∈ S and that there exists a time T < ∞ such that lim sup t→T − sup x∈S u(x, t) = ∞ holds, while u is a classical solution for t < T . The solution is said to blow up in finite time in a classical sense and T is referred to as the blow up time. An x * ∈ S for which u(x, t) → +∞ as t → T − and x → x * holds is called a blow up point.
Single point blow up is the generic structure for the semilinear heat equation in one space dimension, see Herrero and Velásquez (1993) [22]. For convenience, only this blow up structure will be considered, and in the sequel blow up will refer to classical blow up in finite time at x * = 0. The blow up time for solution u may be denoted T u . Some convenient terminology: suppose S ⊆ R and f : In this paper, for all m ≥ 0, the focus will be on p < 3 + m. This prevents blow up from being initial data dependent, see Theorem (F) pg. 263 in Levine (1990) [23] which also contains a citation for the proof. The PDE is quasilinear when m > 0 and semilinear when m = 0. The proofs here will be focused on the case m = 0, due to the presence of the logarithmic correction. However, the general m case will be used in notation, since the behavior when m = 0 is a special case of the behavior when m > 0.
For convenience, fix m and p and denote L m f ≡ (f 1+m ) and A p f ≡ f p . Consider these as a basis for the operator space T = {(D, λ) ≡ DL m + λA p : D, λ ∈ R} with the subspace T + when D ≥ 0 and λ > 0. Given positive sequences, these generate a sequence whose range is in T + , denote by (D, λ) i ≡ (D i , λ i ) with convergence in T + defined by convergence of both D i and λ i . Given a sequence (D, λ) i , denote B i = B Di,λi . In addition, an IVP from (1) can be denoted as ((D, λ), f ) ∈ T + × B.
A local (around x = x * ) supersolution can be shown to both exist and be the limiting solution for vanishing diffusive strength. To this end, denote the Gaussian regular distribution by G(x, 1) ≡ (4π) −1/2 e −x 2 /4 and for q > 0, let . Suppressing x dependence, define G ≡ G(Dt) or G(D(t − s)), whose choice will be clear in context, and G is the fundamental solution for the heat equation. When variation of parameters is applied to (1) there is an integral representation, see [28], and using * to denote convolution wrt the spatial variable, and so u t (x, t) ≥ 0 for x ∈ R, 0 ≤ t < T and u(x, t), u t (x, t) are 0-max if f ∈ B.
Proof. Given an IVP ((D, λ), f ) ∈ T + × B, then a simple calculation shows that the is a supersolution to u(x, t) for some x-interval containing x = 0 (see e.g. Theorem 2, page 4 in Samarskii et al. (1995) [26]) and thus Given some D ≤ D, denote the solution of u t = (D, λ)u by u(x, t; D). It follows that u t = Du xx + λu p < Du xx + λu p for t ∈ (0, T ) and this implies u(x, t; D) is a subsolution to u(x, t; D), so u is monotonically increasing as D decreases for fixed x. Since u(0, t; D) is bounded above by v(0, t) for all D > 0, there exists u * (t) such that u(0, t; D n ) → u * (t) from below if D n → 0 as n → ∞ and u * (t) ≤ v(0, t). In addition, u * (t) is at least continuous over (0, T ), since a jump in u * (t) implies u t (0, t; D n ) → ∞, however u t (0, t; D n ) ≤ v t (0, t) < ∞ for t < T , hence u(0, t; D n ) → u * (t) uniformly by Dini's.

SIMULATION OF FINITE TIME BEHAVIOR VIA RG 3463
For semilinear problems, a logarithmic correction to the spatial rescaling arises, generating approximate self-similarity and a different group structure occurs. It has been shown in [5] and [22] that for m = 0, p > 1 with α = 1 p−1 and β = 1 which is analytic and solves βxΨ + αΨ − λΨ p = 0 with Ψ → 0 as |x| → ∞. In the next section, the result in (3) will be proved using RG maps. There is also non-logarithmic behavior for initial data with special structures, see [9]. For m > 0, the logarithmic correction disappears and β = p−1−m 2(p−1) , while the scaling profile is not known in closed form.
3. RG maps for parabolic finite time blow up. Geometric RG leverages group invariance of solutions to PDEs to rigorously establish self-similar behavior. The presence of self-similarity is usually mitigated by the spatial boundary, hence IAB provides time and space intervals over which the approach to self-similarity occurs. For a Cauchy problem though, the domain is invariant under scale transformations and so are homogeneous Dirichlet conditions. As a consequence, the approach to self-similarity persists in time and IAB becomes finite or global in time behavior, depending on the problem. An IBVP with homogeneous Dirichlet conditions whose domain approaches R under the action of the scale transformation group, will have its IAB given by the behavior of the Cauchy problem. However, the region of validity for the IAB will not be established here.
The important step in the current approach is that the limit with respect to time is replaced by iterations of a map based on the group, in this case, t → T − via iteration of the scale transformation group. If at least one stable fixed point exists, then this iterative approach can rapidly drive a system towards a stable fixed point. The fixed point would represent the finite time behavior for the system which is universal, since all points in the basin of attraction would achieve the same behavior.
In general, establishing convergence of an RG map would imply that the behavior, either finite time or global in time, must be self-similar and of the form given by the maps. Conversely, if the RG map fails to converge then either the parameters in the RG maps are incorrect, or the behavior is not self-similar at all, i.e. a different group should have been considered. It is noteworthy that the numerical version can be set up to calculate all aspects of the self-similar behavior including parameters, rather than use a priori values. So, if the numerical approach fails to converge, then this would imply behavior which is not self-similar. Note that being able to compute the various parts of the self-similarity would be extremely helpful in a case like an anomalous exponent, where the value is not known in advance.
For an IVP from (1) with m ≥ 0, denote (D, λ) = (D, λ) 0 ∈ T + and f = φ 0 ∈ B 0 , with solution u(x, t) = u 0 (x, t). Then construct the following sequences of initial data and strengths inductively for n ≥ 1, using some 0 < L < 1 where the choice of 0 < τ n−1 < T will discussed below. In addition, u n−1 (x, t) , handling the essential difference between the m = 0 and m > 0 cases. The exponents α and β are parameters for the RG maps which happen to be known for these problems.
The focus will now be on the m = 0 case. Note that properties of (1) based on the initial data will hold for each u i , i = 1, . . . , n and that Hence, if convergence occurs, then u n (0, t) → v Ψ (t) from below and T un → T vΨ = 1 from above. The iteration in (4) will first be shown to be well defined.
Given φ n−1 ∈ B n−1 then u n−1 (x, τ n−1 ) satisfies P1, P3 and simple calculations show that these are preserved by the rescaling Since Q β n L β x ∈ R and u t (x, t) is 0-max, it follows that (D, λ) n φ n (x) ≥ 0 and is 0-max as well. In turn, φ n ∈ B n so D n L : B n−1 → B n for each n ≥ 1.
Following [6] the stopping time τ n−1 can occur at a convenient height, which can never be zero starting from f (0) > 0 so that τ n−1 is well defined. Upon rescaling, φ n (0) = Ψ(0) ≡ Ψ 0 which avoids corruption due to large data. This choice is denoted τ BK . Another possible choice, denoted τ I , stops at the convenient time τ n−1 = (1 − L)T and allows the numerical implementation to compute α. Unfortunately, this approach does not converge due to a slowly accumulating error. However it provides an accurate computation of α prior to contamination, see Experiment 5 in Section 5 for more details. Using τ BK , the choice for the strength map S n L and uniqueness of solutions to (1) implies the solution u n of ((D, λ) n , φ n ) is which is an alternate way to represent the coupled sequences in (4). One can evaluate at t = 0 for D n L and substitute (5) into u t = (D, λ) n−1 u for S n L . The maps in (5) are called RG maps for parabolic finite time blow up.
By adopting (5), the original solution's information is passed to later iterations through an identity. For the quasilinear case, Q n = 1 and the choice (5) preserves the semigroup property so that property when m = 0, it is necessary to expand the group to two parameters for n ≥ 1. Hence, an important result has been achieved: the rescaled behavior at a time near blow up is replaced by an iteration of problems away from blow up. Let n * = n if m = 0, and n * = 1 if m > 0 along with t n = n 1 L i−1 τ i−1 . Applying (5) n times and letting t = 0 yields an important identity which can be used to show that t n → T Lemma 3.2. For m = 0 and τ = τ BK , the sequence t n generated by (4) converges to the blow up time from below and L −n (T − t n ) → 1 as n → ∞.
Proof. By construction, τ BK implies φ n (0) = Ψ 0 so that u(0, t n ) = L −nα Ψ 0 < ∞, and so t n < T and since t n is increasing, this implies t n → T * ≤ T . Evaluating (6) at Noting that the PDE is autonomous, denote by (τ n ) 0 the stopping time for initial data f = φ n , with (τ n ) i as the stopping time of the subsequent i th iterate, then Since D n → 0 as n → ∞ and φ n (0) is constant for all n, then Lemma 2.1 implies that |T un −T v φ | → 0. Further, since φ n (0) = Ψ 0 by construction, there is T v φ = T vΨ and so T un → 1 so that (7) implies L −n (T − t n ) → 1.
For the choice τ = τ I , then t n = n 1 L i−1 τ i−1 = T (1−L n ), and t n → T − as well. Note that if φ n converges then using (6) and Lemma 3.2, the behavior of u near blow up would be to leading order u ∼ ( and so limiting behavior will be written in terms of ξ, which does not reside on an interval which shrinks with time. where the limit in (ii) is uniform on compacts in ξ surrounding the origin.
(ii) For a fixed a > 0, consider E n ≡ max ξ∈[−a,a] |u n+1 − u n | (ξ| ln L| β , t). In addition, denote t n = t n + L n t and L β,n ln = (L n | ln L n |) β and G n ≡ G(Dt n ) or G(D(t n − s)), based on context, and F n (ξ, s) ≡ G n * u p L β,n ln ξ, s . Using the integral representation (2), there is and applying Hölder's inequality to the non-integral terms with the choice ||G n || 1 Using Hölder's on F n+1 as before and since [28] implies ||L αn u p (L β,n ln ξ − (·), s)|| ∞ < ∞, this shows the last integrand is finite, and so its integral goes to zero since t n+1 − t n → 0 + . Expanding F n+1 in a Taylor series around L β,n ln ξ, there is for some X ∈ [L β,n ln ξ, L β,n+1 ln ξ]. Expanding in the first integrand around t = t n and using G t = DG xx yields a new integrand of Hence, using Hölder's and [28] as before on both convolutions yields finite integrands, and L β,n ln along with t n+1 − t n force both integrals to zero. It follows that To verify that the limit of φ n is Ψ, i.e. to establish (3), let t = 0 which then implies φ n converges. Upon substitution of (5) into the PDE and simplifying, it follows that φ ∞ satisfies the ODE that Ψ does and the result follows.
If lim(T − t) α (I 2 ) iξ did not exist, then lim u iξ would not exist, which can only occur at (x * , T ). It will now be shown that this limit does indeed exist. The limit is one sided, so there is no possibility of a jump discontinuity. Since u t ≥ 0 at x = 0, there is no oscillating discontinuity in the solution u. Using Hölder's and [28], it follows that for any t ≤ T there is a constant (wrt time) C > 0 such that which is finite, and so the limit of (T − t) α (I 2 ) iξ , and hence e iξ , exists.
In addition, no oscillating discontinuity in time at any fixed ξ implies lim t→t * u (ξ, t) can be computed from lim k→∞ u(ξ, t k ) for any sequence t k → t * as k → ∞. Applying Severini-Egorov's shows e iξ is uniform in the limit on a set whose complement has an arbitrarily small measure. In this case, lim e iξ (ξ, t) = (lim e(ξ, t)) iξ = 0 on this set and hence u * iξ → Ψ i ae. Since G is analytic in the spatial variable at each t ∈ [0, T ], then u iξ is well defined and continuous. If u * iξ Ψ i for some ξ then continuity of (u * ) iξ forces u * (i+1)ξ → ∞ at the same ξ, a contradiction. Hence (u * ) iξ → Ψ i uniformly on compacts.
It follows from Lemma 3.3 that |u iξ | → ∞ at ξ = 0 and so, there exists t ξ < T so that t > t ξ implies max R×[0,t] |u iξ |(ξ, s) = |u iξ |(0, t) for i = 2, 4. Using the PDE, it can be shown that (3) may be differentiated wrt t as well as using λΨ p = αΨ + βξΨ and Lemma 3.3 for i = 2, which produces the logarithmic correction and eliminates the u xx term in the limit. There is a similar calculation for u tt . Hence, max R×[0,t] |u it |(x, s) = |u it |(0, t) for i = 1, 2 for some t t < T . Let t 3.3 = max{t ξ , t t }. 1. IAB is associated with fixed points of the map in (5). Self-similarity is invariant under the standard rescaling by L α and L β , and is associated with S i L = S L . For various fixed pt. approaches, see [10], [11], [8] and Galaktionov & Vazquez (2004) [20]. 4. Numerical RG procedure for parabolic blow up. For this entire section m = 0. It should be emphasized that for the upcoming numerical procedure, the number of grid points in space at any time step will be constant. In particular, only the locations that the grid points represent will be rescaled.
With numerical experiments in mind, a result of Velázquez (1993) [27] is important. Suppose the spatial domain in the Cauchy problem is changed to x ∈ (−a, a), for some a > 0, and that this interval properly contains the blow up point. Then the behavior of this new problem is the Cauchy problem behavior restricted to (−a, a). Properties of the initial data imply the solution will have one sided limits at x = ±a. So the domain when simulating (1) numerically will be [−a, a], whose choice will be motivated later in this section. In practice, one would like to use the smallest a possible.
The focus will now be on the forward Euler method for the short time problem. Given initial data in f j ∈ B * D,λ , let f j = U 0 j and discretize the temporal domain recursively via t 0 = 0 and t k+1 = t k + ∆t k where k ≥ 0 and with M k = min((U k 0 ) 1−p , 1) and B ≥ 3 max{D, 1}. In particular, using the given initial data's height, ∆t 0 can be computed, which in turn allows the approximation at t 1 to be computed. In addition, if the approximation has been computed up to step k, then ∆t k may be computed, which represents the next time step, in other words, the one needed to compute the next step k + 1. Formula (8) relies only on information computed at step k, namely U k 0 , making (8) an explicit calculation.
when |j| < N , and with extrapolation at the boundaries, fourth order accurate for the proofs and with local truncation error (LTE) 601 The computed U k j will approximate u(x, t) in the following sense Lemma 4.1. Using (9), (8) and (10) with a fixed 0 < ∆x < 1 and B ≥ 3 max{D, 1} andk > 0 such that tk < T , on a fixed ((D, λ), f ) ∈ T + × B with m = 0, and denoting E 0 = max j |U 0 j −u 0 j |, then the following holds, provided the global maximum of the error occurs only for |j| < N : there exists Ck and ck independent of ∆x such that for each 0 ≤ k ≤k and U k j blows up at Proof. Writing u k j − U k j = e k j and using (9), there is for 0 ≤ k ≤k − 1 where ν, µ ∈ [0, 1] and σ ∈ [−1, 1]. Concerning (10), it follows that the error can be written e k Let max j |e k j | ≡ E k over |j| < N and since V k < v k then Letting λp(v k ) p−1 ≡ m k , which is independent of ∆x, it follows that the global error obeys E k+1 ≤ (1 + m k ∆t k )E k + Mk 2,4 ∆x 2 ∆t k . This will be true provided the global maximum of the error occurs in the interior so that δ 2 j e k j < 0. Hence the choice of N should be large enough to properly contain the location of the maximum.
Using the discussion immediately after Lemma 3.3 allows Mk 2,4 to be written independent of x Concerning the boundary, and noting that the previous inequality for the global error can be used in the interior where the last inequality follows from v monotonically increasing in time.
The numerical domain size a should be chosen to keep the boundary away from the location of the error's global maximum and to minimize the contribution of u xxxx at the boundary. If so, then the boundary's LTE need not be introduced and a lower order extrapolation may be used in the implementation. A linear extrapolation has actually yielded satisfactory results in practice. In addition, the interior global maximum allows the discrete Laplacian δ 2 j e k j to be ignored in the estimate.
The aim is to now improve the obvious bound given by λpv p−1 . This is too coarse because it is slightly larger than the reciprocal of the rescaling for the nRG procedure, and would produce a diverging error bound. Noting that I T ≡ T 0 v p−1 (s)ds = ∞, this fact can be leveraged to improve the bound to be contractive under the nRG map.
To facilitate proving the contraction, a new comparison between the true and numerical solution will be introduced. Consider u and U at the respective times when both achieve the same height at j = 0, i.e. given U k 0 define s k so that u(0, s k ) = U k 0 for all 0 ≤ k ≤k, with ∆s k = s k+1 − s k and note s 0 = 0. Denote u k j = u(j∆x, s k ) and e k j = u k j − U k j . Also introduce ρ k = ∆s k ∆t k and note that Lemma 4.1 implies that ρ k → 1 as ∆x → 0 + for each k.
Further, it can be established that for small enough D then ρ k = 1 − λp 2B ∆x 2 + O(∆x 4 ). To see this, note that u(0, s k ) → v(s k ) and U k 0 → V k as D → 0 + , and so there exists small enough D so that there is for any k. Writing u(0, s k ) = U k as v(s k ) = V k + O(∆x 4 ) and solving yields Using this in ∆s k = s k+1 − s k and after factoring and expanding in a Maclaurin series in ∆x, the result follows. Note that the series expansion alternates. It can be shown that the new error e k j has a bound similar to e k j , noting that Proof. Modifying (12) to account for the time evaluation difference, there is Hence, letting E k = max j |e k j | and Remark 2. Note that in the proof of Lemma 4.2, j m = 0 since by definition e k 0 = 0 for each k. In addition, u(x, t) ∞ for a given x = 0 and I T = ∞, and it follows that there exists t 1 large enough so that t 0 λpu p−1 (x, s)ds < I t when t > t 1 . Then, for a given ∆x and recalling t 3.3 from the discussion following Lemma 3.3, there is a k 4.2 so that min{s k , t k } ≥ max{t 3.3 , t 1 } and c k < exp(I s k ) when k > k 4.2 .
A numerical procedure based upon RG maps, dubbed nRG, is designed to compute aspects of self-similar behavior a priori. Computations of Ψ and T * are natural based on (5), and so nRG using τ BK will be a discrete version of (5). Discretizing (1), let (D, λ) * 0 = (D, λ) * and (φ * 0 ) j = f j , with solution (U 0 ) k j = U k j . It is best to pick the final point in time to coincide with the approximation U having height L −α Ψ 0 . To this end, let k i−1 − 1 be the largest time level for which ( for each i = 1, . . . , n. Note that the spatial rescaling is handled by rescaling the gridpoint locations, and projecting back onto an integer grid via a polynomial interpolation to keep ∆x fixed. Although this does introduce an additional error, it is controllable and keeps the method simple to implement. It should be emphasized that the number of grid points in space at any time step is constant, only the locations they represent are being rescaled. The analog of Lemma 3.1 holds for (9) and (8), see [1]: a discrete maximum principle holds and f j satisfies P1 , P2 , ⇒ U k j satisfies P1 , P2 , which also holds for the choice (10). The rescaling preserves P1 , P2 as well, and so D i L : B * i−1 → B * i . In addition, since δ k and δ 2 j transform like u t and u xx respectively, ∆t k . Preserve the semigroup property by having U i ≈ u i so that is computed, which is dubbed the solution's mass. If tk 0 λpu p−1 (j m ∆x, s)ds < I sk could be established in advance, e.g. specific initial data for which the detailed structure of u(j m ∆x, s) can be estimated, then there is no restriction on L for convergence of nRG. However, for this condition to hold in more generality, one must be imposed upon L. Theorem 4.3. Fixing ((D, λ), f ) ∈ T + ×B, and > 0, so that for each N ∈ Z + and U k j generated via (9), (8) and (10), there is for small enough ∆x > 0 an 0 < L * < 1, so that when 0 < L < L * , there is Proof. The error at iterate n is decomposed as follows where Proposition 1 implies h n → 0 + as n → ∞. Note that the second term in the first maximum is the same as the first term in the second maximum, only written with different notation. The evolution's LTE from Lemma 4.2 is now denoted by M kn 1,2,4,ζ . Using a quadratic interpolation for the spatial rescaling produces a second order LTE, which can be bound by C|u xx |k 0 ∆x 2 ≡ M kn 2 ∆x 2 for some C > 0. Since T un → 1, then max n T un ≡ T M exists and t kn < T M < ∞. Let n 4.2 be such that D n < D 4.2 but D n−1 ≥ D 4.2 , which exists since D n → 0 as n → ∞. Apply Lemma 4.2 to g n with n > n 4.2 ,k = k n−1 and E 0 = g n−1 , and upon including the spatial rescaling process, there is Note that M kn 1,2,4,ζ and M kn 2 depend on the partials (u n ) t , (u n ) tt , (u n ) ξξ , (u n ) 4ξ , along with (u n ) p−1 , which all converge as n → ∞ upon using Proposition 1, Lemma 3.3 and the discussion thereafter. Hence, the maximum over n for each of these functions exists, so denote M 1,2,4,ζ ≡ max n M kn 1,2,4,ζ and M 2 ≡ max n M kn 2 . Re-applying Lemma 4.2 and denoting P * i = L iα i j=1 c kn−j with an empty product being unity, there is The constraint on L now arises from Remark 2. For τ * n to be large enough so that k n ≥ k 4.2 then L must be small enough so that the rescaling height L −α Ψ 0 is large. Shortening s kn to s n , it follows that ln(c kn ) where κ = k n4.2 . Choose ∆x so g ∞ < and with h n → 0 + as n → ∞, (15) holds.

5.
Results of numerical experiments. An implementation of the nRG procedure was enacted to show that accurate results will occur for reasonable choices of ∆x, n, a and L. For all numerical experiments, a mesh size of ∆x = 0.01 was used, and for most N was chosen to be 5000, so that the spatial domain was [−50, 50]. The nRG parameters were L = 0.25 and n chosen to achieve at least a 1% relative error from known benchmarks. Using D = 1, λ = 1 and only linear extrapolation for the boundary conditions, as well as linear interpolation for the spatial rescaling, the following experiments were treated numerically: 1. Semilinear: m = 0 and p = 1.05, 2, 3, 4 with f (x) = e − x 2 2 . The m = 0, p = 1.05 case required n = 1000, and note that in this case α = 20.
These semilinear cases will be used to verify the convergence results of the previous section. Initial data is in B. 2. Critical β = 0: m = 1 and p = 2, 3 and f (x) = max{1 − 0.3x 2 , 0}.
A gradient diffusion (see page 3478) with exponential source, which arises in turbulent diffusion, is a quasilinear extension of the Frank-Kamenetskii equation. It is considered with small m, where asymptotic results should be valid. 7. Finite time extinction: m = 0, p = 0.5 and f (x) = max{1 − 0.3x 2 , 0}.
When p < 1, the sign of α changes, implying finite time blow up becomes finite time vanishing, i.e. extinction. There are no modifications to make to the blow up case for this nRG implementation. For the quasilinear case, the profiles are not known in closed form and will be denoted by Φ(x). An interesting case arises when p = 1 + m so that β = 0, and assuming compactly supported initial data. The profile has bounded support for all t ≤ T with Θ(x) = 0 for |x| > L * 2 and L * ≡ 2π(m + 1)m −1 , the fundamental length for the given diffusion strength, see Chapter 4, Section 1 of [26].
Heuristically, starting with compact support and knowing that the solution will approach no spatial rescaling (β = 0), the behavior would be expected to have compact support as well. The ODE governing Θ can be integrated to yield a periodic function with a period given by the fundamental length. So behavior with compact support can be achieved by defining Θ to be one period of the ODE's solution, centered around the plow up point, and then zero outside that period.
• Experiment 1-2 The results computed via nRG are summarized in the table below. Displayed in the first row are the relative convergence errors for φ * n , with grid dependence j suppressed The second row has relative accuracy errors, defined for m = 0 and m = 1, p = 2 as For m = 1, p = 3, the correct profile is not known in closed form, and so the accuracy errors are then defined by applying a discretized version of (18), (see below) denoted M, to the computed profile, and (17) is replaced by The last two rows of the table contain T * n from nRG along with T * d from DNI. The convergence errors for Experiment 1 are much bigger than for Experiment 2, due to the presence of the logarithmic term. However for all three experiments, they are small enough to indicate possible convergence of the implementation. The histories of the convergence errors in these cases show strong evidence that the computed sequences have stabilized. In addition, the accuracy errors indicate that the implementation is indeed arriving at the correct blow up profile. Several arbitrarily chosen experiments were run for n = 5000 with no change in behavior after stabilizing.    The asymptotic results for m > 0 are due to Budd et al. (1998) [12]. These are based on formal calculations, but have been established numerically by independent means, and will be used as benchmarks for small m. Denoting and as t → T − , there is u(x, t)(T − t) α → Φ(ξ) uniformly on compact subsets of ξ. If 0 < m 1, there is an 'inner' behavior where C i q = mα 2p for m − 1 2 y 1, and an 'outer' behavior where C o q = m 4pλ for y 1, • Experiment 4 Results for the a priori detection of the logarithmic correction are given in Figures 3 and 4 below, see Li and Qi (2003) [24] concerning logarithmic corrections. In Figure 3, the semilinear's mass history is plotted using the proper Q i = i i−1 , i ≥ 2 versus the improper Q i = 1. The lack of convergence in the mass history implies a lack of convergence for naïve nRG, indicating to the implementer that a correction is present.
For the quasilinear case, now Q i = 1 is proper, and both M n and ∆φ ∞ 's histories are plotted. The exponent on the number of iterations n for ∆φ ∞ is approximately −1 from n = 40 onward, casting doubt on convergence of the nRG maps when Q i = i i−1 .
• Experiment 5 Displayed in Figure 5 below are the results from computing the blow up rate, which uses the exponent computing choice of τ I = (1 − L)T . The RG maps now will have a similar structure R i L u i = L αi u i−1 (Q β i L β x, τ I + Lt)  for i ≥ 2. Iterating the maps, letting t = 0 and denoting A n = L nα− n 1 αi , there is φ n (x) = A −1 n L nα u(n β L nβ x, (1 − L n )T ) when the semigroup property is preserved. Setting x = 0 in this identity, it follows that 1 = φ n (0) = A −1 n L αn u(0, t n ). Denoting t n = (1 − L n )T there is L n = 1 − tn T = T −tn T . In addition, there is t n = (1 − L n )T → T since L n → 0, and so A n → Ψ 0 T −α by the results in Section 3, which implies the RG maps converge and α i → α.
The nRG map becomes n β * L nβ j with α * i computed by resetting the height of the data back to 1 and the following scaling factor is produced by iterating this map. To be an a priori computation, the value of α is replaced by α * n in the numerical calculation, dubbed A * * n . In particular, since α * i α i as i → ∞ when ∆x > 0, then A * n = L n(α−α * n ) A * * n . This may be small for small n if the method is accurate, but this factor will increase as n increases. Hence a stopping rule for the calculation of A * n is in order. When |α * n − α n | ≤ tolerance, then A * n will no longer be computed at each iteration, i.e. if α * has given evidence of converging numerically, then A * will be considered to have numerically converged as well.
To highlight this issue, the nRG maps here are implemented without the stopping rule for A * . The instability is seen in the mass and L ∞ histories after 500 iterations and the method fails to approximate the profile as well as the τ BK version. Starting at top left and going clockwise, histories of α * n , the scaling factor A * n and ln ∆φ ∞ n vs. ln n are given, along with a plot of the history of M n . The latter two plots clearly show the error accumulation. Note that the α * computation is not affected by this situation.  Figure 6 below contains the results for gradient diffusion equation u t = (|u x | m u x ) x + e u with m = 0.1, references can be found in [12]. Like most quasilinear cases, the spatial profile is only known asymptotically, and its height at x = 0 is m 2 u(x, t) → − log(T − t) +Φ(x(T − t) − 1 m+2 ) Following this lead, define the RG maps with 0 < L < 1 as The numerical method used is still explicit; a backward difference is used for u x , and then a forward difference is used for the final derivative (·). Starting at the top right and working clockwise, the first two smaller subplots contain a history of T * n and ln ∆φ ∞ n vs. ln n (16), where the dotted line is a line of slope -1. The next two smaller subplots contain a history of M n (14) along with a spatial plot of φ * n superimposed with the parabolic asymptotic result that is valid near the origin. These four plots surround a large spatial plot of φ * n . In the center plot, the logarithmic asymptotic result valid for |x| → ∞ is used as the benchmark. If λ < 0 instead of λ > 0, then M → 0 + , which is achieved by α's reversal in sign, so that t * → 0 implying finite time behavior. Hence, α and λ's sign reversal implies growth becomes decay, or extinction in finite time. The nRG maps for blow up need no modification to handle the extinction case. Results are given in Figure 7 for a semilinear extinction problem, with p = 1 2 . This uses the same arrangement as gradient diffusion, except the upper left hand spatial plot is replaced by a comparison plot since the benchmark is known.