4.5 Solution of the time-stepping equation

The properties we have just discussed are strongly related to the implicit nature of mid-point rule. As consequence of this implicit nature, we have to solve the time-stepping Eq. (4.33) for the unknown $\underline{\textbf{m}}^{n+1}$ at each time step which amounts to solve the following system of $3N$ nonlinear equations in the $3N$ unknowns $\underline{\textbf{m}}^{n+1}$:

$$\mathbf{F}(\underline{\textbf{m}}^{n+1})=\mathbf{0}\quad,$$ (4.43)

where $\mathbf{F}(\underline{{\textbf{y}}}):\mathbb{R}^{3N}\rightarrow \mathbb{R}^{3N}$ is the following vector function:

$$\mathbf{F}(\underline{{\textbf{y}}})=\left[{\underline{\textrm{I}... ...}\left(\frac{\underline{{\textbf{y}}}+\underline{\textbf{m}}^{n}}{2}\right) ,$$ (4.44)

and where

$$\mathbf{f}(\underline{\textbf{m}})=-{\underline{\Lambda}}(\underline{\textbf{m}}) \cdot \underline{\textbf{h}}_ (\underline{\textbf{m}})= {\underline{\Lambda}}(\underline{\textb... ...t \frac{\partial \underline{{\text{g}}}}{\partial \underline{\textbf{m}}} ,$$ (4.45)

is the right-hand-side of the conservative LLG equation. It is interesting to notice that the damping is present in only one term in the function $\mathbf{F}(\underline{{\textbf{y}}})$ and it introduce only a slight modification of the function. The solution of the system of equation (4.43) can be obtained by using Newton-Raphson (NR) algorithm. To this end, we derive the jacobian matrix ${\underline{\textrm{J}}}_$F$(\underline{{\textbf{y}}})$ of the vector function $\mathbf{F}(\underline{{\textbf{y}}})$ which, after simple algebraic manipulations, can be written in the following form

$${\underline{\textrm{J}}}_ (\underline{{\textbf{y}}})={\underline{\textrm{I}}}-\alpha{\under... ...(\underline{\textbf{m}}^{n})- \frac{\Delta t}{2} {\underline{\textrm{J}}}_ \left(\frac{\underline{{\textbf{y}}}+\underline{\textbf{m}}^{n}}{2}\right)$$ (4.46)

where ${\underline{\textrm{J}}}_$f is the jacobian matrix associated to $\mathbf{f}(\underline{\textbf{m}})$. By using the Eqs.(4.25), (4.45), one obtains:

$${\underline{\textrm{J}}}_ (\underline{\textbf{m}})=\frac{\partial \mathbf{f}}{\partial \und... ...ne{\textrm{C}}}\cdot\underline{\textbf{m}}+\underline{\textbf{h}}_a\right) .$$ (4.47)

The main difficulty in applying NR method is that the Jacobian ${\underline{\textrm{J}}}_$F$(\underline{{\textbf{y}}})$ of $\mathbf{F}(\underline{{\textbf{y}}})$ is a full matrix, due to the long-range character of magnetostatic interactions which reflects in the full nature of the matrix ${\underline{\textrm{C}}}$. In this connection, let us observe that the damping term affect only a small sparse component of the jacobian ${\underline{\textrm{J}}}_$F$(\underline{{\textbf{y}}})$ and thus does not introduce any basic difficulty. Anyhow, due to the full nature of ${\underline{\textrm{J}}}_$F$(\underline{{\textbf{y}}})$, the use of the plain NR method would require an unpractical computational cost. However, as it is usual in solving field problems with implicit time stepping, we can use a quasi-Newton method by considering a reasonable approximation of the Jacobian. In order to have a sparse Jacobian one can consider the following expression ${\underline{\tilde{\textrm{J}}}}_$F in which magnetostatic interactions are not included:

$${\underline{\tilde{\textrm{J}}}}_ (\underline{{\textbf{y}}})={\underline{\textrm{I}}}-\alpha{\under... ...ine{\textbf{m}}^{n})- \frac{\Delta t}{2} {\underline{\tilde{\textrm{J}}}}_ \left(\frac{\underline{{\textbf{y}}}+\underline{\textbf{m}}^{n}}{2}\right) \quad,$$ (4.48)

where the matrix ${\underline{\tilde{\textrm{J}}}}_$f is

$${\underline{\tilde{\textrm{J}}}}_ (\underline{\textbf{m}})=-{\underline{\Lambda}}(\underline{\textbf{m}})\cdot (-({\underline{\textrm{C}}}_ +{\underline{\textrm{C}}}_ ))+{\underline{\Gamma}}\cdot \left[-({\underline{\textrm{C}}}_\t... ...}}_\text{an})\cdot\underline{\textbf{m}}+\underline{\textbf{h}}_a\right] \quad.$$ (4.49)

Basically, the latter equation is obtained by substituting the full matrix ${\underline{\textrm{C}}}$ with its sparse component ${\underline{\textrm{C}}}_$ex$+{\underline{\textrm{C}}}_$an in Eq. (4.47). Thus, the iterative procedure can be summarized as follows:

$$\underline{{\textbf{y}}}_0=\underline{\textbf{m}}^n , ... ...xtbf{y}}}_{k+1}=\underline{{\textbf{y}}}_k+\Delta\underline{{\textbf{y}}}_{k+1} \quad {\underline{\tilde{\textrm{J}}}}_ (\underline{{\textbf{y}}}_k)\Delta\underline{{\textbf{y}}}_{k+1}=-\mathbf{F}(\underline{{\textbf{y}}}_k) .$$ (4.50)

At each iteration, the linear system defined by the matrix (4.49) has to be inverted. Since this matrix is non-symmetric, we have found appropriate to use generalized minimum residual (GMRES) method [91]. The iteration (4.50) is repeated until the norm $\Vert\mathbf{F}(\underline{{\textbf{y}}}_k)\Vert$ is under a prescribed tolerance. The iterative technique we developed to solve Eq. (4.43) belongs to the main category of quasi-Newton methods. In this respect, it has been proven [92] that this kind of iterative procedure is convergent and the order of convergency is the first order, provided that the initial guess is sufficiently close to the `true' solution of the system.