A posteriori error estimates for wave maps into spheres

We provide a posteriori error estimates in the energy norm for temporal semi-discretisations of wave maps into spheres that are based on the angular momentum formulation. Our analysis is based on novel weak–strong stability estimates which we combine with suitable reconstructions of the numerical solution. We present time-adaptive numerical simulations based on the a posteriori error estimators for solutions involving blow-up.

1. Introduction.This paper is concerned with the numerical approximation of wave maps, i.e., semi-linear wave equations with the point-wise constraint that the solution takes values in some given target manifold.They arise as critical points of a Lagrange functional for manifold valued functions and serve as model problems in general relativity [3] and in particle physics [1].We refer to [31,35,24] and the introduction in [30] for an overview on the general theory of wave maps.The monograph [19] contains a detailed introduction to the recent development in the analysis of wave maps.
There are two key challenges in numerically approximating wave maps: One is the point-wise constraint that make the function spaces in which solutions are sought non-linear and, the second, is gradient-blow up that leads to highly localized phenomena in space and time that need to be suitably resolved by numerical methods.A variety of different numerical methods that deal with the point-wise constraint using different approaches such as projections, penalties and Lagrange multipliers has been proposed [7,9,10,11,13,15,25] and for the methods of Bartels and coworkers a priori convergence analysis is available, in the sense that stability estimates are proven that imply convergence of subsequences to weak solutions.
It seems desirable to obtain more quantitative information on the accuracy of numerical approximations and, given the highly localized dynamics of gradient blowup, we aim to provide a posteriori error estimates.We will focus on a scheme whose a priori analysis was studied in [8] and [23].For the schemes at hand, convergence results are available even beyond gradient-blow up, i.e. limits of subsequences of numerical solutions are weak solutions, but quantitative estimates beyond singularity formation seem to be out of reach due to discontinuity of the solution operator [16], see Section 2 for details.Thus, we focus on estimates for the solution up to the blow up time.
Similarly to what was done in [21,26], we study errors entering via temporal discretization and, indeed, we restrict our study to semi-discretization in time.While the development of estimators for spatial discretization errors is certainly an important task in its own right, it is beyond the scope of this work.It will probably require a suitable extension of elliptic reconstruction techniques to harmonic maps, that does not seem available yet.Indeed, little has been proven concerning convergence of numerical schemes for harmonic maps.In certain situations, uniqueness and regularity of hamonic maps can be ensured, which allows to show high order convergence of so-called geodesic finite elements [22].In the general case, only existence of weak harmonic maps is guaranteed and for general triangulations and minimal regularity solutions, an a posteriori criterion is needed in order to guarantee weak convergence of numerical solutions [6].To the best of our knowledge no quantitative a posteriori error bounds are available for numerical approximation schemes for harmonic maps and harmonic map heat flows.
For a long time, a posteriori error control for (linear) wave equations was limited to first order schemes [12,20].Earlier works on adaptivity for wave equations can be found in [2,5,34].Quite recently, a posteriori error estimates for second order multi-step time discretisations of the linear wave equation were derived [21].
Due to the appearance of singularities, our goals are similar to those pursued in [14,26] that study blow-up for semi-linear parabolic equations.We also refer to these papers for earlier works on numerical approximation of blow-up solutions of nonlinear PDEs such as nonlinear Schrödinger equations or semi-linear parabolic problems.Let us mention that the blow-up mechanism in wave maps is rather different from the blow-up mechanisms in nonlinear Schrödinger equations or semi-linear parabolic problems where the L ∞ -norm blows up in finite time.
In order to derive the desired a posteriori error estimates, we use two main ingredients: Firstly, a suitable reconstruction of the numerical solution that can be understood as the exact solution of a perturbed version of the angular momentum formulation for wave maps into spheres and, secondly, a novel weak-strong stability principle.One key feature of our reconstruction is to ensure that the estimator is formally of optimal order, i.e., that it converges to zero with the same rate as the true error on equidistant meshes.
The remainder of this work is organized as follows: We review some facts from the analysis of wave maps in Section 2 and introduce the problem and basic notation in Section 3. Section 4 provides two stability estimates, one is based on a first order reformulation of the problem that is available when the target manifold is S 2 and the other Theorem covers the general case.In Section 5, we provide a posteriori error estimates for a numerical scheme that is based on the first order reformulation of the problem.The main contribution of this section is the construction of suitable reconstructions of the numerical solution; whereas computable bounds for the residuals, that appear when the reconstruction is inserted into the wave map problem, are postponed to the appendix.Finally, in Section 6, we report on numerical experiments using adaptive time stepping based on the a posteriori error estimators derived before.

Background on wave maps.
A specific feature of wave maps is that depending on the size of initial data (in suitable Sobolev norms) and the dimension of the target manifold either strong solutions may exist on arbitrarily long time intervals or solutions may exhibit gradient blow-up in finite time, see [30].But there are also larger classes of solutions, namely distributional or weak solutions, and in particular finite energy weak solutions.Note that weak solutions for example exist as accumulation points of subsequences of numerical schemes in [8].Existence of global weak solutions for finite energy initial data in 2 + 1 dimensions was established in [29].We will introduce our precise notion of finite energy weak solution in the next section.
If a (global) finite energy weak solution and a non-global strong solution exist, then uniqueness results as in [33,36] guarantee that both solutions agree until the appearance of the singularity.This is, under certain circumstances, strong solutions can be extended as weak solutions through the singularity.Conditions such as an energy inequality are needed to get this uniqueness result because, in general, weak solutions are not unique [36].
We provide weak-strong stability results Theorems 4.4 and 4.8 that can be seen as more quantitative versions of the (finite energy) weak-strong-uniqueness results in [36,33].Indeed, Theorems 4.4 and 4.8 assert that as long as there exists a strong solution for certain initial data, there are explicit bounds for the difference between this solution and solutions to perturbed problems even if those are only finite energy weak solutions.Note that this quantitative control requires the same regularity as the uniqueness result [33], see Remark 4.7 for more details.
We employ Theorem 4.4 (that assumes that the target manifold is S 2 ) in our a posteriori error analysis but since weak-strong stability results are interesting in their own right, we present a weak-strong stability result for general target manifolds in Theorem 4.8.While it can be thought of as an extension of Theorem 4.4 (where the target manifold is S 2 ) there are some significant differences on the technical level that will be discussed in Remark 4.9 and Remark 4.10.Those are the reason why we base our a posteriori error analysis on Theorem 4.4.This is discussed in more detail in Remark 4.11.
Uniqueness and stability properties of wave maps do not only depend on the target manifold but also on the dimension of the domain.Uniqueness for wave maps in 1 + 1 dimensions was shown in [37] and non-uniqueness in the supercritical dimension 3 + 1 was shown in [32,36].To the best of our knowledge, in the critical dimension 2 + 1 uniquenss of weak solutions to finite energy data is unknown and it is unclear whether imposing an energy inequality restores uniqueness in 2 or more space dimensions.An interesting observation in 2 + 1 dimensions is that solutions (even if they are unique) do not depend continuously on the initial data in the energy norm [16].It should be noted that our weak-strong stability results use the energy norm and are valid in arbitrarily many space dimensions.It follows from [16] that any such stability result, i.e. any bound for the difference of two solutions, measured in the energy norm, needs to involve a stronger norm (than the energy norm) of at least one of the solutions.This is a fundamental obstacle to deriving a posteriori error estimates (in the energy norm) that are convergent (i.e.go to 0 for τ, h 0, where τ, h denote spatial and temporal mesh width respectively) in case the exact solution does not have any additional regularity (beyond being a finite energy weak solution).
3. Problem statement and notation.For some bounded Lipschitz domain Ω ⊂ R m , some final time T > 0 and some n-dimensional submanifold without boundary N ⊂ R a wave-map is a map × Ω where T u N denotes the tangent space of N at u. Equation (3.1) needs to be complemented with initial and boundary data.To this end, maps u 0 : Ω → N and u 1 : Ω → R such that u 1 (x) ∈ T u0(x) N for all x ∈ Ω are fixed and one requires and homogeneous Neumann boundary conditions Strong solutions of the wave map equation satisfy an energy conservation principle and are critical points of the Lagrangian We reformulate the wave map equation (3.1) in order to see that this is a semi-linear wave equation.Let A p (•, •) : T p N × T p N → (T p N ) ⊥ be the second fundamental form of the compact submanifold N at a point p ∈ N .We denote the variables on [0, T)×Ω by (t, x) = (x α ), 0 ≤ α ≤ m.We raise and lower indices with the Minkowski metric (η αβ ) = diag(−1, 1, ..., 1) and we sum over repeated indices.Then a (strong) wave map is a map u = (u 1 , ..., u where A[u](Du, Du) stands for A i jk u ∂ α u j ∂ α u k 1≤i≤ , see [31].
A significant part of our analysis will consider the case that the target manifold N is the 2-sphere S 2 ⊂ R 3 and in this case (3.6) reduces to Let us also mention that, using angular momentum ω := ∂ t u × u, the wave map equation can be phrased as [23] (3.9) This variant is the one underlying the numerical scheme that we will study.

4.
A quantitative stability estimate.In this section, we establish a weakstrong stability result that complements weak-strong uniqueness results in [36,33] by providing bounds for differences between solutions.In particular, our estimates quantify the the impact of residuals that is crucial for the use of stability results in proving a posteriori error estimates.We prove weak-strong stability in the general case as well as in the special case of N = S 2 since we believe that the former nicely highlights the general geometric structure while the latter proof uses very elementary techniques and does not require any background in differential geometry.
for almost all 0 < t < T. We say that u satisfies a local energy condition provided for almost all 0 < s < t < T Remark 4.2.Note that the weak formulation in Definition 4.1 actually holds for all ψ ∈ W 1,1 (0, T; L 2 (Ω)) ∩ L 1 (0, T; L ∞ ∩ H 1 (Ω)) with ψ(T, •) = 0, due to a density argument where pass to the limit in A[u](Du, Du), ψ using majorised convergence, see the appendix of [18] for details.
We are going to compare a finite energy weak solution u to a strong solution (ũ, w) of the perturbed problem Lemma 4.3.For any two sufficiently regular functions ũ : Ω → S 2 , w : Ω → R 3 the following identities hold: (a) Theorem 4.4.Let u be a finite energy weak solution of (3.1)-(3.3)satisfying the local energy condition.Let p = m for m ≥ 3 and p ∈ (2, ∞] for m = 2. Consider (ũ, w), a solution of the perturbed problem (4.1).We ask for the following regularity: . Then, the difference at time t can be controlled via the difference in initial data and perturbation terms.Indeed, H defined by for almost all s < t ∈ (0, T) with where c q is the squared constant of the Sobolev embedding Remark 4.5 (Constants in the estimate).Explicit upper bounds for c q , for a variety of domains, can be found in [28].
Remark 4.6 (Regularity of weak solution).Note that if u is a finite energy weak solution, not necessarily satisfying the local energy condition, then (4.2) still holds in the special case s * = 0. Remark 4.7 (Regularity of strong solution).Note that our regularity assumptions with k = 1 correspond to the conditions in Struwe's result [33,Thm 2.2].The case studied in [33,Thm 2.2] corresponds to m = 2 and r u = r ω = r g = 0 in our notation so that Dũ =(ũ × w, ∇ũ).Thus, choosing k = 1, p = 2 + , we assume ∇ũ, w ∈ L 1 ((0, T); Struwe assumes Dũ ∈ X which also holds under our assumptions mainly because ũ ∈ L ∞ ((0, T) × Ω).There is the difference that we consider Neumann boundary data and Struwe has a solution on full R 2 .
Note furthermore that the conditions in our theorem allow us to use ũ × w as test function in (2) of Definition 4.
Proof.Let us fix some (arbitrary) s * < t * ∈ (0, T) and let us define for any 0 We will study lim ε 0 On the one hand, for every pair of Lebesgue On the other hand, we may decompose the integral at hand as Concerning E 1 ε , we observe that for every pair of Lebesgue points of t → E[u(t), ∂ t u(t)] Next, we consider E 2 ε for fixed ε and observe (4.8) We also note that (4.9) We insert (4.4) and (4.9) into (4.8) and obtain (4.10) where we have used integration by parts in the last equality.Equation (4.10) allows us to conclude Concerning E 3 ε , we use (4.4) and integration by parts to obtain (4.12) where we have used point-wise orthogonality of ũ to ũ× w in the last equality.Equation (4.12) allows us to conclude lim ε 0 Finally, we find (4.14)lim We combine (4.5), (4.7), (4.11), (4.13) and (4.14) to obtain This implies, for any pair p, q > 1 such that 1 p + 1 q = 1 2 , (4.16) where we have used that (4.17) and u, ũ ∈ S 2 .We insert (4.16) into (4.15) and obtain, for those q for which H 1 (Ω) embeds into L q (Ω) According to [17,Thm 21] equation (4.18) implies (4.2).Equation (4.3) follows by taking the square-root of (4.2) and induction in j.
We now come to the computations for a general closed target manifold N .Since for general N , we use one second order equation instead of a system of first order equations we have only one residual, but we split this into two parts, assuming that one part can be expressed as a time derivative.
We assume that there is a k ∈ [1, 4  3 ] such that We define the energy Then we have the inequality The constant c 3 is the diameter of N , c 3 = diam(N ), c 2 is the squared constant from the Sobolev embedding Proof.We pretend that u is a strong solution in order to keep the computations simpler.For a weak solution, we use the cutoff-trick with φ , that was introduced in the proof of Theorem 4.4.Let us denote Then we compute As A has normal values, we get that Using this orthogonality we get that We consider the term In order to estimate the term (A[u](Du, Du)−A[ũ](Dũ, Dũ), R 2 ) Ω we do the following: where we used the equation Using Hölder's inequality and generalized Hölder's inequality we get that for 2 p + 1 q = 1.For m ≥ 3, we use this inequality for p = 2m m−2 and get by Sobolev embedding where q = m 2 and c 2 is the squared constant from the Sobolev embedding.Thus, in the case m ≥ 3, we get for s = m the inequality For m = 2, we can choose 1 ≤ p < ∞ arbitrarily large in (4.22), which implies (4.24) holds for 2 < s arbitrary, in particular for s close to two.We put (4.21) and (4.24) together and get that Coming back to (4.20), we again use the formula ) for the terms with the second fundamental form there. Choosing any We define c 3 = c 3 (N ) = diam(N ) and c 2 to be the squared Sobolev embedding constant from H 1 → L s.We also use |∂ t ũ||Dũ| ≤ |Dũ| 2 and get that Together with where . From inequality (4.27), we get that and obtain, as the right hand side of (4.28) is monotone increasing in time Gronwall's inequality now implies (4.19).
It remains to explain that the conditions on ũ are strong enough for the above computations.We use ∂ t ũ as a test function in the weak formulation for u.This is allowed because the assumptions on ũ imply ∂ t ũ ∈ W 1,1 ((0, T); L 2 (Ω))∩L 1 ((0, T); L ∞ ∩ H 1 (Ω)).The conditions on ũ imply that all appearing terms are finite.Note particularly that and the latter is finite because k ≤ 4 3 implies k k−1 ≥ Remark 4.9.As one can see, the estimate for general N is slightly different than the one for spheres.The reason is that there is no such formulation as the angular momentum formulation (3.9) for general N .The formulation (3.9) has the advantage that it transforms the wave map equation into a system of equations that are only first order in t.Our numerical scheme is based on that formulation.Our reconstruction will not be W 2,1 but only W 1,∞ , a time derivative is a weak derivative.Furthermore, in the formulation (3.9) the term u × ω is the time derivative ∂ t u, which makes sense for u and ω being only continuous.This is also the reason that the two energies H and I differ in the time derivative.
Our motivation for splitting the residual in Theorem 4.8 is that if ũ in Theorem 4.4 happens to be twice weakly differentiable in time, then it satisfies where we we have used (4.4).This equation provides a connection between R 1 , R 2 and r u , r g and r ω .We formulated the estimates in Theorem 4.4 and Theorem 4.8 as similar as possible -in particular there is no norm of a time derivative of the residuals involved on the right hand sides of the inequalities -but in detail the estimates look a little bit different.
Remark 4.10.We note that our regularity assumptions for general N are only silghtly stronger than the one for spheres in Theorem 4.4.On the one hand, we need a second time derivative of ũ, see the remark above.On the other hand, we need the additional assumption Dũ ∈ L k k−1 ((0, T); L s (Ω)) where k k−1 ≥ 4 (and not k ≤ 2 as in the sphere case).The reason is the appearance of the term (4.19).Due to the different structure of the proof, this term did not appear in Theorem 4.4.
Note furthermore that by choosing k = 1 and considering the case R 1 = R 2 = 0 we have the same regularity assumptions as Struwe has in his result [33, Thm 2.2], see Remark 4.7.
Remark 4.11.In the next chapter, we derive an a posteriori error estimate based on the stability framework of Theorem 4.4 and not on Theorem 4.8.Since an angular momentum formulation is not available for general target manifolds, basing our analysis on Theorem 4.4, requires us to restrict our analysis to numerical methods for wave maps with values in S 2 .The reason why we use this, less general, stability analysis is its ability to handle reconstructions ũ, w that are in W 1,∞ in time.In contrast, a posteriori results based on 4.8 would require some reconstruction of the numerical solution that is in W 2,1 in time.Such a reconstruction was derived for the linear wave equation in [20] but it is not clear how to extend this construction to the wave map case.
5. Numerical scheme and a posteriori error estimates.We are concerned with a semi-discretization in time devised for wave maps into spheres in [23].The scheme is based on the reformulation (3.9) of the wave-map equation and reads (5.1) where d t is a backward difference quotient in time for a step size τ k and the fractional superscript denotes an average in time, i.e.
Note that (5.1) preserves the point-wise constraints |u k (x)| = 1 and u k ⊥ ω k for all k provided they are satisfied for the initial data.Indeed, this follows by multiplying (5.1) by u k+1/2 .A fully discrete version of this scheme, using a finite element discretization in space, was investigated in [8].There, it was shown that the scheme conserves energy, in the sense that and that in the limit of vanishing time step size subsequences of the numerical solution converge to weak solutions of the wave map problem.In contrast to the results of Bartels, we are interested in a posteriori error estimates, i.e., computable error bounds that can be evaluated once the numerical solution has been computed.While weak solutions can be defined beyond times at which a gradient blow-up has occurred there does not seem to be uniqueness beyond these times and, the results of [16] show that numerical schemes cannot be expected to converge with respect to the energy norm, once the exact solution has no additional regularity.As a consequence, we are only able to provide useful estimates up to the blow up time.Beyond apparent gradient blow-up in the solution the error bounds keep converging in τ but blow up for h → 0.
Our error analysis follows the approach outlined by Makridakis [27] in that it combines the 'energy type' stability result derived in Theorem 4.4 with a suitable reconstruction of the numerical solution.
The scheme (5.1) is very close to a Crank-Nicolson scheme and consequently we apply a reconstruction that is close to the reconstruction proposed in [4].A specific feature is that in order to employ the stability result from Theorem 4.4 we need to project the reconstruction into the target manifold.

Reconstruction.
In the sequel, we will define suitable reconstructions of the numerical solutions assuming that a sequence of numerical approximations at different points in time 0 = t 0 < t 1 < .... < t N is given: Firstly, we define preliminary, globally continuous, and piecewise linear interpolants using local Lagrange polynomials In addition, for any g ∈ C 0 ([0, t N ], L 2 (Ω, R 3 )) we define piecewise constant and piecewise linear interpolants by This allows us to rewrite the numerical scheme as Next, we define piecewise quadratic reconstructions via Note that u * , w are globally continuous in time since the trapezoidal formula is exact for linear functions and, in particular, u * (t n ) = u n and, thus, |u * (t n , x)| = |u n (x)| = 1 for all n and all x ∈ Ω.Finally, we define such that ũ is a map into S 2 and ũ(t n ) = u n for n = 0, . . ., N .
A straightforward computation gives (5.5) where a u | (tn,tn+1) := a n u and a ω | (tn,tn+1) := a n ω .5.2.Computable bounds for residuals.It should be noted that, due to the projection onto the sphere, ũ is no longer piecewise quadratic.Thus, it is not straightforward how to compute (norms of) the residuals from (5.5).In addition, the method at hand is formally second order, so that we should strive for a reconstruction making the a posteriori error estimator second order as well.It might seem obvious that ũ × w − I 1 (ũ × w) + a u and ∆ũ × ũ − I 1 [∆ũ × ũ] + a ω are second order in time, but this property does not seem obvious for Let us begin by decomposing the residuals into several parts: (5.6) r ω,2 = a ω .
Let us note that r u,2 , r ω,2 are both computable from the numerical solution, without computing any reconstruction, and are both going to converge to zero as τ 2 as long as the scheme is at least second order convergent.We will also give easily computable (though lengthy) bounds for the other parts of the residuals but we postpone their proof to the appendix.We will state the bounds for a representative time interval (t n , t n+1 ) and use the following abbreviations: x , B u xx to scale like τ n := t n+1 − t n as long as the exact solution is regular enough for the true error to be proportional to τ 2 n .In addition, we expect C ω , C u x , C ω x , C u xx to be bounded (uniformly in τ n ) as long as the exact solution is regular.Both expectations are confirmed by our numerical experiments.Lemma 5.1.Let us denote τ n := t n+1 − t n and let the time-step be chosen sufficiently small such that (A u ) 2 +τ n B u < 1 4 holds point-wise.Then, the residuals defined in (5.5) satisfy the following point-wise estimates.
An a posteriori estimate for the difference between the exact solution and the reconstructions (ũ, w) is easily obtained by combining Lemma 5.1 and Theorem 4.4.It should be noted that the definitions of the residuals and the way they enter into the error estimate is rather similar to [21,Theorem 3.1].However, we cannot avoid exponential-in-time growth of the error due to the non-linearity of the problem, even if the solution does not exhibit gradient blow-up.
6. Numerical experiments.In this section, we perform numerical experiments in order to demonstrate the scaling behaviour of the error estimator and to investigate time stepping strategies.We focus on the scaling (in τ ) of the error and of the error estimator.In order to obtain a fully practical scheme, we used a finite difference discretisation in space (together with the temporal discretisation discussed above), i.e. the scheme from [23].Whenever derivatives occur in the error estimator, they are approximated using finite differences.We consider the problem that was also studied in [7].These are smooth initial data for which a singularity forms in the numerical experiments.We believe that this reflects formation of a singularity in the exact solution, but, strictly speaking, we cannot be sure that this is what happens since the error estimator stops converging for h → 0 after a singularity seems to have formed.The fact that the error estimator does not converge once a singularity (seems to) have formed is related to the fact that our error estimator is an upper bound for the difference between the numerical solution and any finite energy weak solution in the sense of Definition 3.1 and these solutions are expected not to be unique once a singularity has formed.Note that this is not a contradiction to convergence of numerical solutions to some weak solution.In order to obtain an error estimator that converges after singularity formation one would need some technique to compare the numerical solutions to only one weak solution.It is currently unclear how to achieve this.
The problem data are given by N = S 2 , Ω = (− 1 2 , 1 2 ) 2 and (6.1) u(0, x) = The numerical solution develops a singularity at some time t > .2 to be precise.Thus, we expect the scheme to converge with order τ 2 in the energy norm Table 1: Error in energy norm and eoc up to time T = 0.2.We use a mesh constant h = 2 −9 and a reference time step τ = 2 −14 .Since the exact solution is unknown in this case, we use a numerical solution on a very fine mesh (h = 2 −9 , τ = 2 −14 ) as reference solution in order to approximate the error.We observe that the error indeed scales as τ 2 , see Table 1.
One of the goals we had in constructing the error estimator, i.e. the upper bound provided by combining Theorem 4.4 and Lemma 5.1, was that it is formally of optimal order, i.e. that it converges to zero with the same rate as the true error on equidistant meshes.This is guaranteed as soon as the quantities x are uniformly bounded in τ .This is indeed what we observe in the experiments we carried out.We do not report these numbers in detail but plot t 0 α(s) ds and t 0 δ(s) ds for several equidistant step-sizes τ and for different spatial mesh widths h, in Figure 2.Here α, δ are computable upper bounds for α and δ based on Lemma 5.1.These results show that α and δ are independent of the spatial mesh width h, as they are supposed to be since we are quantifying time discretisation errors.We also observe that α is proportional to τ 2 and δ is (essentially) independent of τ .This confirms our theoretical predictions, i.e., optimal convergence order of the error estimator in the smooth regime.
Fig. 2: Left: Temporal evolution of the integral of computable upper bounds for α based on Lemma 5.1.Note that, as long as the numerical solutions is "smooth", α is independent of the spatial mesh width h and proportional to τ 2 .Once numerical solutions show a singularity, α does depend on h.Right: Temporal evolution of the integral of computable upper bounds for δ based on Lemma 5.1.Note that the bounds for δ are independent of h and τ as long as the numerical solutions are "smooth", whereas δ is h dependent once the numerical solutions display singularities.It should be noted that δ contains finite difference approximations of ∇u L 2p and it is anticipated that if the gradient of the exact solution blows up the finite differences of the numerical solution approximating ∇u ∞ scale like 2/h.Due to the expected non-uniqueness of finite energy weak solutions we focus on time-stepping before and up-to singularity formation.Standard algorithms for mesh adaptation aim at equi-distributing 1 τj tj +τj tj α(s) ds amongst all time steps, i.e. the local step size is reduced if 1 τj tj +τj tj α(s) ds exceeds a given tolerance.As observed in [14], this strategy leads to excessive over-refinement close to blow-ups.Moreover, equation (4.3) shows that errors incurred in different time-steps are amplified by different factors.In order to obtain a scheme that is efficient up to blow-up, this needs to be taken into account and we follow a strategy that is similar to the one proposed in [14], i.e. from the j-th to the (j + 1)-th time step we increase the tolerance by a factor of exp( 1 2 tj tj−1 δ(s) ds).In Figure 3, we compare values of the error estimator and time step sizes for the two time-stepping strategies described above.It turns out that the scheme using the updated error tolerances uses much smaller step sizes in the beginning and larger step sizes at later times than the equidistribution strategy.We observe that at time t = 0.17 the strategy using updated error tolerances has used fewer time steps but leads to a smaller error bound.
We observe in Figure 3 that for time stepping with updated tolerances time step sizes stop increasing beyond a certain point in time.This is due to the fact that for time step sizes larger than some threshold value τ thresh (h) > 0 the fixed point iteration that we use for solving the non-linear system (5.1) in each time step, see [8, Algo.2] for details, does not converge.This results in the time step being automatically reduced independently of the error tolerance.
Appendix.This appendix is devoted to the proof of Lemma 5.1.Here, each quantity that is defined continuously in time is to be understood as its restriction to (t n , t n+1 ).Let us begin by giving explicit formulae for certain parts of r u and r ω .The We compare simulations with time step adaptation using a fixed tolerance 10 −4 to ones using updated tolerances starting at 10 −6 , so that at time t = 0.17 the number of computed time steps is comparable (approx.21000).
first such expressions make explicit the difference between some piecewise quadratic, globally continuous function and its piecewise linear interpolations: (6.2) Controlling | u|, |u * |.Let us now study how far away u, u * are from maps into the sphere: If u n+1 − u n is sufficiently small, a geometric argument implies and, therefore, Thus, the conditions in Lemma 5.1 imply Estimating r u,1 .Since ũ(t n ) = u(t n ) and w(t n ) = ω(t n ) we may rewrite r u,1 as (6.5) This proves (5.7).We obtain (5.8) by applying the product rule to (6.5) since | u(t, x)| ≤ 1.We also use the fact that ∇ũ is the projection of ∇u * onto the tangent space of the sphere, so that, provided |u n+1 − u n | < 1/2 holds, we have the point-wise estimate (6.7) |∇ũ| ≤ 2|∇u * |.
Estimating r u,3 .The key to estimating r u,3 is to control ∂ t u * • u * .We notice that (6.8) and, thus, (6.9) so that it remains to understand I 1 [ u × ω] • u.Using orthogonality, we obtain (6.10) We insert (6.10) into (6.8) and obtain (6.11) Moreover, due to (6.4) we arrive at (6.12) we have suitable bounds for all terms on the right hand side of (6.We recall |u n | = 1, which implies ∂ xj u n • u n = 0, for all n, so that (6.16) ).We insert (6.16) into (6.15) and obtain (6.17) Thus, we obtain (6.18) This completes providing bounds for the different components of r u and ∇r u .We insert (6.3) and (6.2) into (6.19) and obtain  Estimating r g .We note that orthogonality u n ⊥ ω n implies (6.25) Thus, using the definition of r g

Fig. 1 :
Fig. 1: Snapshots of numerical solution using h = 1/60 and τ = 2 −10 at different times.Arrows depict first two components of u h , color indicates direction of third component.

2 − 8 h = 2 − 6 Fig. 3 :
Fig.3: Time step sizes (left) and values of the combined error estimator from Theorem 4.4 and Lemma 5.1 (right) as functions of time for different (equidistant) spatial meshes.We compare simulations with time step adaptation using a fixed tolerance 10 −4 to ones using updated tolerances starting at 10 −6 , so that at time t = 0.17 the number of computed time steps is comparable (approx.21000).