Numerische Mathematik

Numer. Math. (2010) 115:395–432 DOI 10.1007/s00211-009-0282-y

Finite element approximations of harmonic map heat flows and wave maps into spheres of nonconstant radii L’ubomír Banas ˇ · Andreas Prohl · Reiner Schätzle

Received: 13 June 2008 / Revised: 2 October 2009 / Published online: 20 January 2010 © Springer-Verlag 2010

Abstract We prove the existence of weak solutions to the harmonic map heat flow, and wave maps into spheres of nonconstant radii. Weak solutions are constructed as proper limits of iterates from a fully practical scheme based on lowest order conforming finite elements, where discrete Lagrange multipliers are employed to exactly meet the sphere constraint at mesh-points. Computational studies are included to motivate interesting dynamics in two and three spatial dimensions. Mathematics Subject Classification (2000)

65M12 · 65M60 · 35K55 · 35Q35

1 Introduction Let  ⊂ Rd (for d ≥ 1) be a bounded domain, and Sm−1 ⊂ Rm (for m ≥ 2) be the unit sphere. The energy of a map w :  → Sm−1 is defined as  1 |∇w|2 dx. (1.1) E(w) = 2 

Critical points are called weakly harmonic maps into the sphere [29], where also the related L 2 -gradient flow is motivated from topological grounds (‘homotopy problem’): L’. Baˇnas Department of Mathematics and the Maxwell Institute for Mathematical Sciences, Heriot-Watt University, Edinburgh EH14 4AS, UK e-mail: [email protected] A. Prohl (B) · R. Schätzle Mathematisches Institut der Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany e-mail: [email protected] R. Schätzle e-mail: [email protected]

123

396

L’. Baˇnas et al.

Find u : T → Rm , and a Lagrange multiplier λ H : T → R, such that ut − u = λ H u, and |u| = 1 ∂u =0 ∂n u(0, ·) = u0

in T := (0, T ) × ,

(1.2)

on ∂T := (0, T ) × ∂,

(1.3)

in .

(1.4)

Below, we refer to equations (1.2)1 , (1.3), (1.4), together with an algebraic equation for |u| as Problem A; the present case where (1.2)2 holds is prototypical for extended models in micromagnetics [4,5,8,23], liquid crystal theory [1], or color image denoising [9,31,32,34]. A different evolutionary problem in this setting is the wave map into the sphere Sm−1 , with applications in general relativity [13,14,27], utt − u = λW u and |u| = 1 ∂u =0 ∂n ∂t u(0, ·) = u0 u(0, ·) = u0 ,

in T ,

(1.5)

on ∂T ,

(1.6)

in .

(1.7)

In the following, we refer to (1.5)1 , (1.6), (1.7), with an algebraic equation for |u| as Problem B; again, (1.5)2 is a possible special constraint in this more general setup of the problem. For the above situation where |u| = 1 in T , we have λ H = |∇u|2 , and λW = |∇u|2 − |ut |2 , and stationary solutions to (1.2)–(1.4), and (1.5)–(1.7) are harmonic maps to Sm−1 . The goal in this work is to study mappings u : T → Rm , where the target is a smoothly varying sphere in space–time, i.e., |u(t, x)| = M(t, x) > 0 for all ( t, x ) ∈ T , and a given M : T → R+ \ {0}. Our study is motivated by the recent work [4] on ferromagnetic modelling, where shrinking, extension, and conservation of the modulus of the magnetization depend on locally applied temperature fields. So far both analytical and numerical studies for Problems A and B have focussed on fixed target manifolds: For Problem A, existence of weak solutions for targets Sm−1 , and possible (finite time, finite energy) singularities, i.e., lim supt→T − u(t, ·)W 1,∞ = ∞ are shown in [15,16,18,29,30]. Weak solutions, and their finite-time blow-up behavior is studied in [12,15,16,18,29,30]. Hence, it is of interest to study evolutions from Problems A and B where the target (smoothly) changes in space and time, and analytical difficulties hence exceed those from above. We refer to Sect. 4, where computational studies show different qualitative behaviors of solutions. Numerically, it is not easy to realize the constraint |u| = M > 0 in a fully discrete scheme, due to restricted flexibility of standard finite element functions, and restricted regularity properties of weak solutions. In fact, penalization strategies to approximate the constraint converge, but the resolution of the solution, e.g., close to singularities, crucially depends on chosen finite scales of penalization and discretization parameters. Projection schemes are available, where iterates are projected back to the sphere at every node z ∈ Nh , but iterates do not satisfy a (discrete) energy law, and convergence is only known in the case of strong solutions; cf. [23]. It is only recently, that different

123

Finite element approximations of harmonic map heat flows and wave maps

397

strategies are developed to handle the case M ≡ 1, in such a way that the constraint holds for involved finite-element functions at all nodes z ∈ Nh of the triangulation, that iterates satisfy a (discrete) energy law, and convergence of the scheme is known for general data. One approach is based on the following equivalent formulation of (1.2), ∂t u ∧ u − div (∇u ∧ u) = 0

in T ,

(1.8)

where a ∧ b = (a i b j − a j bi ), for 1 ≤ i < j ≤ m, and a, b ∈ Rm . Likewise the higher products are defined a1 ∧ · · · ∧ ak ∈ k Rm for a1 , . . . , ak ∈ Rm , whose linear combinations are called k−vectors. More precisely, for an orthonormal basis {e1 , . . . , em } of Rm the set {ei1 ∧ · · · ∧ eik | 1 ≤ i 1 < · · · < i k ≤ m} forms an orthonormal basis for k Rm . On using reduced spatial integration for element-wise affine, globally continuous finite element functions, and the trapezoidal rule for integration in time, as well as a lumped discrete Laplacian, an implicit discretization of (1.8) in space and time is constructed in [7,8], whose solutions unconditionally (sub-)converge to weak solutions of Problem A for vanishing discretization parameters. Rather than enforcing the sphere constraint at all points in T in (1.2), we use approximate discrete Lagrange multipliers to validate |Un | = 1 for iterates 1 ≤ n ≤ N at all mesh-points z ∈ Nh of the triangulation Th of . Here, ‘approximate’ refers to the fact that a proper (time-)discretization for the term on the right-hand side of (1.2)1 is needed to validate a discrete energy identity for corresponding iterates. Again, element-wise affine, globally continuous finite element functions are used. On employing reduced spatial integration, and trapezoidal rule, an implicit discretization of (1.2) in space and time is proposed in [12], where iterates satisfy a discrete energy identity, and (sub-)converge to weak solutions of Problem A for vanishing discretization parameters, provided that k ≤ Ch 2 . This mesh constraint is sufficient for solvability for finite time and space discretization parameters k, h > 0, whereas neither the discrete energy law nor the sphere constraint require mesh constraints. Correspondingly, a convergent discretization for Problem B is obtained in [12]. The situation is more complicated for the target being a sphere with nonconstant radius given by M : T → R+ \ {0}, where the related Lagrange multiplier in (1.2)1 to guarantee |u| = M in T is λH =

1 2 [∂t

− ]M 2 + |∇u|2 M2

in T .

(1.9)

To numerically solve Problem A for |u| = M in T , the numerical approach based on the reformulation (1.8) of (1.2) that crucially benefits from |u| = 1 in T , is not

123

398

L’. Baˇnas et al.

promising. Instead, we propose a discretization strategy that involves approximate discrete Lagrange multipliers to ensure for iterates that |Un (z)| = M(tn , z) at nodes z ∈ Nh , for all 1 ≤ n ≤ N . A major problem to overcome is then to properly account for the modified right-hand side of (1.2)1 where (1.9) is used in the time-discrete setting in particular, such that iterates are of length M : T → R+ \ {0} at mesh-points z ∈ Nh , and to obtain a formula for related approximate discrete Lagrange multipliers {λn+1 H }n≥0 ; this problem is solved by the nontrivial construction (2.9) in Scheme A, where temporal increments of M are properly reconciled with those of the computed iterates, and contain algebraic combinations thereof. We may then employ an inductive argument from [12] to conclude solvability for the resulting scheme, provided a meshconstraint is valid; moreover, corresponding iterates satisfy a discrete energy estimate. This property, in turn, then allows to construct weak solutions of the heat flow harmonic map to spheres with nonconstant radii, as proper limits of existing convergent sequences of solutions of Scheme A for vanishing discretization parameters, provided that certain regularity assumptions for the given positive function M : T → R+ \{0} are met. A similar program is then possible for Problem B with |u| = M in T , for a given smooth function M : T → R+ \ {0}. In this case, the related Lagrange multiplier is λW =

1 2 2 [∂t

− ]M 2 + |∇u|2 − |ut |2 M2

in T .

(1.10)

If compared to Problem A, or [12] where the sphere of constant radius is considered, the main difficulty is to construct a scheme, where increments of M and those of the modulus of the solution are properly adjusted, so that a discrete energy estimate may be obtained. However, (combinations of) increments are now related to second time-derivatives, and the construction (3.3) of the approximate discrete Lagrange multipliers {λn+1 W }n≥0 in Scheme B involves several coupling terms of increments of M, and those of the solution. A key result for the convergence analysis showing convergence of iterates from Scheme B to a weak solution of Problem B is the stability result stated in Lemma 3.1. The remainder is organized as follows: In Sect. 2, we discuss solvability and stability properties (discrete energy inequality) of Scheme A, and construct weak solutions of Problem A as proper limits (for ( k, h ) → 0) of iterates of Scheme A. Section 3 correspondingly shows properties of solutions of Scheme B, and their convergence to weak solutions of Problem B. Computational experiments to motivate interesting dynamics for both problems are reported in Sect. 4.

2 Harmonic map heat flow into a sphere of nonconstant radius A special case of Problem A is M ≡ 1, which is the heat flow of harmonic maps to the unit sphere. In the first part of this section we clarify the notion of weak solutions of Problem A where the target is nonconstant in space–time, and possible singular behaviors of them.

123

Finite element approximations of harmonic map heat flows and wave maps

399

2.1 Singularities, blow-up of classical solutions, and weak solutions In this section, we study possible singularities and a blow-up procedure for classical solutions of Problem A. First, for d ≥ 1, we multiply (1.2) by ut , integrate by parts and use (1.9),    1 d |ut |2 dx + E(u) = λ H u, ut dx = λ H ∂t (M 2 ) dx dt 2       ≤ C M 1 + |∇u|2 dx ≤ C M (|| + E(u)). (2.1) 

Integrating in time and using Gronwall’s lemma, we get  esssupt∈[0,T ] E (u(t, ·)) +

|ut |2 dxdt

T

≤ exp (C M T ) E(u0 ) + (exp (C M T ) − 1) || =: E,

(2.2)

and conclude that the energy remains bounded in finite time. Next, we localize the energy and write for x ∈  and ρ > 0, 

1 E ρ (w, x) = 2

|∇w(y)|2 dy. Bρ(x)∩

  We choose a cut-off function ϕ ∈ C0∞ B2ρ (x) with 0 ≤ ϕ ≤ 1, such that ϕ ≡ 1 on Bρ (x), and |∇ϕ| ≤ Cρ . Again we multiply (1.2) by ut ϕ 2 , integrate by parts and obtain 



1 |∇u|2 ϕ 2 dx 2     2 = λ H u, ut ϕ dx − ∇u, ut 2ϕ∇ϕ dx |ut |2 ϕ 2 dx +

d dt







C M (1 + |∇u|2 )ϕ 2 dx +

≤ 

1 2



|ut |2 ϕ 2 dx + Cρ −2



 |∇u|2 dx. 

Absorbing and integrating in time over [t, t + τ ], we get by Gronwall’s lemma and recalling (2.2) 1 E ρ (u(t + τ, ·), x) + 2

t+τ |ut |2 dxdt t



  ≤ exp (C M τ ) E 2ρ (u(t, ·), x) + (exp(C M τ ) − 1) C M + Cρ −2 E .

123

400

L’. Baˇnas et al.

Now for any ε > 0 there exist τ = τ (E, M, ε) > 0 such that for 0 < ρ ≤ 1, E 2ρ (u(t, ·), x) ≤ ε

⇒

E ρ (u(s, ·), x) ≤ 2ε

for all t ≤ s ≤ t + τρ 2 .

(2.3)

Now we restrict ourselves to two dimensions  ⊂ R2 . There the energy E is scale invariant. Then using the multiplicative Sobolev inequality 1/2

ϕ L 4 (B1 ) ≤ Cϕ L 2 (B

1)

  1/2 1/2 ∇ϕ L 2 (B ) + ϕ L 2 (B ) 1

1

∀ ϕ ∈ W 1,2 (B1 ),

and similarly at the boundary, we can estimate solutions u : T → R3 of Problem A satisfying E 2ρ (u(t, ·), x) ≤ ε for some 0 < ρ ≤ 1 and appropriate ε > 0, under certain smoothness requirements  2 k for M and ∂ in the space C Bρ (x × (t + τ ρ2 , t + τρ 2 ) for all k ≥ 0, see [30]. Lemma 2.1 Let  ⊂ R2 , and u : T → R3 be a weak solution of Problem A with smooth M : T → R+ \ {0}. Then a point (t, x) ∈ (0, T ] × , for T < ∞ is a singularity if and only if lim sup E ρ (u(s, ·), x) > ε for all ρ > 0, st

and there are at most finitely many singularities in finite time. Moreover, we can set up a blow-up procedure by rescaling   um (t, y) := u tm + ρm2 t, xm + ρm y

for 0 ≤ t ≤

τ , 2

for appropriate τ = τ (M, , T, ) > 0, with E(u 0 ) ≤ , and appropriate tm  t, xm → x, such that um (t, ·) → u¯ smoothly in compact subsets of R2 and similarly at the boundary to a non-constant harmonic map u¯ : R2 → M(t, x) · Sm−1 , at the boundary after reflexion, with finite energy and which therefore extends smoothly to S2 . This study regarding occurrence of singularities justifies to look for weak solutions for Problem A.   Definition 2.1 Let M ∈ C 2 T ; R+ \ {0} be given, and u0 ∈ W 1,2 (, Rm ), such  that |u0 | = M(0, ·) almost everywhere in . We call a function u ∈ L ∞ 0, T ; W 1,2  (, Rm ) ∩ W 1,2 (T , Rm ) a weak solution of Problem A, if for all T > 0 there hold

123

Finite element approximations of harmonic map heat flows and wave maps

401

(i) u(0, ·) = u0 in the sense of traces, (ii) |u| = M almost everywhere in T , and (iii) for all φ ∈ C ∞ (T , Rm ) there holds   φ dxdt = 0, ∂t u ∧ u φ dxdt + ∇u ∧ u ∇φ T

(2.4)

T

in short ∂t u ∧ u − div(∇u ∧ u) = 0. As will turn out below, this definition of a weak solution for Problem A is more suitable than the following, where requirement (iii) in the above definition is replaced by (iii’): for all ψ ∈ C ∞ (T , Rm ) there holds 

  ψ − ∂t u, ψ + ∇u, ∇ψ

1 2 [∂t

T

− ]M 2 + |∇u|2 u, ψ M2



dxdt = 0, (2.5)

in short ∂t u − u =

1 2 [∂t

− ]M 2 + |∇u|2 u. M2

We have the following result. Lemma 2.2 Let M : T → R+ \ {0} and u : T → Rm be given, such that all requirements of Definition 2.1 apart from item (iii) are fulfilled. Then u is a weak solution to Problem A if and only if item (iii’) is valid. φ , we obtain (2.4). Proof Wedging (2.5) by ψ = uφ Conversely, we assume that (2.4) holds. We use the Hodge operator ∗ : k Rm → m−k Rm , which is a linear map that acts for an oriented orthonormal basis {ei1 , . . . , eik, e j1 , . . . , e jm−k } of Rm by ∗(ei1 ∧ · · · ∧ eik ) = e j1 ∧ · · · ∧ e jm−k . We observe for ψ ∈ L 2 (0, T, W 1,2 (, Rm )) ∩ L ∞ (T , Rm ) that ψ = ∇u ∧ ∇ ∗ ψ . ∗ ∇u, ∇ψ Extending um+1 ≡ 0, and ψ m+1 ≡ 1, we may assume that u, ψ are linearly independent everywhere. If moreover u ⊥ ψ , we see for ψ ∧ u)/M 2 ∈ L 2 (0, T, W 1,2 (, n−2 Rm )) ∩ L ∞ (T , n−2 Rm ) ω := ∗(ψ that

123

402

L’. Baˇnas et al.

ψ = u ∧ ω. ∗ψ Wedging (2.4) by ω, which is justified by approximation, we calculate  0=

 ∂t u ∧ u ∧ ω dxdt +

T

T



=

∇u ∧ u ∧ ∇ω dxdt 

∂t u ∧ u ∧ ω dxdt + T





∇u ∧ ∇(u ∧ ω) dxdt

T

ψ dxdt + ∂t u ∧ ∗ψ

= T

∇u ∧ ∇ ∗ ψ dxdt

T





∗ ∂t u, ψ dxdt +

= T

ψ dxdt, ∗ ∇u, ∇ψ

T

hence 

 ∂t u, ψ dxdt +

T

ψ dxdt = 0 for ψ ⊥ u. ∇u, ∇ψ

(2.6)

T

For any ψ ∈ C ∞ (T , Rm ), we write ψ = αu + ψ ⊥ with ψ , u /M 2 ∈ L 2 (0, T, W 1,2 (, Rm )) ∩ L ∞ (T , Rm ) α := ψ and ψ ∈ L 2 (0, T, W 1,2 (, Rm )) ∩ L ∞ (T , Rm ), ψ ⊥ ⊥ u . By (2.6), we get 

 ∂t u, ψ dxdt +

T

T

 =

ψ dxdt ∇u, ∇ψ 

∂t u, αu dxdt + T

 

= T

 

= T

T

 1 1 ∂t |u|2 + |∇u|2 α dxdt + ∇|u|2 · ∇α dxdt 2 2

T

1 [∂t −]M 2 + |∇u|2 α dxdt = 2

which is (2.5).

123

∇u, ∇(αu) dxdt

 1

T

2 [∂t

− ]M 2 + |∇u|2 u, ψ dxdt, M2  

Finite element approximations of harmonic map heat flows and wave maps

403

Let us consider time-independent data M(t, x) = M(x) . In this case, we get from (2.1) that  d |ut |2 dx + E(u) = 0, dt 

and (2.2) turns into  esssupt∈[0,T ] E (u(t, ·)) +

|ut |2 dxdt = E(u0 ).

T

More generally, let us consider a weak solution of Problem A satisfying the energy inequality  (2.7) esssupt∈[0,T ] E (u(t, ·)) + |ut |2 dxdt ≤ E(u0 ). T

For a global-in-time solution u : [0, ∞) ×  → Rm , we choose tk → ∞ and put uk (t, ·) := u(tk + t, ·) . We see for any 0 < T < ∞ that uk  L 2 (0,T,W 1,2 ()) , ∂t uk  L 2 (0,T,L 2 ()) ≤ C(E(u0 ), T, M) < ∞, hence for a subsequence uk uk ∇uk ∂t uk

→ u∞ → u∞ → ∇u∞ → ∂t u∞

  strongly in L 2 0, T, L 2 (, Rm ) , pointwise almost  everywhere in   T, weakly in L 2 0, T, L 2 (, Rdm ) , weakly in L 2 0, T, L 2 (, Rm ) .

As |uk | ≤ sup M is bounded, we see

∂t uk ∧ uk , → ∂t u∞ ∧ u∞ , weakly in L 2 (0, T, L 2 (, Rm )), ∇uk ∧ uk → ∇u∞ ∧ u∞ , hence passing to the limit in (iii’) ∂t u∞ ∧ u∞ − div(∇u∞ ∧ u∞ ) = 0, and u∞ is a weak solution of Problem A. Moreover, by the energy inequality (2.7) ∞ 

∞  |∂t uk (t, x)| dx dt =

|∂t u(t, x)|2 dx dt → 0,

2

0 

tk 

hence ∂t u∞ = 0, and u∞ is a time-independent, hence stationary solution.

123

404

L’. Baˇnas et al.

We have proved the following lemma. Lemma 2.3 Let M :  → R+ \ {0} be time-independent. For a weak global solution u : [0, ∞) ×  → Rm of Problem A satisfying for all T > 0 the energy inequality  esssupt∈[0,T ] E (u(t, ·)) +

|ut |2 dxdt ≤ E(u0 ),

T

every sequence tk → ∞ admits a subsequence such that {u(tk , ·)} converges weakly in W 1,2 (, Rm ) to a stationary solution. The goal in the following section is to construct weak solutions of Problem A by suitable limits of finite element functions, which solve Scheme A given below. This result, in turn, justifies computational experiments in Sect. 4.2 to motivate possible finite time finite energy blow-up behavior of solutions in the case M : T → R \ {0}. We need some notation to propose Scheme A: Let Vh ⊂ W 1,2 () be the lowest order conforming finite element space subordinate to a triangulation Th of , and Vh = [Vh ]m . By Nh , we denote the set of all nodes associated with Th , and {ϕz ; z ∈ basis for Vh ; we define the nodal interpolation operator Ih : C() → Vh Nh } is nodal by Ih w := ∈Nh w(z)ϕz , for all w ∈ Vh . Below, (·, ·)h denotes the discrete version (reduced integration) of the L 2 -inner product on Vh , i.e.,  (v, w)h :=

Ih (vw) dx = 



βz v(z)w(z)

∀ v, w ∈ C(),

z∈Nh

 where βz =  ϕz dx, for all z ∈ Nh . We use dt ψ n := k −1 {ψ n − ψ n−1 }, and ψ n+1/2 := 21 {ψ n+1 + ψ n } with n ≥ 1, for a sequence {ψ n }n≥0 , and for an equidistant time-step size k > 0. For a given continuous function ξ : T → R+ \ {0}, and a n sequence {ψ n }n≥0 , we further denote ψξn := ξ(tψn ,·) , as well as dt ψξn := k −1 {ψξn − n+1/2 := 21 {ψξn+1 + ψξn }, with n ≥ 1. Let M n := M(tn , ·). By C, C˜ > 0 ψξn−1 }, and ψξ we denote bounded, generic constants, dependent on  and T , where the second additionally depends on M.

2.2 Scheme A: construction of weak solutions of Problem A We propose a practical numerical scheme, whose iterates construct weak solutions of (1.2)–(1.4), where |u| = M in T , for given M : T → R+ \ {0}. A main task here is to provide a formula for the approximate discrete Lagrange multipliers n+1 } {λn+1 n≥0 satisfy a discrete energy identity. We remark H }n≥0 such that iterates {U that, in contrast to the continuous Lagrange multiplier in (1.9), the computation of each n+1/2 at neighboring λn+1 H at a single node z ∈ Nh requires to consider values of U ones—which is due to finite parameters k, h > 0 in the scheme.

123

Finite element approximations of harmonic map heat flows and wave maps

405

Scheme A. Let U0 ∈ Vh , with |U0M (z)| = 1 for all z ∈ Nh . For n = 0, 1, 2, .., find (Un+1 , λn+1 H ) ∈ Vh × Vh such that for all  ∈ Vh , and all z ∈ Nh there holds       n+1 n+1/2  = λn+1 M U ,  , (2.8) dt Un+1 ,  + ∇Un+1/2 , ∇ H M h h ⎧ n+1/2 ⎪ if U M (z) = 0, ⎨0     n+1 n+1 n+1/2 n+1/2 dt M (z) n n+1/2 λ H (z) = βz M n (z) U (z),U M (z) + ∇U ,U M (z)⊗∇ϕz ⎪ ⎩ else, n+1/2 βz M n+1 (z)|U M

(z)|2

(2.9) n+1/2

where U M

:=

1 2

  n , for Un := Un+1 + U M M M

Un Mn

and Un+1 M :=

Un+1 . M n+1

Formula (2.9) determines the approximate discrete Lagrange multiplier λn+1 ≡ H n+1 n+1 n n+1 λ H (U , U ) to enforce the discrete sphere constraint |U M (z)| = 1 for all z ∈ Nh : n+1/2 putting  = U M (z)ϕz in (2.8), and using the algebraic identity n+1 M n+1 dt Un+1 − M = dt U

dt M n+1 n U Mn

(2.10)

for the leading term, together with 

M n+1 dt Un+1 M ,

 h

=

 βz n+1  n+1 M (z) |U M (z)|2 − |UnM (z)|2 = 0 2

∀ z ∈ Nh

to meet the requirement |Un+1 M (z)| = 1 for all z ∈ Nh then leads to definition (2.9). The method is devised in such a way that uniform bounds hold for solutions {Un }n ⊂ Vh : choose  = I h [M n+1 dt Un+1 M ] in (2.8), use the algebraic identity (2.10), and 2 reduced integration for the first and the last term, together with dt |Un+1 M (z)| = 0 for all z ∈ Nh to obtain    1 2 n+1 2  L 2 = ∇Un+1/2 , ∇[Id − I h ] M n+1 dt Un+1 M n+1 dt Un+1 M h + dt ∇U M 2   n+1  n+1 d tM n+1/2 n dt M n ⊗U + + ∇U ,∇ ∇U Mn Mn  dt M n+1 n n+1 n+1 + U , M d U . (2.11) t M Mn h By a standard interpolation result, we may control the first term on the right-hand side by ⎧⎛ ⎫ ⎞1/2 ⎪ ⎪ ⎨    ⎬ 2   2 n+1 ˜ ⎠ + C(d U  +1) . ≤ Ch ∇Un+1/2  L 2 ⎝ 1,2  ∇ M n+1 dt Un+1 t M W M ⎪ ⎪ L 2 (K ) ⎩ K ∈Th ⎭

123

406

L’. Baˇnas et al.

% For the sum, we employ the fact that ∇ 2 Un % K = 0, and an inverse estimate, to obtain the further estimate   ≤ C˜ ∇Un+1/2  L 2 dt Un+1  L 2 + ∇Un  L 2 . For k ≤ k0 (M) sufficiently small, we obtain after summation over all iteration steps 0 ≤ n ≤ N in (2.11), by using Young’s inequality and a discrete version of Gronwall’s Lemma, E(U N +1 ) +

N   k 2 ˜ E(U0 ) + 1 . M n+1 dt Un+1  ≤ C h M 2

(2.12)

n=0

This energy bound is crucial to ensure (subsequence) convergence of iterates from Scheme A to weak solutions of (1.2)1 with (1.9), |u| = M in T , (1.3)–(1.4). Remark 2.1 Note that λn+1 H in Scheme A is not a Lagrange multiplier to the discrete n+1/2 n+1 sphere constraint |U M (z)| = 1 for all z ∈ Nh , since M n+1 U M is used on the right-hand side of (2.8) to obtain the bound (2.12) for the energy E(U N +1 ), uniformly in ( k, h ). The next lemma provides sufficient conditions, which ensure solvability of Scheme A. The proof uses a regularization in a first step to apply Brouwer’s fixed-point theorem; a mesh constraint is then needed to exclude the first case in (2.9). Below, we use the identities, dt Un+1 =

2 n+1/2 {U − Un }, and k

n+1/2

M n+1 U M

= Un+1/2 +

k dt M n+1 n U . 2 Mn (2.13)

Thanks to (2.12), this shows part of the following result. Lemma 2.4 Let Th be a quasiuniform triangulation of  ⊂ Rd , and U0 ∈ Vh such that |U0M (z)| = 1 for all z ∈ Nh , and n ≥ 1. For sufficiently small C˜ ≡ ˜ 2  1, there exists ˜ C(, M, Th ) > 0 independent of k, h > 0 such that k ≤ Ch n+1 n+1 ∈ Vh which satisfies (2.8)–(2.9)2 , enjoys |U M | = 1 for all z ∈ Nh , and U N +1 )+ E(U M

N   k dt Un+1 2h ≤ C M 1 + E(U0M ) 2

(N ≥ 0),

(2.14)

n=0

where C M > 0 depends on T > 0. Proof Step 1. Fix n ≥ 0. For every 18 < ε ≤ continuous mapping Fε : Vh → Vh , with (Fε (W),  ) :=

123

1 4,

and all  ∈ Vh we define the

 2 ) + I (W) + I I (W), (2.15) W − Un ,  h + (∇W, ∇ k

Finite element approximations of harmonic map heat flows and wave maps

407

where

I (W) := −

 z∈Nh

⎛ ⎝

βz

&

dt M n+1 (z) n k dt M n+1 (z) n M n (z) U (z), W(z) + 2 M n (z) U (z)

 βz max |W(z) +

k dt M n+1 (z) n 2 2 M n (z) U (z)| , ε

'



⎞   n+1 k dt M Un ,  ⎠ , ×ϕz W + 2 Mn h  ⎛ k dt M n+1 (z) n ∇W, [W(z) + 2 M n (z) U (z)] ⊗ ∇ϕz  ⎝   I I (W) := − n+1 (z) n (z)|2 , ε βz max |W(z) + 2k dt M U n z∈Nh M (z) ⎞   k dt M n+1 n U , ⎠ . ×ϕz W + 2 Mn h n ˜ Choose  = W,and exploit reduced integration to verify IW (W) ≤ CW h U h . z z Let I I (W) = z I I (W). For every z ∈ Nh , then I IW (W) satisfies

& ' n+1 (z) n    W(z)+ k2 dt M M n (z) U (z), W(z) k dt M n+1 (z) n   ∇W, W(z)+ U (z) ⊗ ∇ϕz n+1 (z) 2 M n (z) n 2 max |W(z)+ k2 dt M M n (z) U (z)| , ε  % k %%  % n+1 n (z)U M (z)ϕz % (2.16) ≤ |∇W|, |∇[W(z)ϕz ]| + %∇ dt M 2 supp(∇ϕz )  ≤ Ch −1 |∇W|, |W(z)| + k max |dt M n+1 | . supp(∇ϕz )

supp(∇ϕz )

  Set  = W in (2.15), with Wh ≥ 1 + C˜ 1 k Un h , and values 0 ≤ k ≤ C˜ 2 h 2 , for some existing 0 < C˜ i ≡ C˜ i (, M), for i = 1, 2. On using Young’s inequality, and the fact that the number of nodes y ∈ Nh such that (∇ϕy , ∇ϕz ) = 0 is bounded independently of h > 0, we obtain for sufficiently small k ≤ k0 (M), ) 2( Wh − Un h Wh + ∇W2L 2 k   −C˜ h −1 ∇W L 2 + Un h Wh + *+ , , ˜ 1 C˜ 1 k Ck 1 n ≥ Wh 1− 2 Wh − 1− U h + ∇W2 ≥ 0, k h 2 2

(Fε (W), W) ≥

and a corollary to Brouwer’s fixed-point theorem [28, p. 37] then implies the existence of Un+1 ∈ Vh such that Fε (Un+1/2 ) = 0, for Un+1/2 := 21 (Un+1 + Un ).

123

408

L’. Baˇnas et al.

Step 2. We proceed by induction to show that Un+1/2 ∈ Vh solves F0 (Un+1/2 ) = 0, ˜ 2. for k ≤ Ch Let n ≥ 1. For all 0 ≤  ≤ n, suppose that U ∈ Vh verifies |UM (z)| = 1 ∀ z ∈ Nh ,

  k dt U j+1 2h ≤ Ct 1 + E(U0M ) . 2 −1

E(UM ) +

j=0

(2.17) In order to satisfy F0 (Un+1/2 ) = 0, it suffices to show for all z ∈ Nh , using parallelogram identity n+1/2

|U M

(z)|2 = 1 −

k2 2 ! 1 |dt Un+1 M (z)| > . 4 2

(2.18)

Thanks to (2.15), the iterate Un+1 := 2Un+1/2 − Un satisfies for all  ∈ Vh 

 ) dt Un+1 ,  + (∇Un+1/2 , ∇ h ' ⎛ & d M n+1 (z) n+1/2 n+1/2 (z) + (∇Un+1/2 , U M (z) ⊗ ∇ϕz )  βz t M n (z) Un (z), U M ⎝ = n+1/2 βz max{M n+1 (z)|U M (z)|2 , ε} z∈Nh n+1/2

× ϕz M n+1 U M

,

 h

.

(2.19)

n+1 Let  = dt Un+1 M (z)ϕz , where z = argmaxy∈Nh |dt U M (y)|. We use identity (2.10), properties of reduced integration, the parallelogram identity, and an inverse estimate to obtain for the first term in (2.19)

dt M n+1 n n+1 U , d U (z)ϕ t z M Mn n+1 2 n ∞ ˜ ≥ βz M n+1 (z)|dt Un+1 M (z)| − Cβz U  L |dt U M (z)|. 

2 βz M n+1 (z)|dt Un+1 M (z)| +

(2.20)

To obtain an upper bound for the second term in (2.19), we use an inverse estimate and (2.10), % % % % % % % n+1 n+1/2 ∞ % ≤ ∇U (z) ⊗ ∇ϕ  U (z) % % ∇Un+1/2 , dt Un+1 %d % ∇ϕz  L 1 z L t M M %  % k % % ∇dt Un+1  L ∞ + ∇Un  L ∞ %dt Un+1 ≤ M (z)% βz |∇ϕz | 2 %  % % % ≤ Cβz h −2 kdt Un+1  L ∞ + Un  L ∞ %dt Un+1 M (z)%   ˜ z h −2 k|dt Un+1 (z)| + [C˜ + 1][k + 1]Un  L ∞ |dt Un+1 (z)|. (2.21) ≤ Cβ M M

123

Finite element approximations of harmonic map heat flows and wave maps

409

For the last term in (2.19), we obtain the upper bound   ˜ z Un  L ∞ + ∇Un+1/2  L ∞ ∇ϕz  L ∞ |dt Un+1 (z)| ≤ Cβ M     ˜ z h −2 k |dt Un+1 (z)| + C˜ + 1 [k + 1] |dt Un+1 (z)|. ≤ Cβ M M 2

(2.22)

Putting (2.20)–(2.22) together, and using assumption (2.17)1 for  = n, we arrive at     k −2 ˜ . 1 − 2 |dt Un+1 M (z)| ≤ C 1 + h h ˜ 2 and sufficiently small C˜ > 0, (2.18) is valid, and hence As a consequence, for k ≤ Ch n+1/2 ∈ Vh satisfies F0 (Un+1/2 ) = 0. Therefore, |Un+1 U M (z)| = 1 for all z ∈ Nh , and the energy bound (2.17)2 is valid again, for all 0 ≤  ≤ n + 1. This closes the induction argument. Let {Un }n≥0 ⊂ Vh be a sequence of iterates computed from Scheme A. In the next step, we show convergence of a subsequence of functions Uk,h : T → Rm , where for all ( t, x ) ∈ [tn , tn+1 ) × , t − tn n+1 tn+1 − t U (x) + U(x), k k − + Uk,h (t, x) := Un (x), Uk,h (t, x) := Un+1 (x), / . n+1/2 Uk,h (t, x) := Un+1/2 (x), and U M k,h (t, x) := U M (x),

Uk,h (t, x) :=

(2.23)

where the limit for (k, h) → 0 is a weak solution of Problem A. Theorem 2.1 Let the assumptions of Lemma 2.4 be valid, with E(u0 ) ≤ C, and U0 → u0 ∈ W 1,2 (, Rm ). There exists a subsequence {Uk,h } of solutions to Scheme A, which for (k, h) → 0 converges weakly in W 1,2 (T , Rm ) to a weak solution u of Problem A, which satisfies (2.2). Proof Step 1. By Lemma 2.4, there exist a convergent subsequence {Uk,h } ⊂ ˜ 2 , and ( k, h ) → 0, W 1,2 (T , Rm ), and u ∈ W 1,2 (T , Rm ) such that for k ≤ Ch + Uk,h , Uk,h , Uk,h  u

  in L ∞ 0, T ; W 1,2 (, Rm ) ,

+ Uk,h , Uk,h , Uk,h → u

in L 2 (T ; Rm ),

∂t Uk,h  ut

in L (T ; R ).

∗

2

(2.24)

m

+ (z)| = M + (z) ≡ M(tn+1 , z), for all z ∈ Nh and t ∈ [tm , tm+1 ), such Note that |Uk,h   . / + 2 | = Ih |M + |2 in . For all K ∈ Th , by piecewise affine-ness of that Ih |Uk,h

123

410

L’. Baˇnas et al.

used finite element functions,     + 2 |Uk,h | − |M + |2 

L 2 (K )

        + 2 + 2 ≤ |Uk,h | − Ih |Uk,h | + Ih |M + |2 − |M + |2  2 L (K )   + 2 + 2 ≤ Ch ∇|Uk,h |  L 2 (K ) + ∇|M |  L 2 (K )       +  +  ) ∇Uk,h +  M + ∇ M +  L 2 (K ) ≤ Ch (Uk,h  2 L (K )   + ˜ (2.25) ≤ Ch ∇Uk,h  L 2 (K ) + h d/2 .

+ Hence, after summation over all K ∈ Th , we conclude holds |Uk,h | → M + almost + everywhere in , for h → 0. Moreover, by again using |Uk,h (z)| = M + (z) for all z ∈ Nh , (2.25), and Lemma 2.4 we infer

T  T  2 2 2       + 2 2 2 |Uk,h | − M  2 ds ≤ C |Uk,h | − |Uk,h |2  2 + |M + |2 − |M|2  2 L

0

L

L

0

2    + 2 + 2 + |Uk,h | − |M |  2 ds L

˜ 3 ≤ Ck

[T /k+1] 

dt Un 2L 2 + dt M n 2L 2



n=1

  + 2  L ∞ (0,T ;L 2 ) + 1 . +C˜ T (h 2 + h d ) ∇Uk,h Hence, |Uk,h | → M a.e. in T for (k, h) → 0, and thus (ii) of Definition 2.1 holds for u : T → Rm . Step 2. By piecewise affine interpolation of {λn+1 H (z)}0≤n≤N we define accordingly z∈Nh

n+1 [λ H ]k,h ∈ C(T ), such that [λ H ]+ k,h (t, z) = λ H (z), for all ( t, z ) ∈ [tn , tn+1 ) × . m By (2.8), we obtain for Uk,h : T → R ,

  /     . +  = [λ H ]+ ∂t Uk,h (t, ·),  h + ∇Uk,h (t, ·), ∇ k,h (t, ·)M (t, ·) U M k,h (t, ·), 

h

  ∞  , Rm . Thanks to a ∧ b, a = 0, the for  ∈ Vh , andall t ∈ (0, T ). Let  ∈ C T  . / choice  = I h U M k,h (t, ·) ∧  (t, ·) then leads to 

  .  / / .  (t, ·) + ∇Uk,h (t, ·), ∇ U M k,h (t, ·) ∧  (t, ·) ∂t Uk,h (t, ·), U M k,h (t, ·)∧ h   / . / . = ∇Uk,h (t, ·), ∇ Id − I h U M k,h (t, ·) ∧  (t, ·) . (2.26)

123

Finite element approximations of harmonic map heat flows and wave maps

411

We use (2.13)2 and compute    / .  ∂t Uk,h , U M k,h ∧  − ∂t u, u ∧ h M + 

,  + − Uk,h M −M  − = ∂t Uk,h , + U − ∂t u, u ∧ ∧ M+ 2M − M + k,h M h +

, +

,   Uk,h Uk,h =: I + ∂t Uk,h , I h ∧ − ∂t Uk,h , I h ∧ M+ M+ h + 

, . / Uk,h + ∂t Uk,h , I h − Id ∧ M+ + 

,  Uk,h u  + ∂t Uk,h , − ∧  + ∂t [Uk,h − u], u ∧ M+ M M =: I + I I + · · · + V I.

(2.27)

Next, we show that after integration in time, all terms I , I I –I I I , I V, . . . , V I indepen˜ dently tend to zero for ( k, h ) → 0. By Lemma 2.4, there exists 0 < C˜ ≡ C(M) < ∞, such that in the limit k → 0, T

˜ ∂t Uk,h  L 2 ( ) ψ ψ  L 2 (T ) → 0. |I | ds ≤ Ck T

0

The properties of ( ·, ·)h , and W 1,2 ()-stability of Ih yield for some existing 0 ≤ C˜ ≡ ˜ C(M) < ∞,

     Uk,h   Ih |I I − I I I | + |I V | ≤ Ch ∂t Uk,h  L 2 ∇I ∧  +  2  M L  



    U k,h   ˜ ∂t Uk,h  L 2 1 + ∇Uk,h  L 2   W 1,∞ . + ∇ ∧  ≤ Ch +  2  M L

Convergence to zero for ( k, h ) → 0 of term V follows from smoothness properties of M > 0, and (2.24)2 . Term V I also tends to zero for ( k, h ) → 0, thanks to (2.24)3 , ˜ 2, and (2.14). Summing up, we find for k ≤ Ch T  T   /   .  ∂t u, u ∧ dt ∀  ∈ C ∞ T , Rm . ∂t Uk,h , U M k,h ∧  dt = lim k,h→0 h M 0

0

(2.28) In the next step, we show that the limit for ( k, h ) → 0 of the second term in (2.26) ˜ 2 is in the case k ≤ Ch

123

412

L’. Baˇnas et al.

  T  T  .  /  ∇u, ∇ u ∧ dt ∇Uk,h , ∇ U M k,h ∧  dt = lim k,h→0 M 0 0   ∀  ∈ C ∞ T , Rm . (2.29) Thanks to (2.13)2 , we restate the second term in (2.26) as 

, + Uk,h M+ − M− − + U ∧ ∇Uk,h , ∇ M+ 2M − M + k,h       ˆ . ± =: I + ∇Uk,h , ∇ Uk,h ∧ M+ M

(2.30)

After integration in time, the term Iˆ vanishes for k → 0, since T

 + 2 ˜ | Iˆ| ds ≤ Ck 1 + ∇Uk,h  L 2 (

T

− 2 + ∇Uk,h  L 2 ( )

 T)

ψ W 1,∞ (T ) → 0. ψ

0

Moreover, the second difference term vanishes on'the right-hand side of (2.30). ' There& &   fore, on using the identities ∇Uk,h , ∇[Uk,h ∧ M ] = ∇Uk,h , Uk,h ∧ ∇ M , and again & ' ' &  ∇u, ∇[u ∧  M ] = ∇u, u ∧ ∇ M almost everywhere, we may conclude    T    ∇Uk,h , ∇ Uk,h ∧ − ∇u, ∇[u ∧ ] dt M M 0

    T  T  . /   ∇Uk,h , [Uk,h − u] ∧ ∇ ∇ Uk,h − u , u ∧ ∇ dt + dt. = M M 0

0

We may use (2.24)2 , and (2.24)3 , respectively, to establish convergence to zero (for ˜ 2 ) for the last two terms. ( k, h ) → 0, with k ≤ Ch It remains to show convergence to zero of%the last term in (2.26). For this purpose, we employ (2.13)2 , and the fact that D 2 Uk,h % K = 0 for all K ∈ Th , + ,

 T    + − M−   . / U M k,h   − + U ∧   dt ∇Uk,h  L 2  ∇ Id − I h k,h  2  M+ 2M − M + L 0   ˜ ∇Uk,h  L 2 ( ) +1 ∇Uk,h  L 2 ( )   W 2,∞ (T )→0 (h → 0), ≤ Ch T T

(2.31)

˜ for some finite C˜ ≡ C(M) > 0. Putting things together in (2.26) then implies that u : T → R3 satisfies item (iii) of Definition 2.1.

123

Finite element approximations of harmonic map heat flows and wave maps

413

Finally, since the trace operator is bounded and linear, it is weakly continuous from W 1,2 (T ) into L 2 (), and we deduce u(0, ·) = u0 in the sense of traces. Therefore, the limit u : T → Rm is a weak solution to Problem A. By weak lower semicontinuity of norms, properties of (·, ·)h , and U0 → u0 ∈ 1,2 W (, Rm ) we conclude from (2.14) that u satisfies (2.2).

3 Wave map to the moving sphere We construct weak solutions of Problem B as proper limits of iterates from Scheme B. We denote D = ( ∂t , ∇ ).   Definition 3.1 Let M ∈ C 2 T ; R+ \ {0} be given, and ( u0 , v0 ) ∈ W 1,2 (, Rm )× L 2 (, Rm ) such that |u0 | = M(0, ·), and u0 , v0 = 21 dtd M 2 (0, ·) a.e. in  ⊂ Rd . We call a function u : T → Rm , such that Du ∈ L 2 (T , Rm ) a weak solution of Problem B, if for all T > 0 there hold (i) u(t, ·)  u0 in W 1,2 (, Rm ), and ut (t, ·) → v0 in L 2 (, Rm ), for t → 0, (ii) |u| = M almosteverywhere in T , and  (iii) for all  ∈ C0∞ [0, T ); W 1,2 (, Rm ) there holds T

T (ut ∧ u,  t ) dt +

− 0

) dt = (v0 ∧ u0 ,  (0, ·)) . (∇u ∧ u, ∇ 0

Similarly to Lemma 2.2, we can show that a mapping u : T → Rm which satisfies the regularity properties given in Definition 3.1, and properties (i), (ii), and for all   ∈ C0∞ [0, T ); W 1,2 (, Rm ) satisfies 

  1 2 2 2 2 2 [∂t − ]M + |∇u| − |ut | ) − u,  dxdt = 0 −(ut ,  t ) + (∇u, ∇ M2

T

(3.1) is a weak solution to Problem B. We plan to construct weak solutions to Problem B as proper limits of iterates which solve the following. Scheme B. Given Un−1 , Un ∈ Vh , find Un+1 such that for all  ∈ Vh , and all z ∈ Nh there hold     n+1 n+1/2 ) = λn+1 M U ,  , dt2 Un+1 ,  + (∇Un+1/2 , ∇ W M h

h

(3.2)

123

414

L’. Baˇnas et al.

⎧ n+1/2 ⎪ 0 if U M (z) = 0, ⎪   ⎪ . / 0 1 n+1 ⎪ (z) k dt M 1 2 n+1 (z)|2 − d Un (z),d Un+1/2 (z) ⎪ |M [1+ ] d t t ⎪ 2 2 t M n (z) ⎪ ⎪ n+1/2 ⎪ ⎨ |M n+1 (z)|2 |U M (z)|2 ( ) k 2 dt M n+1 (z) λn+1 (3.3) dt |dt Un+1 (z)|2 +k|dt2 Un+1 (z)|2 W (z) = ⎪ 8 M n (z) − ⎪ n+1/2 ⎪ n+1 (z)|2 |U 2 |M (z)| ⎪ M  ⎪ ⎪ n+1/2 ⎪ ⊗∇ϕz ∇Un+1/2 (z),U M ⎪ ⎪ ⎩ + else, n+1 2 n+1/2 2 βz |M

n+1/2

where U M

:=

1 2

(z)| |U M

(z)|

  UnM + Un+1 , for UnM := M

Un Mn

and Un+1 M :=

Un+1 . M n+1

We set U−1 := U0 − kV0 , with the given initial velocity V0 , so that V0 = dt U0 , and use the above algorithm for n ≥ 0. Let Vn := dt Un for n ≥ 1. We introduce the discrete energy    1 n 2 V h + ∇Wn 2 . E h Vn , Wn = 2 Formula (3.3) determines the approximate discrete Lagrange multiplier  n+1 n , Un+1 at nodal points to enforce the discrete sphere constraint λn+1 U ≡ λ W W n+1/2 |Un+1 (z)| = 1 for all z ∈ Nh : choose  = U M (z)ϕz in (3.2), and use (2.13)2 , and M k n n+1/2 n+1 U =U − 2 dt U to restate the first term in (3.2) as follows, βz

 dt2 Un+1 (z), Un+1/2 (z) +

k dt M n+1 (z) n U (z) 2 M n (z)



M n+1 (z)   & ' βz k dt M n+1 (z) 2 n+1 n+1/2 d = 1+ U (z), U (z) t 2 M n (z) M n+1 (z) ' 2 n+1 dt M (z) & 2 n+1 βz k n+1 d − U (z), d U (z) . t t 4M n+1 (z) M n (z)

(3.4)

We compute the leading scalar product on the right-hand side of (3.4). ' 1& dt Un+1 (z) − dt Un (z), Un+1/2 (z) k  1& ' 1  n+1 dt |U (z)|2 − dt Un (z), Un−1/2 (z) + kdt Un+1/2 (z) = 2k k  & ' 1  n+1 2 dt |U (z)| − |Un (z)|2 − dt Un (z), dt Un+1/2 (z) = 2k  & ' 1  = dt2 |M n+1 (z)|2 − dt Un (z), dt Un+1/2 (z) . 2 2 n+1 dt Un+1 (z) on the right-hand side of (3.4) is The second scalar product d . / t Uk 2(z), 1 n+1 2 n+1 equal to 2 dt |dt U (z)| + 2 |dt U (z)|2 . Putting things together then leads to (3.3). The steps to come in this section follow the program in Sect. 2: in a first step, we establish the existence of solutions from Scheme B for restricted choices

123

Finite element approximations of harmonic map heat flows and wave maps n+1/2

k = O(h max{1,d/2} ) to exclude the case U M and show (3.3)2 for all n ≥ 0 and z ∈ Nh .

415

(z) = 0 for any z ∈ Nh in (3.3),

d 0 0 Lemma 3.1 Let Th be a quasiuniform0 triangulation 1 of 1d ⊂ 2R , and ( U , V ) ∈ 0 0 0 Vh × Vh such that |U M (z)| = 1, and U (z), V (z) = 2 dt M (0, z) for all z ∈ Nh . ˜ For n ≥ 0, and sufficiently small C˜ ≡ C(, Th ) > 0 independent of k, h > 0 such max{1,d/2} n+1 ˜ , there exists U ∈ Vh , which satisfies Scheme B, satisfies that k ≤ Ch |Un+1 M (z)| = M(z) for all z ∈ Nh , and N   k2     E h V N +1 , U N +1 + dt Vn+1 2h ≤ C M 1 + E h V0 , U0 , 2

(3.5)

n=1

where C M > 0 depends on T > 0. Proof Step to verify the given energy inequality, choose  1. In order  n+1 n+1 2 and recall (2.10); together with dt |Un+1 dt U M  = Ih M M (z)| = 0 for all z ∈ Nh , we find (analogously to (2.11))  1 1 dt dt Un+1 2h + kdt2 Un+1 2h + dt ∇Un+1 2L 2 2 2  n+1 d M n+1 M t = dt2 Un+1 , Un Mn h   /  n+1 . n+1/2 + ∇U , ∇ Id − I h M dt Un+1 M    n+1 dt M dt M n+1 n+1/2 n n ⊗U + + ∇U ,∇ ∇U . Mn Mn

(3.6)

Here, and repeatedly below, we employ the identity dt ϕ n , ψ n = dt ϕ n , ψ n − ϕ n−1 , dt ψ n

(3.7)

to restate the leading term on the right-hand side as follows,  n+1 d M n+1 t n+1 M n U = dt dt U , Mn h  n+1   nd Mn M dt M n+1 M t n n − dt Un , d U + U d . t t M n−1 Mn h Summing up in (3.6) over 1 ≤ n ≤ N , and proceeding as below (2.11) to control the second term on the right-hand side of (3.6) by an approximation argument, we arrive at

123

416

L’. Baˇnas et al.

dt U N +1 2h + k 2

N 

  dt2 Un+1 2h + ∇U N +1 2L 2 ≤ C˜ 1 + dt U1 2h + ∇U1 2L 2

n=1

˜ +Ck

N  

 ∇Un+1 2L 2 + dt Un 2h .

n=1

The discrete Gronwall’s inequality then verifies the energy inequality (3.5). Step 2. Fix n ≥ 0. For every 18 < ε ≤ 41 , and all  ∈ Vh consider the continuous function Fε : Vh → Vh , such that  2  ) (Fε (W),  ) := 2 W + Un−1/2 − 2Un ,  + (∇W, ∇ h k     . / k dt M n+1 n ε ε Iz (W) + I Iz (W) ϕz W + U , , − 2 Mn h

(3.8)

z∈Nh

where Izε := Iza,ε + Izb,ε , with

Iza,ε (W) := −

Izb,ε (W) := −

'   & n+1 (z)  1 2 n+1 (z)|2 − 1 d Un (z), W(z) − Un−1/2 (z) [1 + k2 dt M k t M n (z) ] 2 dt |M

n+1 (z) n 2 max{|W(z) + k2 dt M M n (z) U (z)| , ε}     k 2 dt M n+1 (z) 1 |2W(z)−2Un (z)|2 − |d Un (z)|2 + 1 |2W(z) − 2Un (z) − d Un (z)|2 t t 8 k k M n (z)

n+1 (z) n 2 max{|W(z) + 2k dt M M n (z) U (z)| , ε}  n+1 (z) n ∇W, [W(z) + 2k dt M M n (z) U (z)] ⊗ ∇ϕz ε I Iz (W) :=   . n+1 (z) n 2 βz max |W(z) + k2 dt M M n (z) U (z)| , ε

On setting  = W in (3.8), we find the following upper bounds for the sum in (3.8) for k > 0 sufficiently small,   % % k2 βz % Iza,ε (W)% |W(z)|2 + |dt M n+1 (z)|2 4 z∈Nh      1 1 n ˜ 1 + |dt U (z)| {|W(z)| + 1} ≤C βz 1 + ε k z∈Nh   1 ≤ C˜ dt Un 2h + 1 + 2 2 W2h , 4k ε 

123

(3.9)

Finite element approximations of harmonic map heat flows and wave maps

417

 % % k2 % % βz %Izb,ε (W)% |W(z)|2 + |dt M n+1 (z)|2 4 z∈Nh     1  ≤ C˜ 1 + |dt Un (z)|2 1 + |W(z)|2 βz k 1 + ε z∈Nh     1 1  1+ 1 + |Un (z)|2 + |Un−1 (z)|2 1 + |W(z)|2 . (3.10) ≤C βz k ε 

z∈Nh

The term which involves I Izε (W) is given  and can be controlled as in (2.16).  ) ( For all W =  such that Wh ≥ G Un−i W 1,2 i=1,2 , dt Un h > 0 suf˜ ficiently large, we then find for sufficiently small, positive C˜ ≡ C(), and values ˜ k ≤ Ch, using Young’s inequality, and the fact that the number of nodes y ∈ Nh such that (∇ϕy , ∇ϕz ) = 0 is bounded independently from h > 0, %  % %  % 2  % % (Fε (W), W) ≥ 2 W2h − % Un−1/2 , W % − 2 % Un , W h % + ∇W2L 2 h k 

  ˜ 1 C n 2 ˜ −1 ∇W L 2 Wh −C˜ dt U h + 1 − + W2h − Ch 4k 2 ε2 kε + 

, ˜ ˜ 2 1 2 Ck Ck + 2 Wh − Un−1/2 h − 2Un h ≥ 2 Wh 1− + k h2 ε 4ε 1 + ∇W2 ≥ 0, 2 thanks to (3.5). Brouwer’s fixed-point theorem then implies the existence of Un+1 ∈ Vh , such that Fε (Un+1/2 ) = 0, for Un+1/2 := 21 (Un+1 + Un ). Step 3. We show that Un+1/2 ∈ Vh solves F0 (Un+1/2 ) = 0. For this purpose, it is again sufficient to validate that for all z ∈ Nh holds |Un+1/2 (z)|2 ≥ 1 −

k2 ! 1 |dt Un+1 (z)|2 > . 4 2

(3.11)

˜ d/2 : Let n ≥ 1. For all We proceed by induction to show this result, for k ≤ Ch   0 ≤  ≤ n, suppose that U ∈ Vh satisfies |U M (z)| = 1, for all z ∈ Nh , and −1   k2    E h V , U + dt Vm+1 2h ≤ Ctn 1 + E h (V0 , U0 ) . 2

(3.12)

m=1

The iterate Un+1 := 2Un+1/2 − Un ∈ Vh from Step 2 satisfies for all  ∈ Vh , + 3 ,    n+1 2 n+1 n+1/2 n+1 n+1/2 ) = dt U ,  + (∇U , ∇ λW,ε,i M UM ,  , (3.13) h

i=1

h

λn+1 W,ε,i (z),

where for i = 1, 2, 3 are each equal to the corresponding terms in (3.3) for n+1 λW (z), where |Un+1/2 (z)| in each denominator is replaced by max{|Un+1/2 (z)|, ε},

123

418

L’. Baˇnas et al.

respectively. We choose  = I h [M n+1 dt Un+1 M ], and use (2.10). Terms from the right-hand side of (3.13) are dealt with separately in the next step. We use kdt Un+1 = 2Un+1/2 − 2Un below. Let k > 0 be sufficiently small.   n+1 n+1 n+1/2 n+1 λn+1 M U , M d U t W,ε,1 M M h , + ( ) n n+1 C˜ (1 + |dt U |) 1 + |dt U | + |dt Un | n+1/2 n , |U ≤ | + |U | k max{|Un+1/2 |, 2ε } h 

  ˜ ˜ C C ≤ 1+ dt Un h dt Un+1 h + dt Un h , k ε   n+1 n+1/2 λn+1 UM , M n+1 dt Un+1 (3.14) W,ε,2 M M h + , . /  Ck |dt Un+1 |2 + |dt Un |2 1 + |dt Un+1 (z) | ≤ max{|Un+1/2 |, 2ε } h    1 n+1 2 n 2 , ≤ C dt U  L 2 + dt U  L 2 1 + ε   n+1 n+1/2 λn+1 UM , M n+1 dt Un+1 W,ε,3 M M h    1 h −1 ∇Un+1/2  L 2 1 + dt Un+1 h . ≤ C˜ 1 + ε The first two terms in (3.13) are may be restated by means of (2.10) as given in the next equality.  k 1  dt dt Un+1 2h + ∇Un+1 2L 2 + dt2 Un+1 2h 2 2     n n+1  n+1 d M d t tM n n n+1 dt M n d U d = dt Un , U + d − d U , U t t t t M n−1 Mn Mn h h    dt M n+1 dt M n+1 n+1/2 n n ⊗U + + ∇U ,∇ ∇U Mn Mn + 3 ,  n+1 n+1 n+1 n+1/2 n+1 + λW,ε,i M UM ,M dt U M . i=1

h

Thanks to (3.14), summation over 1 ≤  ≤ n yields to 



˜ ˜ Ck 1 1 k 1 Ck k n+1 2 1− − dt U h + 1− − ∇Un+1 2L 2 + dt2 Un+1 2h k 8 h k 8 h 2

   ˜ 1 1 n 2 ˜ + C Un  L ∞ dt Un 2h + C˜ k + 1 Un 2 2 + C∇U ˜ ≤ L 2 . 1 + + Ck L k 4 ε k

123

Finite element approximations of harmonic map heat flows and wave maps

419

˜ max{1,d/2} , by inverse Consequently, for sufficiently small C˜ > 0, such that k ≤ Ch estimate, and induction assumption (3.12), there holds k2 dt Un+1 2L ∞ ≤ Ck 2 h −d dt Un+1 2h ≤ Ck 2 h −d . 2

(3.15)

˜ max{1,d/2} . Therefore, F0 (Un+1/2 ) = Hence, assertion (3.11) is valid for values k ≤ Ch n+1 0 holds, and hence |U M | = 1, as well as (3.12) for 0 ≤  ≤ n + 1, and hence the induction argument is complete.   ˜ max{1,d/2} The proof shows that the first case in (3.3) does not apply in case k ≤ Ch holds. In the next step, we verify convergence of iterates {Un } of Scheme B towards weak solutions of Problem B for ( k, h ) → 0. In the sequel, we use the notation (2.23), but drop sub-indices k, h. For  ∈ C0∞ ([0, T ); Vh ), on putting Vn+1 = dt Un+1 , we may rewrite the first equation of Scheme B as T

.     / +  − λ+ − (Vt ,  )h + ∇U, ∇ W M U M ,  h dt = 0.

(3.16)

0

The key to validate Definition 3.1, (iii) is to obtain a reformulation .  uses / of (3.16), which ∞ [0, T ); wedge products: in particular, we choose  = I U ∧  , with  ∈ C h M 0  C ∞ (, Rm ) . Here, we use again the notation from (2.23). Then, the right-hand side of (3.16) which involves the Lagrange multiplier drops out. In fact, on using U M = + − U + k2 MM −−M U− , cf. (2.13)2 , we can easily adopt the proof of Lemma 4.2 in [12] M+ M+ to the present case. Lemma 3.2 Suppose that the assumptions of Lemma 3.1 are valid. There holds for all  ∈ C0∞ ([0, T ); C ∞ (, Rm )), % T % %   %% % .   / % I h [U M ∧  ] dt − V0 , U M (0, ·) ∧  (0, ·) %% Ut , [U M ∧  ]t h + ∇U, ∇I % h% % 0 % T % % %   %  +  % % ˜ 1/2 1 + E h (V0 , U0 )    L ∞ (T ) . ≤% V − V, [U M ∧  ]t dt %% + Ck % % 0

(3.17)    The bound (3.5) yields the existence of u ∈ L ∞ 0, T ; W 1,2 (, Rm ) ∩ W 1,∞ 0, T ;  L 2 (, Rm ) , and a convergent subsequence {Um } (not relabeled), such that for k ≤

123

420

L’. Baˇnas et al.

˜ max{1,d/2} and ( k, h ) → 0, Ch ∗

U, U+ , U  u U, U+ , U → u + ∗

Ut , V, V  ut

  in L ∞ 0, T ; W 1,2 (, Rm ) , in L 2 (, Rm ),   in L ∞ 0, T ; L 2 (, Rm ) .

(3.18)

In order to show that |U| → M a.e. in T for ( k, h ) → 0, we use |U+ (z)| = M + for all z ∈ Nh , and then follow arguments in (2.25). Being provided with these properties, we can follow the arguments given in the proof of [12, Theorem 4.1], namely: (1) to show that the first term in (3.17) may be restated as a controllable, asymptotically vanishing (for ( k, h ) → 0) perturbation of   T  − 0 (Ut , U ∧ M )h dt, and t

    T  T    Ut , U ∧ ut , u ∧ dt lim = k,h→0 M t h M t 0 0   ∀  ∈ C0∞ [0, T ); C ∞ (, Rm ) . (2) to show that the second term in (3.17) may be restated as a controllable, asympT  . / totically vanishing (for ( k, h ) → 0) perturbation of − 0 ∇U, ∇ U M ∧  dt, and T lim

k,h→0

. /  ∇U, ∇ U M ∧  dt =

0

  T   ∇u, ∇ u ∧ dt M 0   ∀  ∈ C0∞ [0, T ), C ∞ (, Rm ) .

Here, we can use the arguments below (2.29) and up to (2.31). (3) to verify assertion (i) of Definition 3.1, we follow the arguments in [12, Step 3 of the proof of Theorem 4.1]. . / Let E (v, w) = 21 v2 + ∇w2 . On using (3.5), the properties (3.18), those of (·, ·)h , and assumptions regarding convergence of initial data (U0 , V0 ) for h → 0 yields for almost every t ∈ (0, T ), E (ut (t, ·), u(t, ·)) ≤ lim inf E (Ut (t, ·), U(t, ·)) k,h→0

= lim inf E h (Ut (t, ·), U(t, ·)) ≤ C M (1 + E (u0 , v0 )) . k,h→0

(3.19)

Thus we have shown the following result. Theorem 3.1 Let the assumptions of Lemma 3.1 be valid, and U0 → u0 ∈ W 1,2 (, Rm ) as well as V0 → v0 ∈ L 2 (, Rm ), for h → 0. There exists a subsequence {Uk,h } of solutions to Scheme B, which for ( k, h ) → 0 converges weakly in W1,2 (T , Rm ) to a weak solution u of Problem B, which satisfies (3.19).

123

Finite element approximations of harmonic map heat flows and wave maps

421

4 Computational studies 4.1 Solution of the nonlinear systems For every n ≥ 0, we solve the nonlinear problems (2.8) resp. (3.2) by a fixed point algorithm with projection. For Scheme A and n = 1, . . . , N it can be described as follows. Algorithm A1 : Let  > 0 be given. Set u0 = u∗0 = Un and iterate for l = 0, . . . (1) Find ul+1 ∈ Vh such that for all  ∈ Vh (i) 

 λ˜ lH ul+1 ,  h   n+1 ) + k2 λ˜ lH MM n Un ,  . = (Un ,  )h − k2 (∇Un , ∇

 l+1  u , h +

k 2

 l+1   − ∇u , ∇

k 2

h

(ii) For all z ∈ Nh set ul+1 ∗ (z) =

M n+1 (z) l+1 u (z). |ul+1 (z)|

(2) Stop the iteration when l max |ul+1 ∗ (z) − u (z)| ≤ .

z∈Nh

(3) Set Un+1 = ul+1 ∗ and proceed to the next time level. The Lagrange multiplier λ˜ lH is defined as in (2.9), with the unknown values of Un+1 replaced by the values of ul∗ for l ≥ 1. We proceed similarly for Scheme B, where a fixed-point iteration scheme for n = 2, . . . , N (with an obvious modification for n = 1) takes the following form. Algorithm B1 : Let  > 0. Set u0 = u∗0 = Un and iterate for l = 0, . . . (1)

(i) Find ul+1 ∈ Vh such that for all  ∈ Vh  k  1  l+1  k  l+1  − u , + ∇u , ∇ λ˜ lW ul+1 ,  h h k 2 2   1 n k n−1 n  2U − U ,  − ∇U , ∇ = h k 2  M n+1 n k U , . + λ˜ lW 2 Mn h (ii) For all z ∈ Nh set ul+1 ∗ (z) =

M n+1 (z) l+1 u (z). |ul+1 (z)|

123

422

L’. Baˇnas et al.

(2) Stop the iterations when l max |ul+1 ∗ (z) − u (z)| ≤ .

z∈Nh

(3) Set Un+1 = ul+1 ∗ and proceed to the next time level. The Lagrange multiplier λ˜ lW is defined as in (3.3), with the unknown values of Un+1 replaced by the values of ul∗ for l ≥ 1. Remark 4.1 (Projection step) In our computational studies, we find that the projection step is not necessary for the convergence of the above fixed-point algorithms for choices k = O(h 2 ). In this case, numerical experiments indicate that the fixed-point algorithms with and without the projection step produce indistinguishable results, for sufficiently small . In contrast, we observe that convergence of the projection algorithm does not require the time step size restriction k = O(h 2 ), unlike the case of fixed point iterations without the projection step. Moreover, the algorithms with projection converge after significantly fewer iterations (typically 4 up to 6 iterations) than the fixed-point algorithms without the projection step, for parameters k chosen such that the latter converges. n+1 (z) Remark 4.2 (Sphere constraint) The discrete sphere constraint |ul+1 ∗ (z)| = M l+1 the constraint is only satisholds trivially for l ≥ 0, ∀z ∈ Nh . In general for u fied approximately with error O(). Therefore, to obtain a good approximation of the nonlinear problems (2.8) resp. (3.2) the stopping threshold  needs to be sufficiently small. We set  = 10−14 which results in maxz∈Nh ||ul+1 (z)| − M n+1 (z)| ≈ 10−15 . Also note, that alternatively one could take Un+1 = ul+1 in the above fixed-point algorithms which would result in slightly higher overall error in the sphere constraint maxz∈Nh ||Un+1 (z)| − M n+1 (z)| ≈ 10−12 .

4.2 2D numerical experiments Our code is based on the finite element toolbox ALBERT [25]. In order to minimize the influence of the error caused by the approximate solution of the linear algebraic systems on the constraint |U| = M, we employ a direct solver, see e.g. [19]. We investigate the influence of the space–time varying constraint |u(t, x)| = M(t, x) on the evolution of the solution. Numerical studies for M ≡ 1 using similar numerical techniques as in this paper can be found in [12]; see also [11] for corresponding numerical experiments for wave maps. Videos of some of the computed evolutions can be downloaded from [3]. 2  We use an initial condition on  = − 21 , 21 which motivates a finite-time blow-up for M ≡ 1, ⎧ for |x| ≤ 1/4 ⎨ (4x A, A2 − 4|x|2 )/(A2 + 4|x|2 ) 0 u (x) = (−4x A, A2 − 4|x|2 )/(A2 + 4|x|2 ) for 1/4 ≤ |x| ≤ 1/2 ⎩ (−x, 0)/|x| for 1/2 ≤ |x|

123

(4.1)

Finite element approximations of harmonic map heat flows and wave maps

423

Fig. 1 Example 4.1: sequence of meshes

Fig. 2 Example 4.1: zoom at the solution near the origin at times t = 0.0215(h min = 1/128), 0.0325(h min = 1/256), 0.0395(h min = 1/512), 0.044(h min = 1/1024) and exact stationary solution

where A = (1−4|x|)4 . For general M(t, x), the initial condition is modified as follows, U0 (x) = M(0, x)u0 (x).

(4.2)

4.2.1 Harmonic maps Example 4.1 In this example we examine the bubbling off of spheres near the singularity for harmonic maps with M ≡ 1. The bubbling off of spheres in the continuous case is shown by an appropriate zooming near the singularity, cf. Sect. 2.1. Naturally, the amount of zooming at a discrete level is restricted by the minimum mesh size h around the point where the blow-up occurs. To be able to show the evolution of the solution at different scales we computed the problem on a sequence of adapted meshes. The meshes were con1 around the origin, structed by a gradual local refinement inside circles of radius 22+i 1 i = 1, . . . , 4, where h min = 26+i , see Fig. 1. Note that in the case of singular solutions, √

we have for a mesh with a given h, that ∇U L ∞ = h2min2 , see Fig. 3. The snapshots of the solution in Fig. 2 reveal that by zooming on the solution at the origin we recover solutions which resemble the continuous stationary axisymmetric one of the harmonic map equation given by x → arctan(c|x|), where c is a constant. While the blow-up time for different meshes are quite different, we observe convergence towards a finite blow-up time t < 0.0463. Example 4.2 We compute the evolution starting from the initial condition (4.2) for the following choices of M: (i) M1 (x) ≡ 1

123

424

L’. Baˇnas et al. 3000 h_min=1/128 h_min=1/256 h_min=1/512 h_min=1/1024

2500

2000

1500

1000

500

0

0

0.01

Fig. 3 Example 4.1: evolution of

0.02

0.03

0.04

0.05

0.06

∇U L ∞

(ii) ⎧ ∗ ⎨m M2 (x) = 1 ⎩ ∗ m + 5(1 − m ∗ )(|x| − 0.2)

|x| < 0.2 |x| > 0.4 0.2 ≤ |x| ≤ 0.4

with m ∗ = 2. (iii) M3 (x) defined as (ii) with m ∗ = 0.001. (iv) ⎧ ⎨ M1 (x) M4 (t, x) = M2 (x) ⎩ M1 (x) +

(t−0.001) (0.03−0.001) [M2 (x) −

M1 (x)]

t < 0.001 t > 0.03 0.001 ≤ t ≤ 0.03

M1 (x)]

t < 0.001 t > 0.03 0.001 ≤ t ≤ 0.03

(v) ⎧ ⎨ M1 (x) M5 (t, x) = M3 (x) ⎩ M1 (x) +

(t−0.001) (0.03−0.001) [M3 (x) −

We use a non-homogeneous Dirichlet boundary condition with boundary values prescribed by the initial condition. We display the evolution of ∇U L ∞ for (i)–(v) in Fig. 4.

123

Finite element approximations of harmonic map heat flows and wave maps

425

400

180 160

M_1, h=1/64 M_1, h=1/128 M_2, h=1/64 M_2, h=1/128 M_3, h=1/64 M_3, h=1/128

140 120

350 300

M_4, h=1/64 M_4, h=1/128 M_5, h=1/64 M_5, h=1/128

250

100

200

80 150

60

100

40

50

20 0

0

0.01

0.02

0.03

0.04

0.05

0.06

0

0

0.01

0.02

0.03

0.04

0.05

0.06

Fig. 4 Example 4.2: evolution of ∇U L ∞

Fig. 5 Example 4.2 (i): solution for h = 1/16 at time t = 0, 0.006, 0.0265, 0.06

Fig. 6 Example 4.2 (ii): detail of U at the origin for h = 1/64 at time t = 0, 0.0055, 0.25, 0.06

The evolution in (i) motivates a finite time blow-up behavior, see Fig. 5. The results are analogous to those obtained in [7]. In the case (ii) no singularities arise, and the solution converges to a steady state where vectors surrounding the origin with |U| = 2 point upwards, see Fig. 6. The solution around the origin has a greater influence on the overall evolution, i.e., the evolution towards a blow-up becomes energetically more expensive. The case (iii) leads to an early formation of singularities, see Fig. 7. Similar to the previous case, the ‘solution part’ with greater magnitude dominates the evolution, the vectors near the origin with smaller magnitude are forced to rotate quickly downwards, towards an energetically favourable state, which motivates an early blow-up occurrence. The evolutions of (iv), (v) are less obvious, and computations depend on the mesh size and the time evolution of M(t, ·). Based on Fig. 4 we can summarize the results for (iv), (v) as follows. In the case (iv), the value of M(t, x) grows in time near the origin, where a singularity formation is expected. For h = 1/64, the evolution toward blow-up prevails over the influence of growing M. At the end of the computation all

123

426

L’. Baˇnas et al.

Fig. 7 Example 4.2 (iii): detail of U at the origin for h = 1/64 at time t = 0, 0.0025, 0.0055, 0.06

vectors point downwards for h = 1/64. However, as we refine the discretization to h = 1/128, the influence of M prevails, and no blow-up of the solutions is observed since the vectors around the origin are forced to turn upwards. The results indicate that for k, h → 0 no blow-up occurs. In the case (v) the M decreases in time around the origin, which accelerates the evolution towards blow-up and switching. Note that since the blow-up time t ∗ differs for different mesh sizes h, the value of ∇U(t ∗ , ·) L ∞ also depends on the value of M(t ∗ , ·). We refer the reader to [4] for similar experiments with space–time varying M. Finally, we conclude that the numerical experiments for (iv), (v) show a complex behaviour of the solution depending on space–time evolution of M and the mesh parameter h. Further numerical and theoretical studies might be needed at this place for a better understanding of the influence of discretization errors. 4.2.2 Wave maps Example 4.3 We compute the evolution starting from the initial condition (4.2) for the following choices of M(t, x): (i) (ii) (iii) (iv)

M1 (x) ≡ 1, M2 (x) defined as in Example 4.2, M3 (x) defined as in Example 4.2, We define 2 = M(t)

⎧ ⎨1 m∗ ⎩1 π π −3 ∗ −3 ∗ 2 [cos( 0.199 (t −10 ))+m cos(π + 0.199 (t − 10 ))+1+m ]

t < 0.001, t > 0.2, 0.01 ≤ t ≤ 0.2.

We take m ∗ = 2 and set ⎧ 2 ⎨ M(t) M4 (t, x) = 1 ⎩1 2 [ M(t) cos( 2

π 0.2 (|x| − 0.2))+cos(π

+

π 2 0.2 (|x| − 0.2))+1+ M(t)]

|x| < 0.2 |x| > 0.4 0.2 ≤ |x| ≤ 0.4

(v) M5 (t, x) defined as in (iv) with m ∗ = 0.05. We use a homogeneous Neumann boundary condition in this example. The evolution of ∇U L ∞ is depicted in Fig. 8; for M1 ≡ 1, the results fit to those in [11]. Computed discrete energies decay at certain instances coming along with a

123

Finite element approximations of harmonic map heat flows and wave maps 400

400 M_1, h=1/32 M_1, h=1/64 M_2, h=1/32 M_2, h=1/64 M_3, h=1/32 M_3, h=1/64

350 300

300 250

200

200

150

150

100

100

50

50 0

0.2

0.4

0.6

0.8

M_4, h=1/32 M_4, h=1/64 M_5, h=1/32 M_5, h=1/64

350

250

0

427

1

0

0

0.2

0.4

0.6

0.8

1

Fig. 8 Example 4.3: evolution of ∇U L ∞

Fig. 9 Example 4.3 (i): solution at time t = 0, 0.14, 0.45, 1

Fig. 10 Example 4.3 (ii): solution at time t = 0, 0.35, 0.45, 1

blow-up and switching process, while the energy is approximately conserved during smooth evolution. We observe some small numerical dissipation caused by the implicit Euler method—which was crucial from an analytical standpoint in Sect. 3 to have sufficient compactness properties for sequences obtained by Scheme B to identify limits as weak solutions to (1.5)–(1.7); see also [12]. As expected, the evolution (i) motivates a blow-up of the solution, see Fig. 9. In the case (ii), the evolution still leads to a formation of singularities; however there is no loss of energy, and the vector at the origin does not change its orientation, see Fig. 10. The case (iii) evolves singularities, but they can not be deduced from the graph of ∇U L ∞ , see Fig. 8, due to the small values of M near the origin. The evolution is displayed in Fig. 11. Eventually, the vector at the origin rotates and points in opposite direction than at t = 0. This process, however, is slow for h = 1/16.

123

428

L’. Baˇnas et al.

Fig. 11 Example 4.3 (iii): solution at time t = 0, 0.12, 0.355, 1

Fig. 12 Example 4.3 (iv): solution at time t = 0, 0.1, 0.3, 1

Fig. 13 Example 4.3 (v): solution at time t = 0, 0.1, 0.5, 1

The case (iv) is defined is such a way that M(0, x) = 1 everywhere and then the values of M increase smoothly in a circle around the origin until M = 2 in the circle. The evolution can be seen in Fig. 12. Case (v) shows a complex evolution, see Fig. 13. The computation indicates formation of a singularity at the origin. We found that the wave maps computations were more sensitive to the smoothness of M in time, otherwise discontinuous jumps in time occurred in the values of kinetic energy. Remark 4.3 From the results for the given set of numerical experiments it can be seen that the evolutions governed by the harmonic map equation and the wave map equation retain some characteristics of their linear counterparts, the parabolic and wave equations, respectively. The harmonic map equation results in a diffusion type evolution while the wave map equation leads to solutions containing waves which are reflected from the origin.

123

Finite element approximations of harmonic map heat flows and wave maps 90

45

M_1, h=1/8 M_1, h=1/16 M_2, h=1/8 M_2, h=1/16 M_3, h=1/8 M_3, h=1/16

80 70

429

40 35

60

30

M_1, h=1/8 M_1, h=1/16 M_2, h=1/8 M_2, h=1/16 M_3, h=1/8 M_3, h=1/16

50 25

40

20

30 20

15

10

10

0

0

0.2

0.4

0.6

0.8

1

1.2

5

0

0.2

0.4

0.6

0.8

1

1.2

Fig. 14 Example 4.4: evolution of the discrete energy and ∇U L ∞

4.3 3D numerical experiments In this section, we study a 3D problem with Dirichlet boundary condition u D (x) = x , and  = (−0.5, 0.5)3 . The infinite energy initial data for all u(t, x)|∂ = M(x) |x| the experiments was (cf. [6]) u0 (x)|\∂ = (M(x), 0, 0) on  \ ∂, and u0 (x) = u D (x) on ∂. Example 4.4 We study long time behaviour of the solution governed by (1.2). We compute the problem for the following choices of M(x): (i) M1 (x) ≡ 1, (ii) M2 (x), with M2 defined as in Example 4.2, (iii) M3 (x), with M3 defined as in Example 4.2. The evolution of computed discrete energies and ∇U L ∞ for (i)–(iii) is displayed in Fig. 14. The snapshots of the computed solution at x3 = 0 for (i) are depicted in Fig. 15. It is known that a steady state solution in the case (i) is the mapping x , with a singularity at x = 0. Our computed solution at the steady state limit u(x) = |x| is similar to the results from [6], where the steady state solution is computed by an energy minimization technique. The steady state solutions at x3 = 0 for the cases (ii), (iii) are displayed in Fig. 16, respectively. Since M is symmetric with respect to 0 and M|∂ ≡ 1 for (ii), (iii), it can be expected that the steady state solution for h → 0 x . The experiments indicate that the retains a symmetry similar to the symmetry of |x| x . steady state solutions for space-varying M(x) converge to the map u(x) = M(x) |x| It is interesting to note that in the case (ii), the distance of the inhomogeneity to 0 (which for h → 0 should move towards 0) increases. This can be explained by the fact that such configuration is energetically preferable on the discrete level, i.e., a homogeneous solution close to 0 results in lower discrete energies. It is expected that for suitably small mesh parameter h, the inhomogeneity moves towards 0. The case (ii) contrasts with the case (iii) were M is very small around 0, thus allowing the singularity to be located exactly at 0, with the orientation of u(0) having little impact on the overall energy.

123

430

L’. Baˇnas et al.

Fig. 15 Example 4.4 (i): solution at time t = 0, 0.05, 0.1, 0.15, 0.2, 10

Fig. 16 Example 4.4 (ii)–(iii): stationary solution

Example 4.5 We study long time behaviour of the solution governed by (1.5). The problem was computed for the following choices of M(x): (i) M1 (x) ≡ 1, (ii) M2 (x), with M2 defined as in Example 4.2, (iii) M3 (x), similar as M3 defined as in Example 4.2, with m ∗ = 0.01. The evolution of computed discrete energies and ∇U L ∞ for (i)–(iii) is displayed in Fig. 17. The snapshots of the computed solution at x3 = 0 for (i) are depicted in Fig. 18; an inhomogeneity quickly forms at the left side of the domain, then travels

123

Finite element approximations of harmonic map heat flows and wave maps 90

55 M_1, h=1/8 M_1, h=1/16 M_2, h=1/8 M_2, h=1/16 M_3, h=1/8 M_3, h=1/16

80 70 60

M_1, h=1/8 M_1, h=1/16 M_2, h=1/8 M_2, h=1/16 M_3, h=1/8 M_3, h=1/16

50 45 40 35

50

30

40

25

30

20

20

15

10

10

0

431

5 0

20

40

60

80

100

0

20

40

60

80

100

Fig. 17 Example 4.5: evolution of the discrete energy and ∇U L ∞

Fig. 18 Example 4.5 (i): solution at time t = 0, 1, 1.5, 2, 2.5, 100

slightly to the right side of the domain and subsequently returns to the origin. We observe a rapid release of the initial (infinite) energy, which stabilizes at larger times— which seems to happen because a singularity is always present here. The observed evolution differs from the one obtained by the harmonic maps equation, cf. Remark 4.3. References 1. Alouges, F.: A new algorithm for computing liquid crystal stable configurations: the harmonic mapping case. SIAM J. Numer. Anal. 34, 1708–1726 (1997) 2. Alouges, F., Jaisson, P.: Convergence of a finite elements discretization for the Landau Lifshitz equations. Math. Models Methods Appl. Sci. 16, 299–316 (2006) 3. Baˇnas, L’.: http://www.ma.hw.ac.uk/~lubomir/research.html 4. Baˇnas, L’., Prohl, A., Slodiˇcka, M.: Modeling of thermally assisted magnetodynamics. SIAM J. Numer. Anal. 47, 551–574 (2008)

123

432

L’. Baˇnas et al.

5. Baˇnas, L’., Bartels, S., Prohl, A.: A convergent implicit finite element discretization of the MaxwellLandau-Lifshitz-Gilbert equation. SIAM J. Numer. Anal. 46, 1399–1422 (2008) 6. Bartels, S.: Stability and convergence of finite-element approximation schemes for harmonic maps. SIAM J. Numer. Anal. 43, 220–238 (2005) 7. Bartels, S., Prohl, A.: Constraint preserving implicit finite element discretization of harmonic map flow into spheres. Math. Comp. 76, 1847–1859 (2007) 8. Bartels, S., Prohl, A.: Convergence of an implicit finite element method for the Landau-Lifshitz equation. SIAM J. Numer. Anal. 44, 1405–1419 (2006) 9. Bartels, S., Prohl, A.: Stable discretization of scalar and constrained vectorial Perona-Malik equation. Interfaces Free Boundaries 9, 431–453 (2007) 10. Barrett, J.W., Bartels, S., Feng, X., Prohl, A.: A convergent and constraint-preserving finite element method for the p-harmonic flow into spheres. SIAM J. Numer. Anal. 45, 905–927 (2007) 11. Bartels, S., Feng, X., Prohl, A.: Finite element approximations of wave maps into spheres. SIAM J. Numer. Anal. 46, 61–87 (2008) 12. Bartels, S., Lubich, C., Prohl, A.: Convergent discretization of heat and wave map flows to spheres using approximate discrete Lagrange multipliers. Math. Comp. 78, 1269–1292 (2009) 13. Bizo´n, P., Chmaj, T., Tabor, Z.: Dispersion and collapse of wave maps. Nonlinearity 13, 1411–1423 (2000) 14. Bizo´n, P., Chmaj, T., Tabor, Z.: Formation of singularities for equivariant (2 + 1)-dimensional wave maps into the 2-sphere. Nonlinearity 14, 1041–1053 (2001) 15. Chang, K.C., Ding, W.Y., Ye, R.: Finite-time blow-up of the heat flow of harmonic maps from surfaces. J. Differ. Geom. 36, 507–511 (1992) 16. Chen, Y.M., Struwe, M.: Existence and partial regularity results for the heat flow for harmonic maps. Math. Z. 201, 83–103 (1989) 17. Chen, Y.: The weak solutions to the evolution problems of harmonic maps. Math. Z. 201, 69–74 (1989) 18. Coron, J.-M., Ghidaglia, J.-M.: Explosion en temps fini pour le flot des applications harmoniques. CR. Acad. Sci. Paris Ser. I 308, 339–344 (1989) 19. Davis, T.A., Duff, I.S.: An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM J. Matrix Anal. Appl. 18, 140–158 (1997) 20. Grotowski, J.F., Shatah, J.: A note on geometric heat flows in critical dimensions. Preprint (2006). http://math.nyu.edu/faculty/shatah/preprints/gs06.pdf 21. Hairer, E., Lubich, C., Wanner, G.: Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edn. Springer, New York (2006) 22. Krieger, J., Schlag, W., Tataru, D.: Renormalization and blow up for charge one equivariant critical wave maps. Invent. Math. 171, 543–615 (2008) 23. Kruzik, M., Prohl, A.: Recent developments in the modeling, analysis, and numerics of ferromagnetism. SIAM Rev. 48, 439–483 (2006) 24. Rodnianski, I., Sterbenz, J.: On the formation of singularities in the critical O(3) σ -model. Preprint (arXiv-series) (2006) 25. Schmidt, A., Siebert, K.G.: ALBERT—software for scientific computations and applications. Acta Math. Univ. Comenian. (N.S.) 70, 105–122 (2000) 26. Shatah, J.: Weak solutions and development of singularities in the SU (2) σ model. Comm. Pure Appl. Math. 41, 459–469 (1988) 27. Shatah, J., Struwe, M.: Geometric Wave Equations. New York University, Courant Institute of Mathematical Sciences, New York (1998) 28. Showalter, R.E.: Monotone Operators in Banach Space and Nonlinear Partial Differential Equations. AMS (1997) 29. Struwe, M.: Geometric evolution problems. IAS/Park City Math. Series 2, 259–339 (1996) 30. Struwe, M.: On the evolution of harmonic maps of Riemannian surfaces. Math. Helv. 60, 558–581 (1985) 31. Tang, B., Sapiro, G., Caselles, V.: Diffusion of generated data on non-flat manifolds via harmonic maps theory: the direction diffusion case. Int. J. Comput. Vis. 36, 149–161 (2000) 32. Tang, B., Sapiro, G., Caselles, V.: Color image enhancement via chromaticity diffusion. IEEE Trans. Image Proc. 10, 701–707 (2001) 33. Tataru, D.: The wave maps equation. Bull. Am. Math. Soc. 41, 185–204 (2004) 34. Vese, L.A., Osher, S.J.: Numerical methods for p-harmonic flows and applications to image processing. SIAM J. Numer. Anal. 40, 2085–2104 (2002)

123