B-Methods for the Numerical Solution of Evolution Problems with Blow-up Solutions Part II: Splitting Methods

B-methods are numerical methods which are especially tailored to solve non-linear partial di(cid:27)erential equation that have blow up solutions. We have presented in Part I a systematic construction of B-methods based on the variation of constants formula. Here, we use splitting methods as a second way to construct B-methods, and we prove several special properties of such methods. We illustrate our analysis with numerical experiments


Introduction
Nonlinear partial dierential equations (PDEs) arise in many important models in science and engineering, and very few of those models have closed form solutions.One therefore has to resort to numerical methods to compute approximations.If the partial dierential equation has further geometric properties, it is often an advantage for the numerical approximation to also have the same geometric property, which led to the research eld of geometric numerical integration.Much progress has been made over the last two decades in this area, see for example [28,40,27,15,7] and references therein.
In specic applications, the nonlinear PDE models can have solutions that blow up in nite time.This is in particular the case for combustion models [23,17,31,33], turbulent ow [38], nonlinear optics [36,41,42] and population dynamics [46,26].This blow-up indicates in general that the model is losing its validity, and it is therefore important to understand the precise behavior of the model when the blow-up time is approached.Studying such blow-up phenomena is necessarily done on a case-by-case basis, see for example [35,49,24,50,25,20,2,19,16], and the reviews [3,21].Blow-up can even happen when rst integrals are conserved, see [9] and references therein.The analysis of blow-up phenomena is an active eld of research, and many results have been obtained over the past two decades, see [8,43,37,44,29,13,14,47] and references therein.
The construction of numerical methods to approximate the blow-up time and rate of such models focuses in general on adaptive techniques.Very successful are moving mesh methods, see [10].In [11] self-similar solution techniques are employed for obtaining a scale invariant adaptive numerical method.A more direct numerical time stepping approach can be found in [30], where a numerical method using arclength ingredients is constructed and analyzed, and in [48], compactication of base spaces is combined with the validation of Lyapunov functions.Adaptive time stepping is also a very successful technique, where the time step is proportional to the inverse of the norm of the solution, see for example [45].
Considering blow-up as a geometric property is a more recent area of research, and so far mostly ad hoc constructions have been used to obtain numerical schemes with good blow-up properties, see for example [39].In a rst paper [6], we have shown how one can systematically construct B-methods using the technique of variation of the constant.The goal of this paper is to present a second systematic way of obtaining B-methods, using splitting techniques.

B-Methods Based on Splitting
To x ideas, we rst show the construction of splitting B-methods for a quasilinear parabolic problem.The construction for a few other scalar or systems of nonlinear partial dierential equations can be found in Section 4.

Model Problem and Assumptions
We consider the quasi-linear parabolic partial dierential equation u t = ∆u m + δF (u), for (x, t) ∈ Ω × (0, T ), u = 0, for (x, t) ∈ ∂Ω × (0, T ), u(x, 0) = u 0 (x), where δ is a positive constant, Ω is a bounded domain of R d and u 0 is a positive continuous function on Ω.In our analysis, we need the following Assumption 1.The function F is assumed to be positive, strictly increasing and strictly convex on (0, ∞), belonging to C 2 ([0, ∞)) and satisfying for b > 0. Then the function g(s) = ∞ s 1 F (σ) dσ is continuous and strictly decreasing on (0, ∞).The function G = g −1 is continuous and strictly decreasing on (0, M ), where M = lim s→0 g(s) ≤ ∞.Note also that g and G are positive with lim s→∞ g(s) = 0 and lim s→0 G(s) = ∞.
In order to be able to construct B-methods, we need to get an explicit form of g and often G. Examples of functions F which satisfy all these conditions are F (u) = e u , g(u) = e −u , G(u) = − ln u, F (u) = (u + α) p+1 , α ≥ 0, p > 0, g(u) = 2 − arctan(u), G(u) = cot(u).In problem (1) the nonlinearity in F is responsible for the nite-time blow-up and becomes increasingly important as we approach the blow-up time.The conditions imposed on F allow us to write explicitly the solution of the nonlinear ordinary dierential equation y t = δF (y).Indeed we get for any S > 0, where K is a constant, for all t satisfying K − δt ∈ (0, M ).It is then natural to seek integrators that exploit this information.In the following we present a new approach to obtain semi-discretizations in time for the semi-linear problem (1) from this exact solution.This approach allows us to derive many new Bmethods which are dierent from the B-methods obtained using the variation of constants approach in [6].

Construction of B-Methods based on Splitting
As suggested in Hairer, Lubich and Wanner [28] 1 , one way to exploit the exact solution of the nonlinear part of the equation is by using splitting methods.We illustrate this construction on the quasi-linear scalar PDE (1); the construction for a system of semi-linear PDEs is given in Section 4.3.
If we decompose u t = ∆u m +δF (u) into f [1] (u) = δF (u) and f [2] (u) = ∆u m , we can make good use of the fact that we know the exact ow φ [1] t of u t = δF (u) (note that φ t does not represent a time derivative) .Indeed, the exact ow of 1 It may happen that the dierential equation ẏ = f (y) can be split according to ẏ = f [1] (y) + f [2] (y), such that only the ow of, say, ẏ = f [1] (y) can be computed exactly.If f [1] (y) constitutes the dominant part of the vector eld, it is natural to search for integrators that exploit this information.
an equation y t = f (y) is the map dened by φ t (y 0 ) = y(t) if y(0) = y 0 , so in this case, using (3), we have φ [1] t (u n ) = G(g(u n ) − δt), for t < g(u n )/δ.Then we can choose any numerical integrator Φ [2] h for u t = ∆u m , and by composing the exact ow and the numerical integrator, we obtain two new methods for where h (see Section II.3 in [28]).Note that the two original methods Φ [2] h and Φ h (z 0 ) = z 0 + hf [2] t is the exact ow of u t = δF (u), so that its Taylor expansion is φ h (y 0 ) = y(h) = y 0 + hf [1] (y 0 ) + O(h 2 ).Therefore, the resulting methods Φ h and Φ * h are of rst order.This construction can only lead to methods of rst order, however as these two integrators are adjoint, we can use them as the basis of the composition method to construct methods of any desired order (see [28]).In particular, by choosing α 1 = β 1 = 1/2 for s = 1, we obtain a second-order symmetric method h ) is the forward (respectively backward) Euler method, the resulting method Ψ h corresponds to the midpoint (respectively trapezoidal) rule.
We saw that the exact ow of u t = δF (u) is given by φ , so we just have to choose a numerical integrator for the second part u t = ∆u m .For example, even though this problem is sti, we start with forward Euler Φ By composing these integrators with the exact ow φ t , we get two B-methods.The rst one is the Splitting Forward Euler B-Method (SpFE) which gives the explicit scheme and requires the condition g(u n + h∆u m n ) ∈ (0, M ).The second one is the which gives the implicit scheme and requires the condition g(u n ) − δh ∈ (0, M ).This scheme is studied in detail in Section 3 for the special semi-linear case, m = 1.

Instead of choosing Φ [2]
h to be forward Euler in (4), we could choose it to be backward Euler; then the adjoint Φ and the Splitting Backward Euler Adjoint B-Method (SpBE) Another possibility would be to choose Φ [2] h to be a second-order method, like the symmetric midpoint rule (SpMid) or the trapezoidal rule (SpTrap).However, the scheme becomes more complicated without necessarily bringing more accuracy, as the resulting scheme is only rst order.In order to get higher order methods, we need to compose rst order methods.The simplest way to obtain a second-order method is thus to construct where Φ [2] h and Φ [2] * h are adjoint rst-order methods.If we choose Φ [2] h to be forward Euler, we obtain the Second order Splitting Forward Euler B-Method (SoSpFE) and if Φ [2] h is chosen to be backward Euler, we get the Second order Splitting Backward Euler B-Method (SoSpBE) )) m .
(13) Similarly we can construct arbitrary high order splitting B-methods.

Truncation Error Analysis
In order to show that B-methods have the potential to be better than standard methods, we need to compare the local truncation errors of both types of methods.To start, we consider the problem u t = F (u) + Υ(u), where Υ can be a function or an operator (like the Laplacian in our example).We denote by φ the function that satises Keeping the notation introduced earlier, we have φ(t, v) = G(g(v) − t).We also consider the numerical method Φ applied to v t = Υ(v), with v(0) = v 0 .If v(t) solves this simplied problem, we have where E represents the local truncation error of the standard method.
We rst consider the B-methods obtained by applying the numerical method rst and use the result in the exact scheme (like SpFE or SpBE): starting with u 0 , we dene v 0 = u 0 and we apply the numerical method Φ to get v 1 = v(h) + E(h), then we set u 1 (h) := φ(h, v 1 ) = φ(h, v(h) + E(h)).To expand u 1 as a series of h, we need to compute its derivatives, , where the derivatives of φ are evaluated at (h, v(h) + E(h)).
From the denition of φ given in (14) (or using The values of E ′ (0) and E ′′ (0) depend on the standard method used, in particular for any consistent method, we have E ′ (0) = 0 and if the method is of second or higher order, we also have E ′′ (0) = 0.The Taylor expansion of the exact solution u is where the derivative Υ ′ (u 0 ) can be an operator, so the local truncation error of the B-methods is given by if a rst-order standard method is used, and for higher order standard methods we get To construct the adjoint methods, we rst use the exact scheme and then apply a numerical methods on the result.In other words, starting with the initial condition u 0 , we dene v 0 = φ(h, u 0 ), where φ satises condition (14), and we compute u 1 = Φ(h, v 0 ), where Φ is dened by (15) (to get a simpler notation, we denote the numerical method by Φ instead of Φ * ).The denition of Φ implies in particular that for all ξ, we have We now expand where we used the properties of Φ in ( 18) and ( 19) and the denition of φ given in (14).As the Taylor expansion of the exact solution u is given by ( 16), the local truncation errors of these B-methods are, as expected, for rst-order standard methods, and for higher-order standard methods we get We now need to show that in case of nite-time blow-up, the local truncation error of B-methods is smaller than that of the corresponding standard methods.We illustrate the dierence in truncation errors by considering the forward Euler method.
The local truncation error of the forward Euler method applied to the general equation y t = f (t, y) is given by which means that if we apply this method to u t = F (u) + Υ(u), we obtain On the other hand, if we apply forward Euler to Going back to (17) and ( 20) we obtain the truncation error of the corresponding B-methods, In order for the function F to be responsible for the nite-time blow-up, it needs to be superlinear at innity, while the remaining part Υ(u) becomes less important as u becomes large.We therefore expect the term F ′ (u 0 )F (u 0 ), which is present in τ s but absent in τ B and τ B * , to be large relative to the other terms.
As an example, let us rst consider the case where Υ(u) is a bounded function of u.We dene F (u) := e u and Υ(u) := sin(u).The local truncation error can then be written as for the standard method and for the specialized SpFE method.We see that the fastest growing term (e u0 ) 2 in τ s does not appear in τ B , while the other terms are of similar order.Given the size of this term compared to the remaining terms, τ B is considerably smaller than τ s .Going back to the case Υ(u) = ∆u, we observe numerically the same phenomenon.Indeed, with F (u) = 3e u and Υ(u) = ∆u, the local truncation errors are for the SpFE method, whereas for the (SpFE) * method, we have In this case, the term e 2u0 of τ s is also absent from τ B and τ B * , but it is not obvious that this term is much larger than the remaining terms.Some numerical experiments using Matlab show that the dierence between e 2u and the other terms is considerable and increases as u gets larger.Using the built-in adaptive method ode45 we computed the solution of u t = 3e u + ∆u on [−1, 1] with u 0 (x) = cos(πx/2), we then evaluated each of the four terms that appear in τ s .When t = 0.1660 (the blow-up occurs approximately at t=0.1664), the norm of the dierent terms are ∥∆(∆u 0 )∥ 2 = 342 439, ∥∆(3e u0 )∥ 2 = 1 466 377, ∥3e u0 (∆u 0 )∥ 2 = 1 542 768, and ∥(3e u0 ) 2 ∥ 2 = 16 544 121.So removing this last term from the local error greatly improves the results in this example.
3. Analysis of (SpFE) * for semi-linear parabolic problems We now analyze the properties of the (SpFE) * scheme applied to the model problem (1) for the special case of m = 1, i.e., when the problem is semi-linear.An explicit formula for the method is given by (8).By letting A := −∆, the scheme ( 8) for m = 1 can be written in the form

Existence and Uniqueness of the Solution
Since the scheme ( 23) is linear, it has a unique solution if and only if G(g(u n ) − δh) is well dened, i.e., g(u n ) ∈ (δh, M + δh).Since g is decreasing, M = lim s→0 g(s) and u n > 0, we have g(u n ) < M + δh, so the only condition is ∥u n ∥ ∞ < G(δh).We will need the following theorem due to Amann [1]: Theorem 1 (Amann).Let f ∈ C α ( Ω × R + ) be given, with α ∈ (0, 1), and assume that f (x, 0) ≥ 0. Then a necessary and sucient condition for the existence of a non-negative solution u ∈ C 2+α (Ω) of the BVP is the existence of a non-negative function Moreover, if this condition is satised, there exist a maximal non-negative solution û ≤ v and a minimal non-negative solution ū ≤ v in the sense that, for every non-negative solution u ≤ v of (24), the inequality ū ≤ u ≤ û holds.Theorem 2. If the function u n is positive in Ω, continuous in Ω, and satises then the scheme (23) has a maximal non-negative solution Remark 1.We can make the bound in condition (25) on the right-hand side as large as desired by choosing h small enough.This condition is necessary for the scheme (23) to be well-dened.
Proof.The constant C n is a supersolution of the scheme, if it satises 25) is satised and positive by denition of G (see Assumption 1), is a supersolution.Moreover, since f (x, 0) = 1 h G(g(u n ) − δh) > 0, we conclude using Theorem 1 from Amann.Since u n ≡ 0 is not a solution of the scheme, this result implies that there exists a non-zero nonnegative solution.Moreover the strong maximum principle applies (see for example [51]) and any nonnegative solution is positive on Ω. Uniqueness of the positive solution can also be obtained using the following result of Keller [34] with m = 0 and M = C n .
Theorem 3 (Keller).If there exist two constants m and M such that for all x ∈ Ω and all u 1 , u 2 such that m ≤ u 1 < u 2 ≤ M , we have f (x, u 1 ) ≥ f (x, u 2 ), then problem (24) has at most one solution u ∈ C 2 satisfying m ≤ u ≤ M .
Since f (x, u) dened in ( 23) is decreasing in u, we get the uniqueness of the solution, and we can show the same minimal time of existence as for the VBE scheme in [6]: Theorem 4. Scheme (23) has a unique positive solution u n for n such that the proof is exactly the same as the proof of [6,Theorem 3.11].
Finally, we recall that the scheme ( 23) is linear, so that no specialized nonlinear solver is required to solve for u n+1 .

Rate of Growth
We now prove some growth rate estimates for the scheme (23).Note that we will do this on a case-by-case basis for the functions F (u) listed in the introduction, since the estimate depends on the particular function at hand.We rst consider the function F (u) = e u , before turning our attention to the case of F (u) = (u + α) p+1 .Theorem 5. Let C 0 be a constant such that C0 , the function u n+1 given by satises for all Remark 2. Note that if Au 0 ≥ 0, we can take C 0 = δe ∥u0∥∞ , so that Proof.We prove this result by induction, using a supersolution approach.First, let us prove that if t 1 = h < T 2 , we have The function Since ln(1 k for x smaller than 1, we have Since 1 T2 = C 0 ≥ δe ∥u0∥ ≥ δe u0 , the bracket is negative and β is decreasing in h so inequality (29) holds for all h ∈ (0, T 2 ) if which is exactly condition (26), and we get (28).To complete the induction we assume that and we show that u n + ln T2−tn T2−tn+1 is a supersolution of (27), that is First, we note that since 1 T2 = C 0 ≥ δe ∥u0∥ , we have T 2 ≤ 1 δe ∥u 0 ∥ , and u 0 ≤ ∥u 0 ∥ ≤ ln( 1 δT2 ), and by induction By denition of u n , we have u n + hAu n = − ln(e −un−1 − δh), and from the induction hypothesis (30), we obtain − ln(e −un−1 − δh) > − ln[e −un ( T2−tn−1 T2−tn ) − δh], so that inequality (31) which simplies to (T 2 − t n )δ ≤ e −un , which is exactly (32).
Theorem 6.Let p > 0. Suppose there exists a constant C 0 that satises If t n+1 < T 2 := 1 C0 , the function u n+1 given by satises for all Proof.Throughout this proof, we will write w n = u n + α for all n.The recurrence can then be written as where we have Au n+1 = Aw n+1 , since A = −∆ annihilates the constant α.For the initial step, we want to show that v v 1 is therefore a supersolution, so by Theorem 1, there exists a solution w 1 of (35) that satises α ≤ w 1 ≤ ( T2 T2−h ) 1/p w 0 and w −p This completes the base case.
For the induction step, suppose the solution w n satises We now need to show that there exists a solution w n+1 such that We start by showing that v n+1 := ( T2−tn T2−tn+1 ) 1/p w n is a supersolution, i.e., we have It is clear from the denition that Substituting into (36) gives In other words, we need to show that ( T2−tn T2−tn+1 To prove the above inequality, we use the induction hypothesis: we have Therefore, the inequality (37) is true if which simplies to (T 2 − t n )pδ ≤ w −p n , which we know is true by the induction hypothesis.Thus, we have shown that v n+1 is a supersolution; by Theorem 1, a solution w n+1 exists and satises α ≤ w n+1 ≤ ( T2−tn T2−tn+1 ) 1/p w n , so that )pδ, which completes the induction step.

Numerical Blow-up
In this section, we want to prove that for values of δ large enough, the (SpFE) * method will blow up before a certain time T * ≤ ∞ 0 ds δF (s)−λs < ∞, with λ being the smallest positive eigenvalue of −∆.The existence of such a blow-up time in the continuous case has been shown by Kaplan in [32].Since we already proved that T * > T 1 = 1 δ (g(∥u 0 ∥ ∞ ), proving this result leads to exactly the same bounds as Kaplan for the discrete case.
To do so, we rst need to dene what we mean by numerical blow-up time.Suppose we use a numerical method of xed time step size h to integrate the model problem (1).We dene the numerical blow-up time T * h to be the smallest multiple of h such that the numerical solution ceases to exist.To estimate T * h , we adapt the approach used by Kaplan for the continuous problem to our semi-discretization: we show that there exists a nite time T * such that for all K > 0 and h small enough, there exists n < T * /h such that ∥u n ∥ ∞ > K, so that T * h ≤ T * for all h small enough. 2 We now state our main result.
Theorem 7. Suppose that δ satises where λ is the rst eigenvalue of −∆φ = λφ, φ = 0 on the boundary.We x some large positive constant K and choose ε ∈ (0, g(K)).Then there exists h * > 0 such that for all h < min(h * , g(K)−ε δ ), the numerical scheme has a numerical blow-up time T * ≤ ∞ 0 ds δF (s)−λs , in the sense that there exists Note that the proof presented in this section is constructive so that one can compute an explicit bound h * .We suppose thereafter that K and ε are xed.Remark 3. The assumption h < g(K)−ε δ implies that K < G(δh + ε) < G(δh) so that as long as ∥u n ∥ ∞ ≤ K, condition (25) is satised and scheme (39) has a unique positive solution.
2 While most of our previous results were following Le Roux's approach in [39], we could not use the same method as hers to prove this result.Indeed, a key element of Le Roux's approach is the use specic functionals and no equivalent functionals could be found for this scheme.

Outline of the Proof
We need to show that there exists n * < T * /h such that ∥u n * ∥ ∞ > K, where K is a xed large constant.Following the eigenfunction methods, we introduce the sequence (a n ), dened by where φ is the eigenfunction corresponding to the rst eigenvalue λ of −∆φ = λφ, φ = 0 on the boundary, with λ > 0, φ ≥ 0 and Ω φ dx = 1 (we can assume φ ≥ 0 since by Courant's theorem, the eigenfunction φ does not change sign in Ω).Our approach consists of nding n * such that a n * > K. Indeed we have We divide our proof into the following steps: 1.We prove that (a n ) is increasing.
2. We dene a(t), solution of a ′ (t) = δF (a(t)) − λa(t), a(0) = a * ∈ (0, a 0 ), which blows up in nite time at T = ∞ a * ds δF (s)−λs if δ satises condition (38).Dening D n = a n − a(nh), we need to bound D n from below in order to prove that for h small enough, D n is positive for all n for which a n and a(t n ) are well-dened.

Growth of the sequence (a n )
To prove that (a n ) is increasing, we need the following lemma.
Proof.Since ∥u n ∥ ∞ < G(δh), scheme (39) is well-dened.We multiply each side by φ and integrate over Ω to get Ω φ u n+1 − hφ∆u n+1 dx = Ω φ G(g(u n ) − δh)dx.Using the fact that u n and φ vanish on the boundary, the left-hand side can be rewritten as a n+1 − h Ω u n+1 ∆φ dx = (1 + hλ)a n+1 , and we obtain a n+1 = 1 1+hλ Ω φ G(g(u n ) − δh)dx.We now prove that the function f (x) := G(g(x) − δh) is convex for x ≥ 0. We have , and then which is positive since F being strictly convex implies that F ′ is increasing and we have G(g(x) − δh) ≥ x.Hence f is convex and we apply Jensen's inequality to get which completes the proof.

Denition of a(t) and D n
From now on, we assume that condition (38) is satised and h < g(K)−ε δ and ∥u n ∥ ∞ ≤ K.This implies that u n+1 is well-dened, thus so are a n+1 and D n+1 dened below.

Denition of a(t).
From Lemma 1, we have an+1−an h ), hence we will compare (a n ) with (a(t n )) where t n = nh and a(t) is the solution of where a * can be any xed number in [0, a 0 ).This limit simplies to By integrating this equation, we note that a(t) is dened on [0, T a * ), where δF (s)−λs ds < ∞, so that a(t) blows up at nite time T a * .Our goal is to show that a n is larger than a(t n ).
Denition of D n .For all n such that a n and a(t n ) are well-dened, we dene To prove Theorem 7, we will prove by induction that there exists h * such that ∀ h ≤ h * , ∀n such that ∥u n ∥ ∞ ≤ K, we have D n+1 > 0. The initial condition a * was chosen such that D 0 is positive, so assuming that D n is positive, we prove that D n+1 is also positive.First, we need to verify that a(t n+1 ) exists so that D n+1 is well-dened and t n+1 < T a * .Lemma 3. If D n > 0, the function a(t n +ξ), with ξ ∈ [0, h], is bounded above by a(t n + ξ) < G(ε), where ε is a xed number belonging to (0, g(K)) (see Theorem 7).
Proof.We introduce for t ≥ t n the function b(t), solution of This function can be written explicitly, b(t) = G(g(a(t n )) + δt n − δt), and we have a(t) ≤ b(t), ∀t ≥ t n .Moreover since δ satises condition (38), a(t) is increasing and we have a( and since a(t n ) < a n ≤ K and h < g(K)−ε δ , we get as required a( Hence D n+1 is well-dened and we rst bound it using Lemma 1, We then take a Taylor expansion of the right hand side function η(h) around h = 0, Thus, we have for some ξ ∈ (0, h), with ψ(x) = δF (x) − λx.Since η(h) is twice continuously dierentiable for all 0 ≤ h ≤ g(K)−ϵ δ , η ′′ (h) is continuous on the same interval, so there exists a (possibly negative) constant C 2 (which depends on δ, K and ϵ, but not on h) such that η ′′ (h) ≥ C 2 for all 0 < h < g(K)−ϵ δ .We are now able to prove Theorem 7.

Proof of Theorem 7
We suppose that ∥u n ∥ ∞ ≤ K and D n > 0. We now show that D n+1 > 0. Indeed, since a(t) blows up at time T a * with T a * ≤ T 0 = ∞ 0 ds δF (s)−λs , there exists ñ < T a * /h, such that a(t ñ) ≤ K and either t ñ+1 ≥ T a * or a(t ñ+1 ) > K.
The rst case implies that ∥u n ∥ ∞ > K for some n ≤ ñ, and in the second case, by the positivity of D n+1 , we have ∥u ñ+1 ∥ ∞ > a(t ñ+1 ) > K with t ñ+1 < T a * .Hence there exists n * < T 0 /h such that ∥u n * ∥ ∞ > K.
We assume that D n > 0 and we go back to (42) to write By induction, we obtain We assume that 1 + hψ ′ (a * ) > 0, so if ψ ′ (a * ) < 0, we need h to be smaller than If C 2 is positive, the positivity of D n+1 follows from (43).We now study the case C 2 < 0. We obtain dierent bounds on h depending on the sign of ψ ′ (a * ): , To prove that β(h) is strictly increasing for h > 0, we consider has exactly one solution h and if h < h we have D n+1 > 0.
The existence and uniqueness results of this section can be generalized to quasi-linear parabolic equations with power-like nonlinearities where Ω is a bounded domain in R d , m > 1 and α ≥ 0, see [5], but the upper bound blow up estimate remains currently open.

Numerical Results
We now test the new splitting B-methods on several non-linear partial dierential equations, and also compare them to the B-methods based on the variation of constants approach from [6] called VCFE, VCBE, VCMR and VCTR.

A Semi-Linear Parabolic Equation
For the rst example, we study the semi-linear parabolic equation u t = ∆u + δe u on the interval Ω = [−1, 1].We discretize the Laplacian operator in space using a fourth order nite dierence method with a ve point stencil and a mesh size of ∆x = 2/30.We set δ = 3 and u 0 (x) = cos(πx/2), which is concave on the whole interval.Using adaptive methods, we can estimate the blow-up time at T b ≈ 0.1664.We show in Table 1 the errors in the computed solutions up to T f = 0.1660 with dierent step sizes.We observe that the error of B-methods is approximately 10 times smaller for rst-order methods (and even more for SpBE and SpBEA) and 30 times smaller for second-order B-methods compared to standard methods.In Figure 1, we show these results graphically.As expected, the slopes of the lines corresponding to rst-order methods are approximately one, whereas the slopes of the lines corresponding to second-order methods are close to two.In Figure 2 we show the behavior in time of the methods as blow-up is approached, using h = 0.0001 and computing the solution up to T f = 0.1663.

A Quasi-Linear Parabolic Equation
Splitting B-methods can also be constructed for more general non-linear PDEs.We illustrate this now with the quasi-linear equation with power-type nonlinearities, u t = ∆u σ+1 + αu β+1 , with β > 0, σ > 0 and α ≥ 0. We split the PDE right hand side into f [1] (u) = αu β+1 and f [2] (u) = ∆u σ+1 , and using that the nonlinear part y t = αy β+1 is solved by y(t) = ( 1 K−αβt ) 1/β , the exact ow of the rst part is φ h to be the forward Euler method, so that Φ [2] * h is the backward Euler method we obtain the corresponding SpFE method and its adjoint, the (SpFE) * , If we choose Φ [2] h to be backward Euler we get SpBE, where v is solution of v − h∆(v σ+1 ) = u n , and its adjoint (SpBE) The second-order methods obtained by composing these methods are quite simple.If Φ [2] h is the forward Euler method, the composed method is SoSpFE where v is the solution of Similarly, the second-order method obtained using the backward Euler method for Φ [2] h is SoSpBE, given implicitly by where ].We show a numerical example for the quasi-linear equation u t = ∆u 2 + 8u 3 , on Ω = (−1, 1) with the same initial condition as above: u 0 (x) = cos(πx/2).The blow-up time is approximately T b ≈ 0.1128.We list in Table 2 the errors we obtained.We observe that the B-methods obtained by variation of the constant are more accurate than those obtained by splitting methods.Compared with standard methods, the errors are 10 times smaller for rst-order methods of the rst type and between 2 and 7 times smaller for rst-order methods of the second type.Among second-order methods, the method obtained by variation of the constant and the midpoint rule (VCMR) is remarkably better than the others, as its error is more than fty times smaller that the error of the standard midpoint rule.In Figure 3 we show the corresponding data graphically.The step-by-step errors are plotted in Figure 4 up to T f = 0.1110, when the solutions are computed using the timestep h = 0.0001.

A Semi-Linear System
In [18] and [19] , Friedman and Giga considered parabolic systems of the form u , where f and g are positive, increasing   and superlinear.They showed that the solutions exhibit a single-point blow-up.More complex systems of the form (u i ) t = ∆u i + f i (u 1 , . . ., u m ) were studied by Bebernes and Lacey [4], Gang and Sleeman [22] and Chen [12].In this subsection, we derive several specialized methods for the simple case We rst solve the nonlinear system of ordinary dierential equations y , where K and D are constants of integration, determined by the initial conditions, K = γe y(0) − δe z(0) and D = z(0) − y(0) − ln γ.Then for each choice of numerical integrator Φ [2] h applied to u t = ∆u, v t = ∆v, we obtain two schemes that are adjoint to each other.The forward Euler method leads to the explicit SpFE scheme We can also compose these methods to construct second-order splitting Bmethods.For these, we rst dene K n := γe un − δe vn , D n := v n − u n − ln γ as before.Then, if we choose Φ [2] h to be the forward Euler method, we dene w 1 and w 2 to be the solutions of If we choose to use the backward Euler method as Φ  3. and shown graphically in Figure 5.
In Figure 6, we show again the evolution of the solution as we approach blow up.We used h = 0.0001 and computed the solutions up to T f = 0.1170.Further examples can be found in [5, Appendix A].

Conclusions
We presented in this paper a systematic approach for deriving numerical integrators which are very accurate for semi-and quasi-linear parabolic and hyperbolic partial dierential equations exhibiting blow-up in nite time.We call this new class of geometric integration methods B-methods, where B stands for blow-up.Our construction is completely general, and can lead to B-methods for many other nonlinear partial dierential equations that were not considered in this paper.Because of their construction, which takes the blow-up behavior into account, all these methods will behave substantially better close to blow-up, while their behavior before blow-up is similar to classical time stepping schemes.

Figure 1 :
Figure 1: Error at T f = 0.1660 for rst-order methods (left) and second-order methods (right) applied to the semi-linear equation with δF (u) = 3e u , with dierent values of h.

Figure 2 :
Figure 2: Error for rst-order methods (left) and second-order methods (right) applied to the semi-linear equation with δF (u) = 3e u , for time steps close to T f = 0.1663.

Figure 3 :
Figure 3: Error at T f = 0.1000 for rst-order methods (left) and second-order methods (right) applied to the quasi-linear equation ut = ∆u 2 + 8u 3 , with dierent values of h.

Figure 4 :
Figure 4: Error for rst-order methods (left) and second-order methods (right) applied to the quasi-linear equation ut = ∆u 2 + 8u 3 , for timesteps close to T f = 0.1110.

Table 1 :
Error at T f = 0.1660 for rst-order methods (top) and second-order methods (bottom) applied to the semi-linear equation with δF (u) = 3e u .
u , with dierent values of h.

Table 2 :
Error at T f = 0.1000 for rst-order methods (top) and second-order methods (bottom) applied to the quasi-linear equation ut = ∆u 2 + 8u 3 .
3, with dierent values of h.