- Email: [email protected]
A numerical scheme for a Maxwell–Landau–Lifshitz–Gilbert system
Applied Mathematics and Computation 158 (2004) 79–99 www.elsevier.com/locate/amc
A numerical scheme for a Maxwell–Landau–Lifshitz–Gilbert system Mari an Slodicka
*,1,
LÕubomır Ba nas
2
Department of Mathematical Analysis, Faculty of Engineering, Ghent University, Galglaan 2, B-9000 Ghent, Belgium
Abstract We consider MaxwellÕs equations together with a nonlinear dissipative magnetic law described by the Landau–Lifshitz–Gilbert equation. The paper is devoted to the numerical study of this problem. We introduce two algorithms for the time discretization. The modulus of magnetization is conserved in both cases. Assuming a sufficiently smooth electromagnetic field, we derive the error estimates for both approximation schemes. 2003 Elsevier Inc. All rights reserved. Keywords: Electromagnetic field; Maxwell–Landau–Lifshitz–Gilbert system; Time discretization; Hysteresis
1. Introduction Simulations of magnetic behavior based on the micromagnetic theory exhibits hysteresis. These systems have a highly nonlinear character. The mean field theory of ferromagnetism was proposed by Weiss [27], where the hysteresis was based on a positive feedback effect. Weiss assumed that ferromagnetic materials have a spontaneous magnetization M, even in the absence of any
*
Corresponding author. E-mail addresses: [email protected] (M. Slodicka), [email protected] (L. Ba nas). 1 Supported by the BOF/GOA-project no. 12 052 499 of Ghent University. 2 Supported by the project no. 174 IW 300 of Ghent University. 0096-3003/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.08.065
80
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
applied magnetic field H. The fact that they can look as demagnetized on a macromagnetic scale was assign to the presence of a microscopic domain structure. This issue was later described by Landau and Lifshitz [16,17]. Let a ferromagnetic material occupy a bounded domain X R3 with a Lipschitz boundary C. The time-dependent magnetic variables are related as follows: B ¼ l0 ðH þ MÞ;
ð1Þ
where l0 denotes the magnetic permeability of free space and B stands for the magnetic induction. The usual Maxwell equations are valid in X e0 ot E r H ¼ J; ot B þ r E ¼ 0:
ð2Þ
Here e0 is the permitivity of free space, E denotes the electric field and J is the current density. The conventional theory of ferromagnetic materials is based on the assumption, following [16,17], that the magnetization M varies with the position x, but it has a fixed temperature-dependent magnitude jMj ¼ M. We assume that temperature is below a critical value (CurieÕs temperature depending on the material). The evolution of M is governed by the Landau–Lifshitz–Gilbert (LLG) equation jcj M ot M ¼ ðH H ðH; MÞ M þ a ðH; MÞ MÞ T T 1 þ a2 jMj ¼ f ðH; MÞ:
ð3Þ
The symbol c represents the gyromagnetic factor and a is a damping constant. The total magnetic field H T ðH; MÞ in the ferromagnet is H T ðH; MÞ ¼ H þ H s þ KP ðMÞ;
ð4Þ
where H s is a given static field. The material is characterized by the constant K (in reality K may depend on position). We consider a ferromagnetic crystal with one distinguished (the easiest magnetization) axis, which is represented by a unit vector p, jpj ¼ 1. The projection of M on p is denoted by P ðMÞ P ðMÞ ¼ ðp MÞp:
ð5Þ
We neglect the exchange magnetic field in (4). This is possible in some situations, like modeling radar absorbing materials for stealth applications. For a more complete discussion of the model see, e.g., [1,8,9,23,24]. The first term in the right-hand side of (3) causes a precession of M around HT ðH; MÞ and the second term is dissipative. In fact, a scalar multiplication of (3) by M implies
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
81
1 2 ot M M ¼ ot jMj ¼ 0 2 from which we after the time integration deduce for any time t > 0 jMðtÞj ¼ jMð0Þj:
ð6Þ
Hence, the length of Mðt; xÞ remains constant during the time for any position x. Eliminating B in (2), we get the following Maxwell-LLG system e0 ot E r H ¼ J; l0 ot H þ r E ¼ l0 ot M; ot M ¼ f ðH; MÞ:
ð7Þ
For simplicity we assume that the boundary C of X is perfectly conducting so that mE ¼0
ð8Þ
on C;
where m stands for the outward unit normal vector on the boundary. 3 We pretend that the initial data Eð0Þ ¼ E 0 , Hð0Þ ¼ H 0 and Mð0Þ ¼ M 0 are sufficiently smooth for our purposes. For physical reasons we assume that r ðH 0 þ M 0 Þ ¼ 0
in X;
which implies that the magnetic induction B is divergence free. There exists a large body of literature devoted to the study of hysteresis. For more detailed references we refer the reader to [3,18,21,25]. Here are also described various models, e.g., Prandtl, Preisach, Duhamel. These books are dedicated to a mathematical discussion of the problem. The literature concerning the LLG model is also very rich. A single LLG equation has been intensively studied by many authors, e.g., [5–7,26,28]. The Maxwell–LLG system in 1D has been considered, for instance, in [10,13,14]; and the more dimensional situation has been studied in [11,12]. Monk and Vacus [20] proposed a family of mass-lumped finite element schemes for the Maxwell–LLG problem and they proved the existence of a class of Liapunov functions for the continuous problem. Monk and Vacus [15] studied the Maxwell–LLG problem. The authors showed how a certain class of finite element methods can be used to approximate the Maxwell–LLG equations while preserving energy decay and the modulus of magnetization M. Monk and Vacus [19] established the error estimates for the space discretization for a sufficiently smooth electromagnetic 3
One can also consider the following boundary conditions on C ¼ CE [ CH : m E ¼ 0 on CE and m H ¼ 0 on CH , where CE \ CH ¼ ;.
82
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
field. The time discretization was done by a finite difference scheme, but no error estimates have been derived. The main goal of this paper is to present two new approximation schemes for the time discretization of the Maxwell–LLG system. We use backward EulerÕs method for the time discretization of E, and H. We design two algorithms (linear and nonlinear) for the computation of the vector field M. The length of M is conserved in both cases. Supposing that the exact solution is sufficiently smooth, we derive the error estimates for the proposed numerical schemes, see Theorems 1 and 2. At the end we present a numerical example in order to demonstrate the efficiency of suggested algorithms.
2. Preliminaries Without loss of generality we put e0 ¼ l0 ¼ a ¼ K ¼ 1, c ¼ 2 and J ¼ H s ¼ 0 in the theoretical part of the paper, but not in the numerical one. We denote by (w; z) be the usual L2 -inner product of any real or pffiffiffiffiffiffiffiffiffiffiffiffi ffi vectorR valued functions w and z in X, i.e., ðwzÞ ¼ X w z and kwk ¼ ðw; wÞ. Further we set QT ¼ ð0; T Þ X and n o 3 3 Hðcurl; XÞ ¼ w 2 ðL2 ðXÞÞ ; r w 2 ðL2 ðXÞÞ ; H 0 ðcurl; XÞ ¼ f/ 2 Hðcurl; XÞ; m / ¼ 0 on Cg: Throughout the rest of the paper we assume the following regularity of an exact solution to the boundary value problem (7) and (8) 3
ot E 2 L1 ðð0; T Þ; ðL2 ðXÞÞ Þ; E 2 L1 ðð0; T Þ; H0 ðcurl; XÞÞ; ot H 2 L1 ðð0; T Þ; ðL2 ðXÞÞ3 Þ;
ð9Þ
H 2 L1 ðð0; T Þ; Hðcurl; XÞÞ \ L1 ðQT Þ; 3
ot M 2 L1 ðð0; T Þ; ðL2 ðXÞÞ Þ; M 2 L1 ðQT Þ: We rewrite (7) into the following variational form ðot E; nÞ ðr H; nÞ ¼ 0; ðot H; uÞ þ ðE; r uÞ ¼ ðot M; uÞ ðot M; wÞ ¼ ðf ðH; MÞ; wÞ
ð10Þ 3
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ and w 2 ðL2 ðXÞÞ . We denote by R a sufficiently large positive real number satisfying jHðt; xÞj 6 R
a:e: in QT :
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
83
Now, we introduce a function gR as 0 for s < 0 gR ðsÞ ¼ minðs; RÞ else and the vector function w wðHÞ ¼
gR ðjHjÞ H: jHj
ð11Þ
Using this cut-off function, we modify the right-hand side f by fR as follows: jcj M ½wðH T ðH; MÞÞ M : wðH T ðH; MÞÞ M þ a fR ðH; MÞ ¼ 1 þ a2 jMj ð12Þ We recall that the assumption H 2 L1 ðQT Þ together with (6) guarantee the Lipschitz continuity (in the L2 ðXÞ-norm) of the right-hand side f in (3)––see [19, Lemma 2.2], thus also the uniqueness of a solution ðE; H; MÞ of (10). Clearly, it also solves ðot E; nÞ ðr H; nÞ ¼ 0; ðot H; uÞ þ ðE; r uÞ ¼ ðot M; uÞ;
ð13Þ
ðot M; wÞ ¼ ðfR ðH; MÞ; wÞ 3
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ and w 2 ðL2 ðXÞÞ . Remark 1. C; e; Ce denote generic positive constants (i.e., they can differ from place to place) independent of the time step s (e is small and Ce ¼ Cðe1 Þ is large). We have modified the vector function f in order to ensure the global Lipschitz continuity for the right-hand side (see Lemma 1). Lemma 1. There exists a positive constant C depending on R such that jfR ðH 1 ; M 1 Þ fR ðH 2 ; M 2 Þj 6 C ðjH 1 H 2 j þ jM 1 M 2 jÞ holds for any H 1 ; M 1 ; H 2 ; M 2 obeying 0 < C0 6 jM 1 j; jM 2 j 6 C1 . Proof. We set | ¼ wðH 1 Þ M 1 wðH 2 Þ M 2 . First, we show the Lipschitz continuity of wðHÞ M, i.e., j|j 6 C ðjH 1 H 2 j þ jM 1 M 2 jÞ:
ð14Þ
84
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
In fact, we successively deduce gR ðjH 1 jÞ gR ðjH 2 jÞ gR ðjH 2 jÞ |¼ ðH 1 M 1 H 2 M 2 Þ H1 M1 þ jH 1 j jH 2 j jH 2 j gR ðjH 2 jÞ H 1 M1 ¼ gR ðjH 1 jÞ gR ðjH 2 jÞ þ ðjH 2 j jH 1 jÞ jH 2 j jH 1 j þ
gR ðjH 2 jÞ ½ðH 1 H 2 Þ M 1 þ H 2 ðM 1 M 2 Þ: jH 2 j
Using the triangle inequality jM 1 j; jM 2 j 6 C and the definition of the function gR , we have gR ðjH 2 jÞ j|j 6 C jgR ðjH 1 jÞ gR ðjH 2 jÞj þ kH 2 j jH 1 k jH 2 j gR ðjH 2 jÞ jH 1 H 2 j þ jgR ðjH 2 jÞjjM 1 M 2 j þ jH 2 j 6 C ðkH 2 j jH 1 k þ jH 1 H 2 j þ jM 1 M 2 jÞ 6 CðjH 1 H 2 j þ jM 1 M 2 jÞ; which proofs (14).Using the elementary inequality M1 M 2 1 1 6 max ; jM 1 M 2 j 6 CjM 1 M 2 j; jM j jM j jM j jM j 1
2
1
we can analogously show that also which concludes the proof. h
2
M jMj
ðwðHÞ MÞ is Lipschitz continuous,
3. Time discretization The time discretization is based on backward EulerÕs method. We use an equidistant partitioning with a time step s ¼ Tn , for any n 2 N. Therefore, we divide the time interval ½0; T into n subintervals ½ti1 ; ti for ti ¼ is. We introduce the following notation zi zi1 zi ¼ zðti Þ; dzi ¼ s for any function z. First, we suggest the following linear recurrent approximation scheme for i ¼ 1; . . . ; n: Algorithm 1 (linear) (1) We start from ei1 ; hi1 and mi1 taking into account e0 ¼ E 0 , h0 ¼ H 0 and m0 ¼ M 0 . (2) We solve the linear ordinary differential equation (ODE) with an unknown mðtÞ on the subinterval ½ti1 ; ti
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
ot m ¼ wðH T ðhi1 ; mi1 ÞÞ m þ
85
m ðwðH T ðhi1 ; mi1 ÞÞ mi1 Þ: jmi1 j ð15Þ
(3) We set mi :¼ mðti Þ. (4) We solve the system of partial differential equations (PDEs) with the unknowns ei ; hi ðdei ; nÞ ðr hi ; nÞ ¼ 0; ðdhi ; uÞ þ ðei ; r uÞ ¼ ðot mðti Þ; uÞ
ð16Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. Inspecting (15) we see that after a scalar multiplication by m, we obtain 1 ot m m ¼ ot jmj2 ¼ 0: 2 Further, we integrate this over ½ti1 ; t and we get jmi1 j ¼ jmðtÞj
for t 2 ½ti1 ; ti ;
i.e., the magnetization m conserves the modulus during the time for any fixed position x 2 X. 3.1. Existence of ei ; hi ; mi The well posedness of (15) for any x 2 X is guaranteed by the following mi1 lemma, where we set a ¼ wðH T ðhi1 ; mi1 ÞÞ wðH T ðhi1 ; mi1 ÞÞ jm . i1 j Lemma 2. Let a and u0 be any vectors in R3 . Then the solution of ot uðtÞ ¼ a uðtÞ
t > 0;
ð17Þ
uð0Þ ¼ u0 is given by k
uðtÞ ¼ u0 þ u? 0 cosðjajtÞ þ k
a u? 0 sinðjajtÞ; jaj
ð18Þ
k
? where u0 ¼ u? 0 þ u0 , u0 is parallel to a, and u0 is perpendicular to a. Moreover, the vector field uðtÞ preserves its modulus, i.e., juðtÞj ¼ ju0 j for any time t > 0.
Proof. The assertion can be readily obtained by differentiation of (18) with respect to the time variable. Therefore we skip the details. h Once we have determined ot mðti Þ 2 ðL2 ðXÞÞ3 for a given i, we have to prove the existence of ei and hi .
86
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99 3
Lemma 3. For any i ¼ 1; . . . ; n, there exist uniquely determined ei ; hi 2 ðL2 ðXÞÞ satisfying (16), where r hi and r ei exist in a distributional sense.
Proof. Suppose that we want to prove the existence and uniqueness of solution for the following system ðe; nÞ ðr h; nÞ ¼ ða; nÞ;
ð19Þ
ðh; uÞ þ ðe; r uÞ ¼ ðb; uÞ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ and some a; b 2 ðL2 ðXÞÞ3 . We approximate the spaces H0 ðcurl; XÞ and Hðcurl; XÞ by their finite dimensional subspaces H0 ðcurl; XÞk ¼ ½n1 ; . . . ; nk and Hðcurl; XÞl ¼ ½u1 ; . . . ; ul , respectively. We assume the following approximation property lim kn PH0 ðcurl;XÞk nkHðcurl;XÞ ¼ lim ku PHðcurl;XÞl ukHðcurl;XÞ ¼ 0
k!1
l!1
ð20Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. Here PH0 ðcurl;Xk ; PHðcurl;Xl denote the projectors on H0 ðcurl; XÞk ; Hðcurl; XÞl , respectively. The approximation of (19) reads as ðek ; PH0 ðcurl;XÞk nÞ ðr hl ; PH0 ðcurl;XÞk nÞ ¼ ða; PH0 ðcurl;XÞk nÞ; ðhl ; PHðcurl;XÞl uÞ þ ðek ; r PHðcurl;Xl uÞ ¼ ðb; PHðcurl;XÞl uÞ
ð21Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. Pk Pl We are looking for ek ; hl of the form ek ¼ i¼1 ei ni and hl ¼ i¼1 hi ui . Setting this into (21) we obtain the following linear algebraic system (where the superscript T denotes transpose) Aek Chl ¼ f ; CT ek þ Bhl ¼ g; along with
A ¼ ðni ; nj Þ i;j¼1;...;k ;
B ¼ ðui ; uj Þ i;j¼1;...;l ;
C ¼ ðr ui ; nj Þ i¼1;...;l;
T f ¼ ða; nj Þ j¼1;...;k ;
T g ¼ ðb; uj Þ j¼1;...;l :
j¼1;...;k
;
This form is possible because of the identity ðr ui ; nj Þ ¼ ðui ; r nj Þ; which is valid for any i ¼ 1; . . . ; l and j ¼ 1; . . . ; k.
ð22Þ
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
87
Matrices A and B are symmetric and regular. Eliminating ek from (22a) and setting the result into (22.b), we get ek ¼ A1 f þ A1 Chl ; B þ CT A1 C hl ¼ g CT A1 f : Due to the fact that the matrix B þ CT A1 C is regular, we have proved the existence of solution of (21). Moreover, this is also uniquely determined. As a next step we need uniform a priori estimates for ek and hl with respect to k; l. Therefore, we set n ¼ ek , u ¼ hl in (21) and sum up both equation. We get 2
2
kek k þ khl k ¼ ða; ek Þ þ ðb; hl Þ: Using the Cauchy and YoungÕs inequalities to the right-hand side, we readily obtain kek k2 þ khl k2 6 kak2 þ kbk2 ; which is valid for any k; l 2 N. Thus, due to the reflexivity of the space 3 ðL2 ðXÞÞ , we can choose such subsequences from fek g and fhl g (denoted by the same symbol again) for which ek * e;
hl * h
3
in ðL2 ðXÞÞ
for some e; h 2 ðL2 ðXÞÞ3 . In view of the approximation property (20), we deduce k!1
ðek ; PH0 ðcurl;XÞk nÞ ¼ ðek ; nÞ þ ðek ; PH0 ðcurl;XÞk n nÞ ! ðe; nÞ; ða; PH0 ðcurl;XÞk nÞ
k!1
! ða; nÞ;
and applying the Green formula we have l!1
ðr hl ; PHðcurl;XÞl nÞ ¼ ðhl ; r nÞ þ ðhl ; r ðPHðcurl;XÞl n nÞÞ ! h; r nÞ ¼ ðr h; nÞ:
Recall that r h exists in a distributional sense.Analogously we pass to the limit for k; l ! 1 in (21b). Hence, e and h solve (20). One can easily check that this solution is uniquely determined. h 3.2. A priori estimates Next step is to derive suitable a priori estimates for ei ; hi , i ¼ 1; . . . ; n. We do it in a few steps.
88
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
Lemma 4. There exist positive constants C and s0 such that j j X X kei ei1 k2 þ khi hi1 k2 6 C kej k2 þ khj k2 þ i¼1
i¼1
holds for all j ¼ 1; . . . ; n and all 0 < s < s0 . Proof. Examining (15), we can easily get jot mðtÞj 6 C for any t. Now, we set n ¼ ei and u ¼ hi in (16), then we sum both equations and get ðdei ; ei Þ þ ðdhi ; hi Þ ¼ ðot mðti Þ; hi Þ: We multiply this by s and sum it up for i ¼ 1; . . . ; j. Using Abel summation for the left-hand side and CauchyÕs and YoungÕs inequalities for the right-hand side, we deduce ! j j j X X X 2 2 2 2 2 kej k þ khj k þ kei ei1 k þ khi hi1 k 6 C 1 þ khi k s : i¼1
i¼1
i¼1
The rest of the proof follows from GronwallÕs argument.
h
For the next lemma we need the following compatibility condition ðot Eð0Þ; nÞ ðr H 0 ; nÞ ¼ 0; ðot Hð0Þ; uÞ þ ðE 0 ; r uÞ ¼ ðot Mð0Þ; uÞ;
ð23Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. This says that the Maxwell equations (2) are satisfied at the time t ¼ 0. We define de0 ¼ ot Eð0Þ;
dh0 ¼ ot Hð0Þ;
ot mð0Þ ¼ ot Mð0Þ:
Let us note that a simple calculation yields kot Mð0Þk 6 Cð1 þ kH 0 kÞ 6 C: Lemma 5. Suppose (23). Then there exists a positive constants C and s0 such that 2
2
kdej k þ kdhj k þ
j X
2
kdei dei1 k þ
i¼1
j X
2
kdhi dhi1 k 6 C
i¼1
holds for all j ¼ 1; . . . ; n and all 0 < s < s0 . Proof. We subtract (16) for i ¼ i 1 from (16), then we set n ¼ dei , u ¼ dhi and we sum up both equations. For i ¼ 0 we use the compatibility condition. We have
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
89
ðdei dei1 ; dei Þ þ ðdhi dhi1 ; dhi Þ ¼ ðot mðti1 Þ ot mðti Þ; dhi Þ: We sum up this equation for i ¼ 1; . . . ; j. Then we apply Abel summation, CauchyÕs and YoungÕs inequalities and we get in a straightforward way 2
2
kdej k þ kdhj k þ
j X
2
kdei dei1 k þ
i¼1
6C 1 þ
j X
j X
kdhi dhi1 k
2
i¼1
! 2
kdhi k s :
i¼1
The desired result follows from GronwallÕs lemma. h 3.3. Error estimates We introduce the piecewise linear in time vector fields en ; hn (i ¼ 1; . . . ; n) given by en ð0Þ ¼ E 0 ; en ðtÞ ¼ ei1 þ ðt ti1 Þdei
for t 2 ðti1 ; ti
and hn ð0Þ ¼ H 0 ; hn ðtÞ ¼ hi1 þ ðt ti1 Þdhi
for t 2 ðti1 ; ti :
Next, we define the step vector fields en ; hn ; mn and ot mn en ð0Þ ¼ E 0 ;
en ðtÞ ¼ ei ;
hn ð0Þ ¼ H 0 ;
hn ðtÞ ¼ hi ;
mn ð0Þ ¼ M 0 ;
mn ðtÞ ¼ mi ;
ot mn ð0Þ ¼ ot Mð0Þ;
ot mn ðtÞ ¼ ot mðti Þ;
for t 2 ðti1 ; ti :
Further, we prescribe the vector field mn as follows: mn ðtÞ ¼ mðtÞ
for t 2 ½ti1 ; ti
and for all i ¼ 1; . . . ; n. Using the new notation we rewrite (16) into a more suitable form for our purposes ðot en ; nÞ ðr hn ; nÞ ¼ 0; ðot hn ; uÞ þ ðen ; r uÞ ¼ ðot mn ; uÞ;
ð24Þ
which holds for any n 2 H0 ðcurl; XÞ and u 2 Hðcurl; XÞ. Now, we are in a position to derive the error estimates for the linear approximation scheme (15) and (16).
90
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
Theorem 1. Let the assumptions of Lemma 5 be satisfied. There exist positive constants C and s0 such that (i) maxt2½0;T kEðtÞ en ðtÞk2 þ maxt2½0;T kHðtÞ hn ðtÞk2 6 Cs; (ii) maxt2½0;T kMðtÞ mn ðtÞk2 þ maxt2½0;T kot M ot mn k2 6 Cs, hold for any 0 < s < s0 . Proof. (i) According to the definitions of the vector fields M and mn we can write for any time t ot MðtÞ ot mn ðtÞ ¼ wðH T ðHðtÞ; MðtÞÞÞ MðtÞ MðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ MðtÞ jMðtÞj
w H T hn ðt sÞ; mn ðt sÞ mn ðtÞ
mn ðtÞ w H T hn ðt sÞ; mn ðt sÞ mn ðt sÞ jmn ðtÞj þ
¼ R1 þ R 2 þ R 3 þ R 4 þ R5 ;
ð25Þ
where R1 ¼ wðH T ðHðtÞ; MðtÞÞÞ MðtÞ
w H T hn ðt sÞ; mn ðt sÞ MðtÞ
R2 ¼ w H T hn ðt sÞ; mn ðt sÞ ðMðtÞ mn ðtÞÞ MðtÞ mn ðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ MðtÞ R3 ¼ jMðtÞj mn ðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ ðMðtÞ mn ðtÞÞ R4 ¼ jMðtÞj mn ðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ mn ðtÞ R5 ¼ jMðtÞj
mn ðtÞ w H T hn ðt sÞ; mn ðt sÞ mn ðtÞ jMðtÞj
mn ðtÞ w H T hn ðt sÞ; mn ðt sÞ ðmn ðtÞ mn ðt sÞÞ : þ jMðtÞj Applying the Lipschitz continuity of w we get in a straightforward way
jot MðtÞ ot mn ðtÞj 6 C jHðtÞ hn ðt sÞj þ jMðtÞ mn ðtÞj þ CjMðtÞ mn ðt sÞj
6 C jHðtÞ hn ðtÞj þ jhn ðtÞ hn ðt sÞj þ C ðjMðtÞ mn ðtÞj þ jmn ðtÞ mn ðt sÞjÞ:
ð26Þ
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
91
The vector fields mn and M are continuous in time. Moreover, they start from the same initial datum M 0 . Hence, we can write Z t ot M ot m n : ð27Þ MðtÞ mn ðtÞ ¼ 0
According to (26) and Lemma 5 we have Z t kMðtÞ mn ðtÞk 6 kot M ot mn k 0 Z t 6C s þ ðkH hn k þ skot hn k þ kM mn kÞ 0 Z t 6C s þ ðkH hn k þ kM mn kÞ : 0
Here, we have used the inequality Z t kot mn k 6 Cs: kmn ðtÞ mn ðt sÞk 6 ts
GronwallÕs lemma gives Z t kMðtÞ mn ðtÞk 6 C s þ kH hn k :
ð28Þ
0
This relation together with (26) imply Z t Z t 2 2 2 kot M ot mn k 6 C s þ kH hn k : 0
ð29Þ
0
Now, we subtract (24a) from 13a and (24b) from (13b). We have ðot ðE en Þ; nÞ ðr H hn ; nÞ ¼ 0; ðot ðH hn Þ; uÞ þ ðE en ; r uÞ ¼ ðot mn ot M; uÞ;
ð30Þ
which holds for any n 2 H0 ðcurl; XÞ and u 2 Hðcurl; XÞ. Further we set n ¼ H hn , u ¼ E en , sum both equations up and integrate over time. We obtain 1 1 2 2 kEðtÞ en ðtÞk þ kHðtÞ hn ðtÞk 2 2 Z t Z t ¼ ðot ðE en Þ; en en Þ þ ðot ðH hn Þ; hn hn Þ 0 0 Z t þ ðot mn ot M; H hn Þ:
ð31Þ
0
For the first term on the right-hand side we deduce using the Cauchy inequality and Lemma 5
92
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z t sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z t Z t 2 2 ðot ðE en Þ; en en Þ 6 kot ðE en Þk ken en k 6 Cs: 0
0
0
ð32Þ The second term on the right in (31) can be estimated in the same way, i.e., Z t ðot ðH hn Þ; hn hn Þ 6 Cs: ð33Þ 0
For any t 2 ½ti1 ; ti we can write kot mn ðtÞ ot MðtÞk 6 kot mn ðti Þ ot Mðti Þk þ kot Mðti Þ ot MðtÞk:
ð34Þ
The difference ot Mðti Þ ot MðtÞ can be estimated using Lemma 1 as follows: kot Mðti Þ ot MðtÞk 6 C ðkMðti Þ MðtÞk þ kHðti Þ HðtÞkÞ 6 Cs: For the term ot mn ðti Þ ot Mðti Þ we apply (26), triangle inequality and (28) kot mn ðti Þ ot Mðti Þk 6 C ðkHðti Þ hn ðti1 Þk þ kMðti Þ mn ðti ÞkÞ þ CkMðti Þ mn ðti1 Þk 6 C ðs þ kHðti Þ hn ðti Þk þ kMðti Þ mn ðti ÞkÞ Z t 6 C s þ kH n ðtÞ hn ðtÞk þ kH hn k : 0
Taking the second power in (34) and integrating over the time, we deduce Z t Z t Z t 2 2 2 2 kot mn ot Mk 6 C s þ kH n hn k þ kH hn k 0 0 0 Z t 6 C s2 þ kH hn k2 : 0
For the third term on the right in (31) we deduce Z t Z t Z t 2 2 ðot mn ot M; H hn Þ 6 ko m o Mk þ kH hn k t n t 0 0 0 Z t 2 6 C s2 þ kH hn k : 0
Summarizing (31)–(33) and (35) we arrive at Z t 2 2 2 kH hn k : kHðtÞ hn ðtÞk þ kEðtÞ en ðtÞk 6 C s þ 0
Applying the Gronwall lemma we get kHðtÞ hn ðtÞk2 þ kEðtÞ en ðtÞk2 6 Cs:
ð35Þ
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
93
Due to the fact that this relation is valid for all t 2 ½0; T , we conclude the proof. (ii) The assertion immediately follows from Theorem 1 (i) and the relations (26) and (27) and GronwallÕs lemma. h
4. Nonlinear approximation scheme In this section we present the following recurrent nonlinear approximation scheme for i ¼ 1; . . . ; n: Algorithm 2 (nonlinear) (1) We start from ei1 ; hi1 and mi1 taking into account e0 ¼ E 0 , h0 ¼ H 0 and m0 ¼ M 0 . (2) We solve the nonlinear ODE with an unknown mðtÞ on the subinterval ½ti1 ; ti m ðwðH T ðhi1 ; mi1 ÞÞ mÞ: ot m ¼ wðH T ðhi1 ; mi1 ÞÞ m þ ð36Þ jmi1 j (3) We set mi :¼ mðti Þ. (4) We solve the system of PDEs with the unknowns ei ; hi ðdei ; nÞ ðr hi ; nÞ ¼ 0;
ð37Þ
ðdhi ; uÞ þ ðei ; r uÞ ¼ ðot mðti Þ; uÞ for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ.
Here, the approximation of M is nonlinear comparing with (15). Also this scheme conserves the modulus of m. This can be shown exactly in the same was as before. Let R be any rotation of the coordinate system. One can easily check that Rðr sÞ ¼ Rr Rs;
8r; s 2 R3 :
Thus, if u solves ot u ¼ a u þ cu ða uÞ t > 0; uð0Þ ¼ u0 ;
ð38Þ
then ot Ru ¼ Ra Ru þ cRu ðRa RuÞ Ruð0Þ ¼ Ru0 :
t > 0;
Therefore, without loss of generality we can assume that the vector a is parallel to the x-axis. This will help us to find the exact formula for the solution. This special choice of a can be obtained by the following suitable rotation.
94
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
Let us choose any x 2 X. The existence of m on ½ti1 ; ti , which solves (36), is proved in the following lemma. Let us set ~ u0 ¼ mi1 ;
~ a ¼ wðH T ðhi1 ; mi1 ÞÞ;
c¼
1 : jmj
Now, we perform a rotation R of the coordinate system in such a way that a :¼ R~ a ¼ j~ ajð1; 0; 0ÞT : At the end we denote u0 :¼ R~ u0 and a ¼ j~ aj. T
Lemma 6. Assume a; c 2 R and u0 ¼ ðx0 ; y0 ; z0 Þ 2 R3 . Then Eq. (38) for T T a ¼ að1; 0; 0Þ admits the solution uðtÞ ¼ ðxðtÞ; yðtÞ; zðtÞÞ , which is given by eactju0 j ðju0 j þ x0 Þ eactju0 j ðju0 j x0 Þ ; eactju0 j ðju0 j þ x0 Þ þ eactju0 j ðju0 j x0 Þ y0 cosðatÞ z0 sinðatÞ yðtÞ ¼ 2ju0 j actju j ; 0 e ðju0 j þ x0 Þ þ eactju0 j ðju0 j x0 Þ y0 sinðatÞ þ z0 cosðatÞ zðtÞ ¼ 2ju0 j actju j : e 0 ðju0 j þ x0 Þ þ eactju0 j ðju0 j x0 Þ xðtÞ ¼ ju0 j
Proof. See [22].
h
Further we follow the same way as in Section 3, i.e., we show the existence of ei ; hi for all i ¼ 1; . . . ; n and we prove the same a priori estimates as in Lemmas 4 and 5. Then we can derive the same convergence rates as we did in Theorem 1. We state the result without the proof. Theorem 2. Let the assumptions of Lemma 5 be satisfied. There exist positive constants C and s0 such that (i) maxt2½0;T kEðtÞ en ðtÞk2 þ maxt2½0;T kHðtÞ hn ðtÞk2 6 Cs, (ii) maxt2½0;T kMðtÞ mn ðtÞk2 þ maxt2½0;T kot M ot mn k2 6 Cs, hold for any 0 < s < s0 .
5. Numerical experiment In this section we show one example which demonstrates some of the features which make ferromagnetic layers interesting. We consider a rectangular parallelepipedon with the length Lx ¼ 16, the width Ly ¼ 23 and the height Lz ¼ 0:4 (see Fig. 1). This domain is occupied by a
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
95
Lz
z y x
Ly Lx Fig. 1. Computational domain X.
ferromagnet. We split the domain into small blocks with Dx ¼ Dy ¼ Dz ¼ 152 . Afterwards we divide each one into 6 tetrahedrons. For the approximation of the fields E and H we use WhitneyÕs edge elements, cf. [2,4]. Recall that the approximation of M can be settled down at any point x from the approximation of the LLG equation––see (36), which is used for computations. We have determined m at the barycenter of each tetrahedron, i.e., the approximation of M is a piecewise constant vector field. We consider the following boundary conditions
T left: m Eðt; x; y; zÞ ¼ m ð0; 0; 1Þ E0 g f t xc , right: m E ¼ 0, bottom: m E ¼ 0, top: m E ¼ 0, front: m H ¼ 0, back: m H ¼ 0, where E0 ¼ 4:0104 and g is the Ricker pulse (see Fig. 2) 8 < 2k ð2kðs 1:35Þ2 1Þ expðkðs 1:35Þ2 Þ if t < 3:48; gðsÞ ¼ 11 : 0 if t P 3:48 p 2 and k ¼ 1:31 . The values for f and c are given in Table 1. Let us note that we are using the following notation for the coordinates: xdashed, y-gray and z-solid. The initial data are given as follows: E 0 ¼ H 0 ¼ 0 T and M 0 ¼ ð0; 0; 100Þ , which ensures that r B ¼ 0. The computed results for the fields E; H and M along the horizontal center line of the domain are shown in Figs. 3–5, respectively. The waves correspond to two different time points t ¼ 6:705 108 s and t ¼ 1:3205 107 s. The incoming wave for E at the left boundary part has vanishing x- and y-components. When the wave enters the ferromagnetic layer, it is rotated, which makes x- and y-components nonzero. Moreover, due to the fact that a 6¼ 0 is the layer absorbing.
96
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
10000 0 –10000 –20000 –30000 –40000 –20
–15
–5
–10
0
s Fig. 2. Ricker pulse g.
Table 1 Parameters Parameter
Value
Parameter
Value
K e0 c s
0 0.88422 · 1011 Fm1 3.0 · 108 ms1 0.5 · 1010 s
a l0 c f
0.2 1.25667 · 106 Hm1 2.2 · 105 Hz mA1 4.0 · 107 Hz
40000
40000
20000
20000
0
0
–20000
–20000
–40000
–40000 2
4
6
8
10 12 14 16
2
4
6
8
10 12 14 16
Fig. 3. Left: Incoming E-wave at t ¼ 6:705 108 s. Right: Reflected E-wave at t ¼ 1:3205 107 s.
The electric field E in Fig. 3 has a dominant z-coordinate. The incident wave is traveling to the right and it is reflected and rotated by the right wall. Its amplitude is decreasing in time. The magnetic field H in Fig. 4 has a dominant
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99 100
100
50
50
0
0
–50
–50
–100
–100 2
4
6
8
Fig. 4. Left: Incoming t ¼ 1:3205 107 s.
2
10 12 14 16
H-wave
at
8
t ¼ 6:705 10
100
100
80
80
60
60
40
40
20
20
0
0 2
4
6
8
10 12 14 16
2
Fig. 5. Left: Incoming M-wave at t ¼ 6:705 108 t ¼ 1:3205 107 s.
4
s.
6
8
Right:
4
6
10 12 14 16
Reflected
8
97
H-wave
at
10 12 14 16
s. Right: Reflected M-wave at
y-coordinate, but here the reflected wave is not rotated. The amplitude is also decreasing. The magnetization M shows only small changes in Fig. 5. The initial datum was Mð0Þ ¼ ð0; 0; 100ÞT . Due to the fact that the magnitude of M remains constant, we can see that the increase of jxj; jyj-components implies the drop of the z-coordinate.
E
H H v
v E
Fig. 6. Position among E, H and the velocity vector v. Left: Incoming wave. Right: Reflected wave.
98
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
The mutual position among E; H and the velocity vector v of the wave is fixed. Therefore, when the reflected wave is coming back, then the wave for H keeps its orientation and E is ‘‘upside down’’ (see Fig. 6).
Acknowledgements The authors would like to express their thanks to Prof. R. Hiptmair (ETH, Zurich, Switzerland) and Dr. L. Dupre (Ghent University, Belgium) for many encouraging and enlightening discussions. References [1] G. Anzellotti, S. Baldo, A. Visintin, Asymptotic behavior of the Landau–Lifshitz model of ferromagnetism, Appl. Math. Optim. 23 (3) (1991) 171–192. [2] A. Bossavit, Computational Electromagnetism, Variational Formulations, Complementarity, Edge Elements, in: Electromagnetism, vol. XVIII, Academic Press, Orlando, FL, 1998. [3] M. Brokate, J. Sprekels, Hysteresis and Phase Transitions, in: Applied Mathematical Sciences, vol. 121, Springer, 1996. [4] M. Cessenat, Mathematical Methods in Electromagnetism, Linear Theory and Applications, in: Advances in Mathematics for Applied Sciences, vol. 41, World Scientific Publishers, Singapore, 1996. [5] Y. Chen, Dirichlet boundary value problems of Landau–Lifshitz equation, Commun. Partial Differ. Eqs. 25 (1–2) (2000) 101–124. [6] Y. Chen, Existence and singularities for the Dirichlet boundary value problems of Landau– Lifshitz equations, Nonlinear Anal., Theory Meth. Appl. 48A (3) (2002) 411–426. [7] W. E, X.-P. Wang, Numerical methods for the Landau–Lifshitz equation, SIAM J. Numer. Anal. 38 (5) (2000) 1647–1665. [8] B. Guo, S. Ding, Neumann problem for the Landau–Lifshitz–Maxwell system in two dimensions, Chin. Ann. Math. Ser. B 22 (4) (2001) 529–540. [9] B. Guo, F. Su, Global weak solution for the Landau–Lifshitz–Maxwell equation in three space dimensions, J. Math. Anal. Appl. 211 (1) (1997) 326–346. [10] H. Haddar, P. Joly, Effective boundary conditions for thin ferromagnetic layers: The onedimensional model, SIAM J. Appl. Math. 61 (4) (2001) 1386–1417. [11] J.L. Joly, G. Metivier, J. Rauch, Electromagnetic wave propagation in presence of a ferromagnetic material, ESAIM, Proc. 3 (1998) 85–99. [12] J.L. Joly, G. Metivier, J. Rauch, Global solutions to Maxwell equations in a ferromagnetic medium, Ann. Henri Poincare 1 (2) (2000) 307–340. [13] P. Joly, A. Komech, O. Vacus, On transitions to stationary states in a Maxwell–Landau– Lifschitz–Gilbert system, SIAM J. Math. Anal. 31 (2) (2000) 346–374. [14] P. Joly, O. Vacus, MaxwellÕs equations in a 1D ferromagnetic medium: existence and uniqueness of strong solutions, INRIA, technical report no. 3052, 1996. [15] P. Joly, O. Vacus, Mathematical and numerical studies of nonlinear ferromagnetic materials, M2AN 33 (3) (1999) 593–626. [16] L.D. Landau, E.M. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies, Phys. Z. Sowjetunion 8 (1935) 153–169. [17] L.D. Landau, E.M. Lifshitz, Electrodynamics of Continuous Media, Pergamon Press, Oxford, X, 1960, p. 417 (J.B. Sykes, J.S. Bell, Translated from the Russian).
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
99
[18] I. Mayergoyz, Nonlinear Diffusion of Electromagnetic Fields with Applications to Eddy Currents and Superconductivity, Academic Press, London, 1998. [19] P.B. Monk, O. Vacus, Error estimates for a numerical scheme for ferromagnetic problems, SIAM J. Numer. Anal. 36 (3) (1998) 696–718. [20] P.B. Monk, O. Vacus, Accurate discretization of a nonlinear micromagnetic problem, Comput. Meth. Appl. Mech. Eng. 190 (40–41) (2001) 5243–5269. [21] A. Prohl, Computational Micromagnetism, in: Advances in Numerical Mathematics, vol. xvi, Teubner, Leipzig, 2001. [22] M. Slodicka, I. Cimrak, Numerical study of nonlinear ferromagnetic materials, Appl. Numer. Math. 46 (1) (2003) 95–111. [23] F. Su, B. Guo, The global smooth solution for Landau–Lifshitz–Maxwell equation without dissipation, J. Partial Differ. Eqs. 11 (3) (1998) 193–208. [24] A. Visintin, On Landau–LifshitzÕ equations for ferromagnetism, Jpn. J. Appl. Math. 2 (1985) 69–84. [25] A. Visintin, Differential Models of Hysteresis, in: Applied Mathematical Sciences, vol. 111, Springer, 1994. [26] X.-P. Wang, C.J. Garcıa-Cervera, W. E, A Gauss–Seidel projection method for micromagnetics simulations, J. Comput. Phys. 171 (1) (2001) 357–372. [27] P. Weiss, LÕhypothese du champ moleculaire et la propriete ferromagnetique, J. Physique 6 (1907) 661–690. [28] J. Zhai, Existence and behavior of solutions to the Landau–Lifshitz equation, SIAM J. Math. Anal. 30 (4) (1999) 833–847.
A numerical scheme for a Maxwell–Landau–Lifshitz–Gilbert system Mari an Slodicka
*,1,
LÕubomır Ba nas
2
Department of Mathematical Analysis, Faculty of Engineering, Ghent University, Galglaan 2, B-9000 Ghent, Belgium
Abstract We consider MaxwellÕs equations together with a nonlinear dissipative magnetic law described by the Landau–Lifshitz–Gilbert equation. The paper is devoted to the numerical study of this problem. We introduce two algorithms for the time discretization. The modulus of magnetization is conserved in both cases. Assuming a sufficiently smooth electromagnetic field, we derive the error estimates for both approximation schemes. 2003 Elsevier Inc. All rights reserved. Keywords: Electromagnetic field; Maxwell–Landau–Lifshitz–Gilbert system; Time discretization; Hysteresis
1. Introduction Simulations of magnetic behavior based on the micromagnetic theory exhibits hysteresis. These systems have a highly nonlinear character. The mean field theory of ferromagnetism was proposed by Weiss [27], where the hysteresis was based on a positive feedback effect. Weiss assumed that ferromagnetic materials have a spontaneous magnetization M, even in the absence of any
*
Corresponding author. E-mail addresses: [email protected] (M. Slodicka), [email protected] (L. Ba nas). 1 Supported by the BOF/GOA-project no. 12 052 499 of Ghent University. 2 Supported by the project no. 174 IW 300 of Ghent University. 0096-3003/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.08.065
80
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
applied magnetic field H. The fact that they can look as demagnetized on a macromagnetic scale was assign to the presence of a microscopic domain structure. This issue was later described by Landau and Lifshitz [16,17]. Let a ferromagnetic material occupy a bounded domain X R3 with a Lipschitz boundary C. The time-dependent magnetic variables are related as follows: B ¼ l0 ðH þ MÞ;
ð1Þ
where l0 denotes the magnetic permeability of free space and B stands for the magnetic induction. The usual Maxwell equations are valid in X e0 ot E r H ¼ J; ot B þ r E ¼ 0:
ð2Þ
Here e0 is the permitivity of free space, E denotes the electric field and J is the current density. The conventional theory of ferromagnetic materials is based on the assumption, following [16,17], that the magnetization M varies with the position x, but it has a fixed temperature-dependent magnitude jMj ¼ M. We assume that temperature is below a critical value (CurieÕs temperature depending on the material). The evolution of M is governed by the Landau–Lifshitz–Gilbert (LLG) equation jcj M ot M ¼ ðH H ðH; MÞ M þ a ðH; MÞ MÞ T T 1 þ a2 jMj ¼ f ðH; MÞ:
ð3Þ
The symbol c represents the gyromagnetic factor and a is a damping constant. The total magnetic field H T ðH; MÞ in the ferromagnet is H T ðH; MÞ ¼ H þ H s þ KP ðMÞ;
ð4Þ
where H s is a given static field. The material is characterized by the constant K (in reality K may depend on position). We consider a ferromagnetic crystal with one distinguished (the easiest magnetization) axis, which is represented by a unit vector p, jpj ¼ 1. The projection of M on p is denoted by P ðMÞ P ðMÞ ¼ ðp MÞp:
ð5Þ
We neglect the exchange magnetic field in (4). This is possible in some situations, like modeling radar absorbing materials for stealth applications. For a more complete discussion of the model see, e.g., [1,8,9,23,24]. The first term in the right-hand side of (3) causes a precession of M around HT ðH; MÞ and the second term is dissipative. In fact, a scalar multiplication of (3) by M implies
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
81
1 2 ot M M ¼ ot jMj ¼ 0 2 from which we after the time integration deduce for any time t > 0 jMðtÞj ¼ jMð0Þj:
ð6Þ
Hence, the length of Mðt; xÞ remains constant during the time for any position x. Eliminating B in (2), we get the following Maxwell-LLG system e0 ot E r H ¼ J; l0 ot H þ r E ¼ l0 ot M; ot M ¼ f ðH; MÞ:
ð7Þ
For simplicity we assume that the boundary C of X is perfectly conducting so that mE ¼0
ð8Þ
on C;
where m stands for the outward unit normal vector on the boundary. 3 We pretend that the initial data Eð0Þ ¼ E 0 , Hð0Þ ¼ H 0 and Mð0Þ ¼ M 0 are sufficiently smooth for our purposes. For physical reasons we assume that r ðH 0 þ M 0 Þ ¼ 0
in X;
which implies that the magnetic induction B is divergence free. There exists a large body of literature devoted to the study of hysteresis. For more detailed references we refer the reader to [3,18,21,25]. Here are also described various models, e.g., Prandtl, Preisach, Duhamel. These books are dedicated to a mathematical discussion of the problem. The literature concerning the LLG model is also very rich. A single LLG equation has been intensively studied by many authors, e.g., [5–7,26,28]. The Maxwell–LLG system in 1D has been considered, for instance, in [10,13,14]; and the more dimensional situation has been studied in [11,12]. Monk and Vacus [20] proposed a family of mass-lumped finite element schemes for the Maxwell–LLG problem and they proved the existence of a class of Liapunov functions for the continuous problem. Monk and Vacus [15] studied the Maxwell–LLG problem. The authors showed how a certain class of finite element methods can be used to approximate the Maxwell–LLG equations while preserving energy decay and the modulus of magnetization M. Monk and Vacus [19] established the error estimates for the space discretization for a sufficiently smooth electromagnetic 3
One can also consider the following boundary conditions on C ¼ CE [ CH : m E ¼ 0 on CE and m H ¼ 0 on CH , where CE \ CH ¼ ;.
82
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
field. The time discretization was done by a finite difference scheme, but no error estimates have been derived. The main goal of this paper is to present two new approximation schemes for the time discretization of the Maxwell–LLG system. We use backward EulerÕs method for the time discretization of E, and H. We design two algorithms (linear and nonlinear) for the computation of the vector field M. The length of M is conserved in both cases. Supposing that the exact solution is sufficiently smooth, we derive the error estimates for the proposed numerical schemes, see Theorems 1 and 2. At the end we present a numerical example in order to demonstrate the efficiency of suggested algorithms.
2. Preliminaries Without loss of generality we put e0 ¼ l0 ¼ a ¼ K ¼ 1, c ¼ 2 and J ¼ H s ¼ 0 in the theoretical part of the paper, but not in the numerical one. We denote by (w; z) be the usual L2 -inner product of any real or pffiffiffiffiffiffiffiffiffiffiffiffi ffi vectorR valued functions w and z in X, i.e., ðwzÞ ¼ X w z and kwk ¼ ðw; wÞ. Further we set QT ¼ ð0; T Þ X and n o 3 3 Hðcurl; XÞ ¼ w 2 ðL2 ðXÞÞ ; r w 2 ðL2 ðXÞÞ ; H 0 ðcurl; XÞ ¼ f/ 2 Hðcurl; XÞ; m / ¼ 0 on Cg: Throughout the rest of the paper we assume the following regularity of an exact solution to the boundary value problem (7) and (8) 3
ot E 2 L1 ðð0; T Þ; ðL2 ðXÞÞ Þ; E 2 L1 ðð0; T Þ; H0 ðcurl; XÞÞ; ot H 2 L1 ðð0; T Þ; ðL2 ðXÞÞ3 Þ;
ð9Þ
H 2 L1 ðð0; T Þ; Hðcurl; XÞÞ \ L1 ðQT Þ; 3
ot M 2 L1 ðð0; T Þ; ðL2 ðXÞÞ Þ; M 2 L1 ðQT Þ: We rewrite (7) into the following variational form ðot E; nÞ ðr H; nÞ ¼ 0; ðot H; uÞ þ ðE; r uÞ ¼ ðot M; uÞ ðot M; wÞ ¼ ðf ðH; MÞ; wÞ
ð10Þ 3
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ and w 2 ðL2 ðXÞÞ . We denote by R a sufficiently large positive real number satisfying jHðt; xÞj 6 R
a:e: in QT :
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
83
Now, we introduce a function gR as 0 for s < 0 gR ðsÞ ¼ minðs; RÞ else and the vector function w wðHÞ ¼
gR ðjHjÞ H: jHj
ð11Þ
Using this cut-off function, we modify the right-hand side f by fR as follows: jcj M ½wðH T ðH; MÞÞ M : wðH T ðH; MÞÞ M þ a fR ðH; MÞ ¼ 1 þ a2 jMj ð12Þ We recall that the assumption H 2 L1 ðQT Þ together with (6) guarantee the Lipschitz continuity (in the L2 ðXÞ-norm) of the right-hand side f in (3)––see [19, Lemma 2.2], thus also the uniqueness of a solution ðE; H; MÞ of (10). Clearly, it also solves ðot E; nÞ ðr H; nÞ ¼ 0; ðot H; uÞ þ ðE; r uÞ ¼ ðot M; uÞ;
ð13Þ
ðot M; wÞ ¼ ðfR ðH; MÞ; wÞ 3
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ and w 2 ðL2 ðXÞÞ . Remark 1. C; e; Ce denote generic positive constants (i.e., they can differ from place to place) independent of the time step s (e is small and Ce ¼ Cðe1 Þ is large). We have modified the vector function f in order to ensure the global Lipschitz continuity for the right-hand side (see Lemma 1). Lemma 1. There exists a positive constant C depending on R such that jfR ðH 1 ; M 1 Þ fR ðH 2 ; M 2 Þj 6 C ðjH 1 H 2 j þ jM 1 M 2 jÞ holds for any H 1 ; M 1 ; H 2 ; M 2 obeying 0 < C0 6 jM 1 j; jM 2 j 6 C1 . Proof. We set | ¼ wðH 1 Þ M 1 wðH 2 Þ M 2 . First, we show the Lipschitz continuity of wðHÞ M, i.e., j|j 6 C ðjH 1 H 2 j þ jM 1 M 2 jÞ:
ð14Þ
84
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
In fact, we successively deduce gR ðjH 1 jÞ gR ðjH 2 jÞ gR ðjH 2 jÞ |¼ ðH 1 M 1 H 2 M 2 Þ H1 M1 þ jH 1 j jH 2 j jH 2 j gR ðjH 2 jÞ H 1 M1 ¼ gR ðjH 1 jÞ gR ðjH 2 jÞ þ ðjH 2 j jH 1 jÞ jH 2 j jH 1 j þ
gR ðjH 2 jÞ ½ðH 1 H 2 Þ M 1 þ H 2 ðM 1 M 2 Þ: jH 2 j
Using the triangle inequality jM 1 j; jM 2 j 6 C and the definition of the function gR , we have gR ðjH 2 jÞ j|j 6 C jgR ðjH 1 jÞ gR ðjH 2 jÞj þ kH 2 j jH 1 k jH 2 j gR ðjH 2 jÞ jH 1 H 2 j þ jgR ðjH 2 jÞjjM 1 M 2 j þ jH 2 j 6 C ðkH 2 j jH 1 k þ jH 1 H 2 j þ jM 1 M 2 jÞ 6 CðjH 1 H 2 j þ jM 1 M 2 jÞ; which proofs (14).Using the elementary inequality M1 M 2 1 1 6 max ; jM 1 M 2 j 6 CjM 1 M 2 j; jM j jM j jM j jM j 1
2
1
we can analogously show that also which concludes the proof. h
2
M jMj
ðwðHÞ MÞ is Lipschitz continuous,
3. Time discretization The time discretization is based on backward EulerÕs method. We use an equidistant partitioning with a time step s ¼ Tn , for any n 2 N. Therefore, we divide the time interval ½0; T into n subintervals ½ti1 ; ti for ti ¼ is. We introduce the following notation zi zi1 zi ¼ zðti Þ; dzi ¼ s for any function z. First, we suggest the following linear recurrent approximation scheme for i ¼ 1; . . . ; n: Algorithm 1 (linear) (1) We start from ei1 ; hi1 and mi1 taking into account e0 ¼ E 0 , h0 ¼ H 0 and m0 ¼ M 0 . (2) We solve the linear ordinary differential equation (ODE) with an unknown mðtÞ on the subinterval ½ti1 ; ti
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
ot m ¼ wðH T ðhi1 ; mi1 ÞÞ m þ
85
m ðwðH T ðhi1 ; mi1 ÞÞ mi1 Þ: jmi1 j ð15Þ
(3) We set mi :¼ mðti Þ. (4) We solve the system of partial differential equations (PDEs) with the unknowns ei ; hi ðdei ; nÞ ðr hi ; nÞ ¼ 0; ðdhi ; uÞ þ ðei ; r uÞ ¼ ðot mðti Þ; uÞ
ð16Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. Inspecting (15) we see that after a scalar multiplication by m, we obtain 1 ot m m ¼ ot jmj2 ¼ 0: 2 Further, we integrate this over ½ti1 ; t and we get jmi1 j ¼ jmðtÞj
for t 2 ½ti1 ; ti ;
i.e., the magnetization m conserves the modulus during the time for any fixed position x 2 X. 3.1. Existence of ei ; hi ; mi The well posedness of (15) for any x 2 X is guaranteed by the following mi1 lemma, where we set a ¼ wðH T ðhi1 ; mi1 ÞÞ wðH T ðhi1 ; mi1 ÞÞ jm . i1 j Lemma 2. Let a and u0 be any vectors in R3 . Then the solution of ot uðtÞ ¼ a uðtÞ
t > 0;
ð17Þ
uð0Þ ¼ u0 is given by k
uðtÞ ¼ u0 þ u? 0 cosðjajtÞ þ k
a u? 0 sinðjajtÞ; jaj
ð18Þ
k
? where u0 ¼ u? 0 þ u0 , u0 is parallel to a, and u0 is perpendicular to a. Moreover, the vector field uðtÞ preserves its modulus, i.e., juðtÞj ¼ ju0 j for any time t > 0.
Proof. The assertion can be readily obtained by differentiation of (18) with respect to the time variable. Therefore we skip the details. h Once we have determined ot mðti Þ 2 ðL2 ðXÞÞ3 for a given i, we have to prove the existence of ei and hi .
86
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99 3
Lemma 3. For any i ¼ 1; . . . ; n, there exist uniquely determined ei ; hi 2 ðL2 ðXÞÞ satisfying (16), where r hi and r ei exist in a distributional sense.
Proof. Suppose that we want to prove the existence and uniqueness of solution for the following system ðe; nÞ ðr h; nÞ ¼ ða; nÞ;
ð19Þ
ðh; uÞ þ ðe; r uÞ ¼ ðb; uÞ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ and some a; b 2 ðL2 ðXÞÞ3 . We approximate the spaces H0 ðcurl; XÞ and Hðcurl; XÞ by their finite dimensional subspaces H0 ðcurl; XÞk ¼ ½n1 ; . . . ; nk and Hðcurl; XÞl ¼ ½u1 ; . . . ; ul , respectively. We assume the following approximation property lim kn PH0 ðcurl;XÞk nkHðcurl;XÞ ¼ lim ku PHðcurl;XÞl ukHðcurl;XÞ ¼ 0
k!1
l!1
ð20Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. Here PH0 ðcurl;Xk ; PHðcurl;Xl denote the projectors on H0 ðcurl; XÞk ; Hðcurl; XÞl , respectively. The approximation of (19) reads as ðek ; PH0 ðcurl;XÞk nÞ ðr hl ; PH0 ðcurl;XÞk nÞ ¼ ða; PH0 ðcurl;XÞk nÞ; ðhl ; PHðcurl;XÞl uÞ þ ðek ; r PHðcurl;Xl uÞ ¼ ðb; PHðcurl;XÞl uÞ
ð21Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. Pk Pl We are looking for ek ; hl of the form ek ¼ i¼1 ei ni and hl ¼ i¼1 hi ui . Setting this into (21) we obtain the following linear algebraic system (where the superscript T denotes transpose) Aek Chl ¼ f ; CT ek þ Bhl ¼ g; along with
A ¼ ðni ; nj Þ i;j¼1;...;k ;
B ¼ ðui ; uj Þ i;j¼1;...;l ;
C ¼ ðr ui ; nj Þ i¼1;...;l;
T f ¼ ða; nj Þ j¼1;...;k ;
T g ¼ ðb; uj Þ j¼1;...;l :
j¼1;...;k
;
This form is possible because of the identity ðr ui ; nj Þ ¼ ðui ; r nj Þ; which is valid for any i ¼ 1; . . . ; l and j ¼ 1; . . . ; k.
ð22Þ
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
87
Matrices A and B are symmetric and regular. Eliminating ek from (22a) and setting the result into (22.b), we get ek ¼ A1 f þ A1 Chl ; B þ CT A1 C hl ¼ g CT A1 f : Due to the fact that the matrix B þ CT A1 C is regular, we have proved the existence of solution of (21). Moreover, this is also uniquely determined. As a next step we need uniform a priori estimates for ek and hl with respect to k; l. Therefore, we set n ¼ ek , u ¼ hl in (21) and sum up both equation. We get 2
2
kek k þ khl k ¼ ða; ek Þ þ ðb; hl Þ: Using the Cauchy and YoungÕs inequalities to the right-hand side, we readily obtain kek k2 þ khl k2 6 kak2 þ kbk2 ; which is valid for any k; l 2 N. Thus, due to the reflexivity of the space 3 ðL2 ðXÞÞ , we can choose such subsequences from fek g and fhl g (denoted by the same symbol again) for which ek * e;
hl * h
3
in ðL2 ðXÞÞ
for some e; h 2 ðL2 ðXÞÞ3 . In view of the approximation property (20), we deduce k!1
ðek ; PH0 ðcurl;XÞk nÞ ¼ ðek ; nÞ þ ðek ; PH0 ðcurl;XÞk n nÞ ! ðe; nÞ; ða; PH0 ðcurl;XÞk nÞ
k!1
! ða; nÞ;
and applying the Green formula we have l!1
ðr hl ; PHðcurl;XÞl nÞ ¼ ðhl ; r nÞ þ ðhl ; r ðPHðcurl;XÞl n nÞÞ ! h; r nÞ ¼ ðr h; nÞ:
Recall that r h exists in a distributional sense.Analogously we pass to the limit for k; l ! 1 in (21b). Hence, e and h solve (20). One can easily check that this solution is uniquely determined. h 3.2. A priori estimates Next step is to derive suitable a priori estimates for ei ; hi , i ¼ 1; . . . ; n. We do it in a few steps.
88
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
Lemma 4. There exist positive constants C and s0 such that j j X X kei ei1 k2 þ khi hi1 k2 6 C kej k2 þ khj k2 þ i¼1
i¼1
holds for all j ¼ 1; . . . ; n and all 0 < s < s0 . Proof. Examining (15), we can easily get jot mðtÞj 6 C for any t. Now, we set n ¼ ei and u ¼ hi in (16), then we sum both equations and get ðdei ; ei Þ þ ðdhi ; hi Þ ¼ ðot mðti Þ; hi Þ: We multiply this by s and sum it up for i ¼ 1; . . . ; j. Using Abel summation for the left-hand side and CauchyÕs and YoungÕs inequalities for the right-hand side, we deduce ! j j j X X X 2 2 2 2 2 kej k þ khj k þ kei ei1 k þ khi hi1 k 6 C 1 þ khi k s : i¼1
i¼1
i¼1
The rest of the proof follows from GronwallÕs argument.
h
For the next lemma we need the following compatibility condition ðot Eð0Þ; nÞ ðr H 0 ; nÞ ¼ 0; ðot Hð0Þ; uÞ þ ðE 0 ; r uÞ ¼ ðot Mð0Þ; uÞ;
ð23Þ
for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ. This says that the Maxwell equations (2) are satisfied at the time t ¼ 0. We define de0 ¼ ot Eð0Þ;
dh0 ¼ ot Hð0Þ;
ot mð0Þ ¼ ot Mð0Þ:
Let us note that a simple calculation yields kot Mð0Þk 6 Cð1 þ kH 0 kÞ 6 C: Lemma 5. Suppose (23). Then there exists a positive constants C and s0 such that 2
2
kdej k þ kdhj k þ
j X
2
kdei dei1 k þ
i¼1
j X
2
kdhi dhi1 k 6 C
i¼1
holds for all j ¼ 1; . . . ; n and all 0 < s < s0 . Proof. We subtract (16) for i ¼ i 1 from (16), then we set n ¼ dei , u ¼ dhi and we sum up both equations. For i ¼ 0 we use the compatibility condition. We have
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
89
ðdei dei1 ; dei Þ þ ðdhi dhi1 ; dhi Þ ¼ ðot mðti1 Þ ot mðti Þ; dhi Þ: We sum up this equation for i ¼ 1; . . . ; j. Then we apply Abel summation, CauchyÕs and YoungÕs inequalities and we get in a straightforward way 2
2
kdej k þ kdhj k þ
j X
2
kdei dei1 k þ
i¼1
6C 1 þ
j X
j X
kdhi dhi1 k
2
i¼1
! 2
kdhi k s :
i¼1
The desired result follows from GronwallÕs lemma. h 3.3. Error estimates We introduce the piecewise linear in time vector fields en ; hn (i ¼ 1; . . . ; n) given by en ð0Þ ¼ E 0 ; en ðtÞ ¼ ei1 þ ðt ti1 Þdei
for t 2 ðti1 ; ti
and hn ð0Þ ¼ H 0 ; hn ðtÞ ¼ hi1 þ ðt ti1 Þdhi
for t 2 ðti1 ; ti :
Next, we define the step vector fields en ; hn ; mn and ot mn en ð0Þ ¼ E 0 ;
en ðtÞ ¼ ei ;
hn ð0Þ ¼ H 0 ;
hn ðtÞ ¼ hi ;
mn ð0Þ ¼ M 0 ;
mn ðtÞ ¼ mi ;
ot mn ð0Þ ¼ ot Mð0Þ;
ot mn ðtÞ ¼ ot mðti Þ;
for t 2 ðti1 ; ti :
Further, we prescribe the vector field mn as follows: mn ðtÞ ¼ mðtÞ
for t 2 ½ti1 ; ti
and for all i ¼ 1; . . . ; n. Using the new notation we rewrite (16) into a more suitable form for our purposes ðot en ; nÞ ðr hn ; nÞ ¼ 0; ðot hn ; uÞ þ ðen ; r uÞ ¼ ðot mn ; uÞ;
ð24Þ
which holds for any n 2 H0 ðcurl; XÞ and u 2 Hðcurl; XÞ. Now, we are in a position to derive the error estimates for the linear approximation scheme (15) and (16).
90
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
Theorem 1. Let the assumptions of Lemma 5 be satisfied. There exist positive constants C and s0 such that (i) maxt2½0;T kEðtÞ en ðtÞk2 þ maxt2½0;T kHðtÞ hn ðtÞk2 6 Cs; (ii) maxt2½0;T kMðtÞ mn ðtÞk2 þ maxt2½0;T kot M ot mn k2 6 Cs, hold for any 0 < s < s0 . Proof. (i) According to the definitions of the vector fields M and mn we can write for any time t ot MðtÞ ot mn ðtÞ ¼ wðH T ðHðtÞ; MðtÞÞÞ MðtÞ MðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ MðtÞ jMðtÞj
w H T hn ðt sÞ; mn ðt sÞ mn ðtÞ
mn ðtÞ w H T hn ðt sÞ; mn ðt sÞ mn ðt sÞ jmn ðtÞj þ
¼ R1 þ R 2 þ R 3 þ R 4 þ R5 ;
ð25Þ
where R1 ¼ wðH T ðHðtÞ; MðtÞÞÞ MðtÞ
w H T hn ðt sÞ; mn ðt sÞ MðtÞ
R2 ¼ w H T hn ðt sÞ; mn ðt sÞ ðMðtÞ mn ðtÞÞ MðtÞ mn ðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ MðtÞ R3 ¼ jMðtÞj mn ðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ ðMðtÞ mn ðtÞÞ R4 ¼ jMðtÞj mn ðtÞ ½wðH T ðHðtÞ; MðtÞÞÞ mn ðtÞ R5 ¼ jMðtÞj
mn ðtÞ w H T hn ðt sÞ; mn ðt sÞ mn ðtÞ jMðtÞj
mn ðtÞ w H T hn ðt sÞ; mn ðt sÞ ðmn ðtÞ mn ðt sÞÞ : þ jMðtÞj Applying the Lipschitz continuity of w we get in a straightforward way
jot MðtÞ ot mn ðtÞj 6 C jHðtÞ hn ðt sÞj þ jMðtÞ mn ðtÞj þ CjMðtÞ mn ðt sÞj
6 C jHðtÞ hn ðtÞj þ jhn ðtÞ hn ðt sÞj þ C ðjMðtÞ mn ðtÞj þ jmn ðtÞ mn ðt sÞjÞ:
ð26Þ
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
91
The vector fields mn and M are continuous in time. Moreover, they start from the same initial datum M 0 . Hence, we can write Z t ot M ot m n : ð27Þ MðtÞ mn ðtÞ ¼ 0
According to (26) and Lemma 5 we have Z t kMðtÞ mn ðtÞk 6 kot M ot mn k 0 Z t 6C s þ ðkH hn k þ skot hn k þ kM mn kÞ 0 Z t 6C s þ ðkH hn k þ kM mn kÞ : 0
Here, we have used the inequality Z t kot mn k 6 Cs: kmn ðtÞ mn ðt sÞk 6 ts
GronwallÕs lemma gives Z t kMðtÞ mn ðtÞk 6 C s þ kH hn k :
ð28Þ
0
This relation together with (26) imply Z t Z t 2 2 2 kot M ot mn k 6 C s þ kH hn k : 0
ð29Þ
0
Now, we subtract (24a) from 13a and (24b) from (13b). We have ðot ðE en Þ; nÞ ðr H hn ; nÞ ¼ 0; ðot ðH hn Þ; uÞ þ ðE en ; r uÞ ¼ ðot mn ot M; uÞ;
ð30Þ
which holds for any n 2 H0 ðcurl; XÞ and u 2 Hðcurl; XÞ. Further we set n ¼ H hn , u ¼ E en , sum both equations up and integrate over time. We obtain 1 1 2 2 kEðtÞ en ðtÞk þ kHðtÞ hn ðtÞk 2 2 Z t Z t ¼ ðot ðE en Þ; en en Þ þ ðot ðH hn Þ; hn hn Þ 0 0 Z t þ ðot mn ot M; H hn Þ:
ð31Þ
0
For the first term on the right-hand side we deduce using the Cauchy inequality and Lemma 5
92
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z t sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z t Z t 2 2 ðot ðE en Þ; en en Þ 6 kot ðE en Þk ken en k 6 Cs: 0
0
0
ð32Þ The second term on the right in (31) can be estimated in the same way, i.e., Z t ðot ðH hn Þ; hn hn Þ 6 Cs: ð33Þ 0
For any t 2 ½ti1 ; ti we can write kot mn ðtÞ ot MðtÞk 6 kot mn ðti Þ ot Mðti Þk þ kot Mðti Þ ot MðtÞk:
ð34Þ
The difference ot Mðti Þ ot MðtÞ can be estimated using Lemma 1 as follows: kot Mðti Þ ot MðtÞk 6 C ðkMðti Þ MðtÞk þ kHðti Þ HðtÞkÞ 6 Cs: For the term ot mn ðti Þ ot Mðti Þ we apply (26), triangle inequality and (28) kot mn ðti Þ ot Mðti Þk 6 C ðkHðti Þ hn ðti1 Þk þ kMðti Þ mn ðti ÞkÞ þ CkMðti Þ mn ðti1 Þk 6 C ðs þ kHðti Þ hn ðti Þk þ kMðti Þ mn ðti ÞkÞ Z t 6 C s þ kH n ðtÞ hn ðtÞk þ kH hn k : 0
Taking the second power in (34) and integrating over the time, we deduce Z t Z t Z t 2 2 2 2 kot mn ot Mk 6 C s þ kH n hn k þ kH hn k 0 0 0 Z t 6 C s2 þ kH hn k2 : 0
For the third term on the right in (31) we deduce Z t Z t Z t 2 2 ðot mn ot M; H hn Þ 6 ko m o Mk þ kH hn k t n t 0 0 0 Z t 2 6 C s2 þ kH hn k : 0
Summarizing (31)–(33) and (35) we arrive at Z t 2 2 2 kH hn k : kHðtÞ hn ðtÞk þ kEðtÞ en ðtÞk 6 C s þ 0
Applying the Gronwall lemma we get kHðtÞ hn ðtÞk2 þ kEðtÞ en ðtÞk2 6 Cs:
ð35Þ
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
93
Due to the fact that this relation is valid for all t 2 ½0; T , we conclude the proof. (ii) The assertion immediately follows from Theorem 1 (i) and the relations (26) and (27) and GronwallÕs lemma. h
4. Nonlinear approximation scheme In this section we present the following recurrent nonlinear approximation scheme for i ¼ 1; . . . ; n: Algorithm 2 (nonlinear) (1) We start from ei1 ; hi1 and mi1 taking into account e0 ¼ E 0 , h0 ¼ H 0 and m0 ¼ M 0 . (2) We solve the nonlinear ODE with an unknown mðtÞ on the subinterval ½ti1 ; ti m ðwðH T ðhi1 ; mi1 ÞÞ mÞ: ot m ¼ wðH T ðhi1 ; mi1 ÞÞ m þ ð36Þ jmi1 j (3) We set mi :¼ mðti Þ. (4) We solve the system of PDEs with the unknowns ei ; hi ðdei ; nÞ ðr hi ; nÞ ¼ 0;
ð37Þ
ðdhi ; uÞ þ ðei ; r uÞ ¼ ðot mðti Þ; uÞ for any n 2 H0 ðcurl; XÞ, u 2 Hðcurl; XÞ.
Here, the approximation of M is nonlinear comparing with (15). Also this scheme conserves the modulus of m. This can be shown exactly in the same was as before. Let R be any rotation of the coordinate system. One can easily check that Rðr sÞ ¼ Rr Rs;
8r; s 2 R3 :
Thus, if u solves ot u ¼ a u þ cu ða uÞ t > 0; uð0Þ ¼ u0 ;
ð38Þ
then ot Ru ¼ Ra Ru þ cRu ðRa RuÞ Ruð0Þ ¼ Ru0 :
t > 0;
Therefore, without loss of generality we can assume that the vector a is parallel to the x-axis. This will help us to find the exact formula for the solution. This special choice of a can be obtained by the following suitable rotation.
94
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
Let us choose any x 2 X. The existence of m on ½ti1 ; ti , which solves (36), is proved in the following lemma. Let us set ~ u0 ¼ mi1 ;
~ a ¼ wðH T ðhi1 ; mi1 ÞÞ;
c¼
1 : jmj
Now, we perform a rotation R of the coordinate system in such a way that a :¼ R~ a ¼ j~ ajð1; 0; 0ÞT : At the end we denote u0 :¼ R~ u0 and a ¼ j~ aj. T
Lemma 6. Assume a; c 2 R and u0 ¼ ðx0 ; y0 ; z0 Þ 2 R3 . Then Eq. (38) for T T a ¼ að1; 0; 0Þ admits the solution uðtÞ ¼ ðxðtÞ; yðtÞ; zðtÞÞ , which is given by eactju0 j ðju0 j þ x0 Þ eactju0 j ðju0 j x0 Þ ; eactju0 j ðju0 j þ x0 Þ þ eactju0 j ðju0 j x0 Þ y0 cosðatÞ z0 sinðatÞ yðtÞ ¼ 2ju0 j actju j ; 0 e ðju0 j þ x0 Þ þ eactju0 j ðju0 j x0 Þ y0 sinðatÞ þ z0 cosðatÞ zðtÞ ¼ 2ju0 j actju j : e 0 ðju0 j þ x0 Þ þ eactju0 j ðju0 j x0 Þ xðtÞ ¼ ju0 j
Proof. See [22].
h
Further we follow the same way as in Section 3, i.e., we show the existence of ei ; hi for all i ¼ 1; . . . ; n and we prove the same a priori estimates as in Lemmas 4 and 5. Then we can derive the same convergence rates as we did in Theorem 1. We state the result without the proof. Theorem 2. Let the assumptions of Lemma 5 be satisfied. There exist positive constants C and s0 such that (i) maxt2½0;T kEðtÞ en ðtÞk2 þ maxt2½0;T kHðtÞ hn ðtÞk2 6 Cs, (ii) maxt2½0;T kMðtÞ mn ðtÞk2 þ maxt2½0;T kot M ot mn k2 6 Cs, hold for any 0 < s < s0 .
5. Numerical experiment In this section we show one example which demonstrates some of the features which make ferromagnetic layers interesting. We consider a rectangular parallelepipedon with the length Lx ¼ 16, the width Ly ¼ 23 and the height Lz ¼ 0:4 (see Fig. 1). This domain is occupied by a
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
95
Lz
z y x
Ly Lx Fig. 1. Computational domain X.
ferromagnet. We split the domain into small blocks with Dx ¼ Dy ¼ Dz ¼ 152 . Afterwards we divide each one into 6 tetrahedrons. For the approximation of the fields E and H we use WhitneyÕs edge elements, cf. [2,4]. Recall that the approximation of M can be settled down at any point x from the approximation of the LLG equation––see (36), which is used for computations. We have determined m at the barycenter of each tetrahedron, i.e., the approximation of M is a piecewise constant vector field. We consider the following boundary conditions
T left: m Eðt; x; y; zÞ ¼ m ð0; 0; 1Þ E0 g f t xc , right: m E ¼ 0, bottom: m E ¼ 0, top: m E ¼ 0, front: m H ¼ 0, back: m H ¼ 0, where E0 ¼ 4:0104 and g is the Ricker pulse (see Fig. 2) 8 < 2k ð2kðs 1:35Þ2 1Þ expðkðs 1:35Þ2 Þ if t < 3:48; gðsÞ ¼ 11 : 0 if t P 3:48 p 2 and k ¼ 1:31 . The values for f and c are given in Table 1. Let us note that we are using the following notation for the coordinates: xdashed, y-gray and z-solid. The initial data are given as follows: E 0 ¼ H 0 ¼ 0 T and M 0 ¼ ð0; 0; 100Þ , which ensures that r B ¼ 0. The computed results for the fields E; H and M along the horizontal center line of the domain are shown in Figs. 3–5, respectively. The waves correspond to two different time points t ¼ 6:705 108 s and t ¼ 1:3205 107 s. The incoming wave for E at the left boundary part has vanishing x- and y-components. When the wave enters the ferromagnetic layer, it is rotated, which makes x- and y-components nonzero. Moreover, due to the fact that a 6¼ 0 is the layer absorbing.
96
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
10000 0 –10000 –20000 –30000 –40000 –20
–15
–5
–10
0
s Fig. 2. Ricker pulse g.
Table 1 Parameters Parameter
Value
Parameter
Value
K e0 c s
0 0.88422 · 1011 Fm1 3.0 · 108 ms1 0.5 · 1010 s
a l0 c f
0.2 1.25667 · 106 Hm1 2.2 · 105 Hz mA1 4.0 · 107 Hz
40000
40000
20000
20000
0
0
–20000
–20000
–40000
–40000 2
4
6
8
10 12 14 16
2
4
6
8
10 12 14 16
Fig. 3. Left: Incoming E-wave at t ¼ 6:705 108 s. Right: Reflected E-wave at t ¼ 1:3205 107 s.
The electric field E in Fig. 3 has a dominant z-coordinate. The incident wave is traveling to the right and it is reflected and rotated by the right wall. Its amplitude is decreasing in time. The magnetic field H in Fig. 4 has a dominant
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99 100
100
50
50
0
0
–50
–50
–100
–100 2
4
6
8
Fig. 4. Left: Incoming t ¼ 1:3205 107 s.
2
10 12 14 16
H-wave
at
8
t ¼ 6:705 10
100
100
80
80
60
60
40
40
20
20
0
0 2
4
6
8
10 12 14 16
2
Fig. 5. Left: Incoming M-wave at t ¼ 6:705 108 t ¼ 1:3205 107 s.
4
s.
6
8
Right:
4
6
10 12 14 16
Reflected
8
97
H-wave
at
10 12 14 16
s. Right: Reflected M-wave at
y-coordinate, but here the reflected wave is not rotated. The amplitude is also decreasing. The magnetization M shows only small changes in Fig. 5. The initial datum was Mð0Þ ¼ ð0; 0; 100ÞT . Due to the fact that the magnitude of M remains constant, we can see that the increase of jxj; jyj-components implies the drop of the z-coordinate.
E
H H v
v E
Fig. 6. Position among E, H and the velocity vector v. Left: Incoming wave. Right: Reflected wave.
98
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
The mutual position among E; H and the velocity vector v of the wave is fixed. Therefore, when the reflected wave is coming back, then the wave for H keeps its orientation and E is ‘‘upside down’’ (see Fig. 6).
Acknowledgements The authors would like to express their thanks to Prof. R. Hiptmair (ETH, Zurich, Switzerland) and Dr. L. Dupre (Ghent University, Belgium) for many encouraging and enlightening discussions. References [1] G. Anzellotti, S. Baldo, A. Visintin, Asymptotic behavior of the Landau–Lifshitz model of ferromagnetism, Appl. Math. Optim. 23 (3) (1991) 171–192. [2] A. Bossavit, Computational Electromagnetism, Variational Formulations, Complementarity, Edge Elements, in: Electromagnetism, vol. XVIII, Academic Press, Orlando, FL, 1998. [3] M. Brokate, J. Sprekels, Hysteresis and Phase Transitions, in: Applied Mathematical Sciences, vol. 121, Springer, 1996. [4] M. Cessenat, Mathematical Methods in Electromagnetism, Linear Theory and Applications, in: Advances in Mathematics for Applied Sciences, vol. 41, World Scientific Publishers, Singapore, 1996. [5] Y. Chen, Dirichlet boundary value problems of Landau–Lifshitz equation, Commun. Partial Differ. Eqs. 25 (1–2) (2000) 101–124. [6] Y. Chen, Existence and singularities for the Dirichlet boundary value problems of Landau– Lifshitz equations, Nonlinear Anal., Theory Meth. Appl. 48A (3) (2002) 411–426. [7] W. E, X.-P. Wang, Numerical methods for the Landau–Lifshitz equation, SIAM J. Numer. Anal. 38 (5) (2000) 1647–1665. [8] B. Guo, S. Ding, Neumann problem for the Landau–Lifshitz–Maxwell system in two dimensions, Chin. Ann. Math. Ser. B 22 (4) (2001) 529–540. [9] B. Guo, F. Su, Global weak solution for the Landau–Lifshitz–Maxwell equation in three space dimensions, J. Math. Anal. Appl. 211 (1) (1997) 326–346. [10] H. Haddar, P. Joly, Effective boundary conditions for thin ferromagnetic layers: The onedimensional model, SIAM J. Appl. Math. 61 (4) (2001) 1386–1417. [11] J.L. Joly, G. Metivier, J. Rauch, Electromagnetic wave propagation in presence of a ferromagnetic material, ESAIM, Proc. 3 (1998) 85–99. [12] J.L. Joly, G. Metivier, J. Rauch, Global solutions to Maxwell equations in a ferromagnetic medium, Ann. Henri Poincare 1 (2) (2000) 307–340. [13] P. Joly, A. Komech, O. Vacus, On transitions to stationary states in a Maxwell–Landau– Lifschitz–Gilbert system, SIAM J. Math. Anal. 31 (2) (2000) 346–374. [14] P. Joly, O. Vacus, MaxwellÕs equations in a 1D ferromagnetic medium: existence and uniqueness of strong solutions, INRIA, technical report no. 3052, 1996. [15] P. Joly, O. Vacus, Mathematical and numerical studies of nonlinear ferromagnetic materials, M2AN 33 (3) (1999) 593–626. [16] L.D. Landau, E.M. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies, Phys. Z. Sowjetunion 8 (1935) 153–169. [17] L.D. Landau, E.M. Lifshitz, Electrodynamics of Continuous Media, Pergamon Press, Oxford, X, 1960, p. 417 (J.B. Sykes, J.S. Bell, Translated from the Russian).
M. Slodicka, L. Banas / Appl. Math. Comput. 158 (2004) 79–99
99
[18] I. Mayergoyz, Nonlinear Diffusion of Electromagnetic Fields with Applications to Eddy Currents and Superconductivity, Academic Press, London, 1998. [19] P.B. Monk, O. Vacus, Error estimates for a numerical scheme for ferromagnetic problems, SIAM J. Numer. Anal. 36 (3) (1998) 696–718. [20] P.B. Monk, O. Vacus, Accurate discretization of a nonlinear micromagnetic problem, Comput. Meth. Appl. Mech. Eng. 190 (40–41) (2001) 5243–5269. [21] A. Prohl, Computational Micromagnetism, in: Advances in Numerical Mathematics, vol. xvi, Teubner, Leipzig, 2001. [22] M. Slodicka, I. Cimrak, Numerical study of nonlinear ferromagnetic materials, Appl. Numer. Math. 46 (1) (2003) 95–111. [23] F. Su, B. Guo, The global smooth solution for Landau–Lifshitz–Maxwell equation without dissipation, J. Partial Differ. Eqs. 11 (3) (1998) 193–208. [24] A. Visintin, On Landau–LifshitzÕ equations for ferromagnetism, Jpn. J. Appl. Math. 2 (1985) 69–84. [25] A. Visintin, Differential Models of Hysteresis, in: Applied Mathematical Sciences, vol. 111, Springer, 1994. [26] X.-P. Wang, C.J. Garcıa-Cervera, W. E, A Gauss–Seidel projection method for micromagnetics simulations, J. Comput. Phys. 171 (1) (2001) 357–372. [27] P. Weiss, LÕhypothese du champ moleculaire et la propriete ferromagnetique, J. Physique 6 (1907) 661–690. [28] J. Zhai, Existence and behavior of solutions to the Landau–Lifshitz equation, SIAM J. Math. Anal. 30 (4) (1999) 833–847.