Fourier-based numerical approximation of the Weertman equation for moving dislocations

TL;DR

This paper presents a Fourier-based numerical method for approximating solutions to the Weertman equation, modeling moving dislocations, using a time-dependent formulation and a preconditioned collocation scheme for efficiency.

Contribution

It introduces a novel Fourier transform-based collocation scheme with large time steps for solving the nonlinear Weertman equation.

Findings

Efficient numerical approximation of the Weertman equation.

Robust method capturing solutions with large time steps.

Applicable to various nonlinearities in dislocation modeling.

Abstract

This work discusses the numerical approximation of a nonlinear reaction-advection-diffusion equation, which is a dimensionless form of the Weertman equation. This equation models steadily-moving dislocations in materials science. It reduces to the celebrated Peierls-Nabarro equation when its advection term is set to zero. The approach rests on considering a time-dependent formulation, which admits the equation under study as its long-time limit. Introducing a Preconditioned Collocation Scheme based on Fourier transforms, the iterative numerical method presented solves the time-dependent problem, delivering at convergence the desired numerical solution to the Weertman equation. Although it rests on an explicit time-evolution scheme, the method allows for large time steps, and captures the solution in a robust manner. Numerical results illustrate the efficiency of the approach for several types of nonlinearities.

Tables

Table 1. Table 1: Errors as a function of Δ ​ t Δ 𝑡 \Delta t ( 2 ​ L = 638 2 𝐿 638 2L=638 , 2 ​ m = 4096 2 𝑚 4096 2m=4096 , h = 0.156 ℎ 0.156 h=0.156 , σ 1 = 0.9921 subscript 𝜎 1 0.9921 \sigma_{1}=0.9921 ).

$\Delta t$	0.25	0.1	0.01	0.001
$\Vert \Delta \eta \Vert$	$5.39992\times {10}^{-5}$	$5.39991\times {10}^{-5}$	$5.39990\times {10}^{-5}$	$5.39990\times {10}^{-5}$
$\Vert \Delta c/c\Vert$	$2.01277\times {10}^{-5}$	$2.01266\times {10}^{-5}$	$2.01260\times {10}^{-5}$	$2.01260\times {10}^{-5}$

Table 2. Table 2: Error on η 𝜂 \eta minimized over L 𝐿 L for σ = σ 2 𝜎 subscript 𝜎 2 \sigma=\sigma_{2} , as a function of m 𝑚 m .

$2m$	$32$	$128$	$512$	$2048$	$8192$
min. value of $\Vert \Delta \eta \Vert$	$2.1\times {10}^{-2}$	$4.9\times {10}^{-4}$	$5.5\times {10}^{-5}$	$6.8\times {10}^{-6}$	$8.5\times {10}^{-7}$
optimum $hL/a{\left({\sigma }_2\right)}^2$	$0.781$	$0.195$	$0.195$	$0.195$	$0.195$

Equations

$$\{\begin{array}{rl}& -\left|{\partial }_x\right|\eta (x)+c_{\eta }{\partial }_x\eta (x)=F_{\sigma }^{\prime }(\eta (x))\text{for}x\in \mathbb{R},\\ & \eta (-\infty )={\eta }_{\mathrm{l}}\text{and}\eta (+\infty )={\eta }_{\mathrm{r}},\end{array}$$

$$|{\partial }_x|\eta (x)=\mathcal{H}({\partial }_x\eta )(x)=\frac{1}{\pi }p.v.{\int }_{-\infty }^{+\infty }\frac{{\partial }_x\eta (x^{\prime })}{x-x^{\prime }}\mathrm{d}{x}^{\prime }=\underset{ϵ\to 0}{\lim }\frac{1}{\pi }{\int }_{|x-x^{\prime }|>ϵ}\frac{{\partial }_x\eta (x^{\prime })}{x-x^{\prime }}\mathrm{d}{x}^{\prime },$$

$$|{\partial }_x|\eta (x)=-\frac{1}{\pi }{\int }_0^{+\infty }\frac{\eta (x+y)+\eta (x-y)-2\eta (x)}{y^2}\mathrm{d}y.$$

$$\mathcal{F}\{f\}(k)=\hat{f}(k)={\int }_{-\infty }^{+\infty }{\mathrm{e}}^{-\mathrm{i}kx}f(x)\mathrm{d}x.$$

$$\{\begin{array}{rl}& -\left|{\partial }_x\right|\eta (x)=F_{\sigma }^{\prime }(\eta (x))\text{for}x\in \mathbb{R},\\ & \eta (-\infty )={\eta }_{\mathrm{l}}\text{and}\eta (+\infty )={\eta }_{\mathrm{r}}.\end{array}$$

$$\{\begin{array}{rl}& {\partial }_{t}u(t,x)-c(t){\partial }_{x}u(t,x)+|{\partial }_x|u(t,x)=-F_{\sigma }^{\prime }(u(t,x)),\\ & u(t=0,x)=u_0(x),\end{array}$$

$$u_0(-\infty )$$

$$F_{\sigma }(\eta ):=\frac{1}{4\pi }\left[1-{\theta }^2-\left(\theta \sqrt{1-{\theta }^2}+\arcsin (\theta )-\varphi \right)\cot \varphi \right]-\sigma \eta ,$$

$$\theta =\sin (\varphi )\cos (2\pi \eta ),\varphi =\arctan r.$$

$$\eta (x)$$

$$|{\partial }_x|\eta (x)$$

$$F_{\sigma }^{\prime }(\eta )=F_{\sigma }^{\prime }({\eta }_{\mathrm{l},\mathrm{r}})+F_{\sigma }^{\prime \prime }({\eta }_{\mathrm{l},\mathrm{r}})(\eta -{\eta }_{\mathrm{l},\mathrm{r}})+\dots \simeq F_{\sigma }^{\prime \prime }({\eta }_{\mathrm{l},\mathrm{r}})(\eta -{\eta }_{\mathrm{l},\mathrm{r}}),$$

$$a_{\mathrm{l},\mathrm{r}}$$

$$\eta (x)$$

$$\overline{a}:=\frac{1}{2}(a_{\mathrm{l}}+a_{\mathrm{r}})=\pi \left[\frac{1}{F_{\sigma }^{\prime \prime }({\eta }_l)}+\frac{1}{F_{\sigma }^{\prime \prime }({\eta }_r)}\right].$$

$$c_{\eta }$$

$${\partial }_t\phi (t,x)+|{\partial }_x|\phi (t,x)=-F_{\sigma }^{\prime }(\phi (t,x))\text{with}\phi (t=0,x)=u_0(x),$$

$$u(t,x)$$

$$\underset{x\in \mathbb{R}}{\sup }|\phi (t,x)-\eta (x-c_{\eta }t+\xi )|\le K{\mathrm{e}}^{-\kappa t}.$$

$$u(t,x)$$

$$\frac{1}{2L}{\int }_{-L}^L\eta (x)\mathrm{d}x=\overline{\eta },$$

$$\overline{\eta }:=({\eta }_{\mathrm{l}}+{\eta }_{\mathrm{r}})/2,$$

$$c(t):=\frac{F_{\sigma }({\eta }_{\mathrm{r}})-F_{\sigma }({\eta }_{\mathrm{l}})}{{\int }_{-\infty }^{+\infty }[{\partial }_{x}u(t,x){]}^2\mathrm{d}x}+\frac{\kappa }{({\eta }_{\mathrm{l}}-{\eta }_{\mathrm{r}})}I(t),\text{where}I(t):={\int }_{-L}^L\left[u(t,x)-\overline{\eta }\right]\mathrm{d}x,$$

$$\frac{F_{\sigma }({\eta }_{\mathrm{r}})-F_{\sigma }({\eta }_{\mathrm{l}})}{{\int }_{-\infty }^{+\infty }[{\partial }_{x}u(t,x){]}^2\mathrm{d}x}\underset{t\to +\infty }{⟶}{c}_{\eta }.$$

$$c(t)\simeq c_{\eta }+\frac{\kappa }{({\eta }_{\mathrm{l}}-{\eta }_{\mathrm{r}})}I(t).$$

$$\frac{1}{2L}{\int }_{-L}^{L}u(t,x)\mathrm{d}x$$

$${\partial }_{t}u(t,x)$$

$$\dot{I}(t)$$

$$F_{\sigma }^{\prime }(\eta )=\sin \left(2\pi \eta \right)-\sigma ,$$

Full text

Weertman equation, Peierls-Nabarro equation, dislocations, Cauchy-type nonlinear integrodifferential equation, reaction-advection-diffusion equation, fractional Laplacian, preconditioned scheme.

I Introduction

This article addresses the numerical approximation of the following nonlinear integrodifferential equation, with Cauchy-type singular kernel:

[TABLE]

where both the real-valued function $\eta$ and the scalar constant $c_{\eta}$ are the unknowns. The potential $F_{\sigma}$ is a nonlinear bistable function of $\eta$ with (at least) two local minima at values $\eta=\eta_{\rm l}$ and $\eta=\eta_{\rm r}$. The meaning of the subscript $\sigma$ is explained below. The operator $|\partial_{x}|$ is linear, and defined in terms of the Hilbert transform $\mathcal{H}$ KING09 as

[TABLE]

where ${\rm p.v.}$ denotes the principal value KANW04 at $x$. In the context of singular integral equations, the above kernel of the Hilbert transform is known as the Cauchy kernel. The operator $-|\partial_{x}|$, also denoted $-(-\Delta)^{1/2}$ by some authors, is diffusive (Brezis, , p. 181). Another useful representation of (2) is (by integration by parts)

[TABLE]

Let the Fourier transform in the continuum (FT) of a function $f$ be defined at wavemode $k$ as

[TABLE]

One has $\mathcal{F}\left\{\text{p.v.}\,x^{-1}\right\}(k)=-{\rm i}\,\pi\,{\rm{sgn}}(k)$ (Gradshteyn, , p. 1118), whence $\mathcal{F}\left\{\,|\partial_{x}|\,\eta\right\}(k)$ $=$ $|k|\widehat{\eta}(k)$. Thus, the non-local operator $|\partial_{x}|$ is symmetric and positive.

Equation (1) is a dimensionless form of the Weertman equation WEER69a ; WEER69b ; Rosakis (simply referred to as ‘the Weertman equation’ in the following), which models straight dislocations traveling with steady velocity, thus generalizing the Peierls-Nabarro (PN) equation for static dislocations Nabarro46 :

[TABLE]

Dislocations are linear defects in crystals, the motion of which is responsible for the plasticity of metals Hirth . Dislocation lines have a non-vanishing sectional area, i.e., they possess a ‘core’ of finite width. The derivative $\partial_{x}\eta(x)$ (the so-called dislocation density) of the unknown function $\eta$ in (1) describes the shape function of a flat finite-width dislocation on its glide plane, along the $x$-direction. The core is the region of space where $\partial_{x}\eta(x)$ develops peaks. From a physical standpoint, the function $\eta$ represents a local relative material displacement discontinuity between the upper and lower half-spaces surrounding the glide plane on which moves the dislocation line; see, e.g., Hirth for details. From a broader perspective, the function $\eta$ can be understood as a moving phase-transformation front between the states $\eta_{\rm l}$ and $\eta_{\rm r}$ (Fig. 1).

In (1) the term $|\partial_{x}|\eta$ accounts for the long-range elastic self-interactions that tend to spread the core. This repulsive interaction is counterbalanced by the nonlinear pull-back force $F_{\sigma}^{\prime}(\eta)$, which binds together the upper and lower half-spaces, thus giving the dislocation core its finite width. Throughout this article, we consider that $F_{\sigma}^{\prime}(\eta)$ includes a constant externally applied loading $\sigma$ (that is, $F_{\sigma}(\eta)=F(\eta)-\sigma\eta$ where $F(\eta)$ is an energy potential intrinsic to the material, tilted by the adjunction of a linear term $-\sigma\eta$). Moreover, the moving dislocation is subjected to various drag mechanisms encoded into the term $c_{\eta}\,\partial_{x}\eta$. As recalled below, the Weertman equation admits an analytical solution Rosakis if $F(\eta)$ is assumed sinusoidal. For more realistic potentials, a numerical approach is required.

Since its inception, the original PN model has been enriched in various directions. For instance, it most often requires being generalized to vector-valued $\eta$ to be quantitatively predictive SydowHartfordWahnstrom ; LuBulatovKioussis ; LuBulatovKioussis2 ; DENO04 ; DENO07 . Also, the model has been extended to two dimensions of space, to study planar dislocation loops DENO04 ; DENO07 ; Xiang ; zhu2015 . Quite generally, methods to compute the shape of static or moving cores encompass variational approaches and involve finite-element and/or phase-field-type implementations DENO04 ; DENO07 ; ZHAN15 ; Mianroodi . Yet, in spite of such a wealth of enrichments of the PN model and its associated numerical methods of solution, the one-dimensional Weertman equation —a comparatively simpler extension— has not been investigated as thoroughly, while the specific problem of determining the allowed velocities of steadily-moving dislocations for general force laws $F_{\sigma}^{\prime}$ remains an open question of major practical interest Rosakis . For this reason, the present work focuses on solving the simplest, scalar, and one-dimensional case. A generalization to the vector case of the Weertman equation is the subject of ongoing work, and will be presented elsewhere YPP2 .

It must be emphasized that, in the dimensionless form (1) of the equation, $c_{\eta}$ is not the physical velocity of the dislocation. The latter is deduced from $c_{\eta}$ in a post-processing step to be explained in YPP2 (it is a direct application of Equations (46), (47) or (48) of YPP ), which in the present scalar case is independent from the numerical task of solving the equation. Therefore the physical velocity is not further considered hereafter. Moreover, we stress that (1) applies only to subsonic motion, for which the coefficient of $|\partial_{x}|$ does not vanish in the original Weertman equation.

Numerical methods for solving integrodifferential equations such as (1) have been proposed by many authors. The method employed in Lejvcek uses properties of the Hilbert transform to recast Equation (5) into a form amenable to fixed-point approaches. In Kurtz the authors consider a simpler version of (1) on a bounded domain, in which the nonlinear term $F_{\sigma}^{\prime}(\eta(x))$ is replaced by some given $\eta$-independent function $g(x)$, and where $c_{\eta}$ is also given. The solution is then obtained by means of a collocation method with quadratic interpolation. Those works make use of the expression of the operator $|\partial_{x}|$ in the direct space. More recently, Karlin et al. presented Movchan a general iterative method for solving (5), based on the expression of $|\partial_{x}|\eta$ in the Fourier space. Our work borrows from the latter reference. The interested reader can also refer to BUEN14 for Fourier-based numerical schemes applied to the fractional Laplacian operator $|\partial_{x}|^{\alpha}$ with $\alpha>0$.

The present article proposes a numerical method to approximate solutions to (1) in the case where $F_{\sigma}$ is bistable. As in Movchan , we build a dynamical system that admits (1) as its long-time limit, namely,

[TABLE]

where $x\in\mathbb{R}$, and $u_{0}$ is a regular initial data taking values in the interval $\left[\eta_{\rm r},\eta_{\rm l}\right]$ such that

[TABLE]

where $\eta_{\rm l,r}$ are zeros of $F^{\prime}_{\sigma}(\eta)$.

The iterative numerical approach introduced below uses (6) to approximate the solution of (1). It is immediate that if $(\eta,c_{\eta})$ solves (1) and if we impose $c(t):=c_{\eta}$, then $u(t,x):=\eta(x)$ satisfies (6) for the initial data $u_{0}=\eta$. It is proved in MathWeert that if $F_{\sigma}$ is bistable, then for any initial data $u_{0}$ with values in $\left[\eta_{\rm r},\eta_{\rm l}\right]$ that satisfies (7) and for any continuous function $c(t)$, the solution $u$ of (6) converges to the solution of (1) in a sense explained in Sec. II.5 below. This leads to the following procedure:

consider an initial condition $u_{0}(x)$ such that (7) holds; 2. 2.

approximate the solution $u(t,x)$ to (6); 3. 3.

while evaluating $u(t,x)$, choose $c(t)$ such that the core of $u(t,x)$ remains within the computational box and that $c(t)$ converges to $c_{\eta}$; 4. 4.

for $t=t_{\text{f}}$ sufficiently large so that $u\left(t,x\right)$ has numerically converged, return $u\left(t_{\text{f}},x\right)$ as a numerical approximation to $\eta$, and deduce the velocity $c_{\eta}$.

As will be shown, this strategy proves efficient in cases of physical interest. However there might exist alternative strategies for solving (1).

Numerical approximation of (1) paves the way to investigating numerically the Dynamic PN equation YPP , which generalizes the Weertman equation to transient regimes. Indeed, the initial conditions and long-time steady-state regimes of the Dynamic PN equation are solutions to the latter equation YPP1 . However, we emphasize that (6) is only an algorithmic tool that has no relationship whatsoever with the actual dynamics of dislocations.

The article is organized as follows. In Sec. II, we briefly study the Weertman equation, and discuss both the uniqueness of its solution and its interpretation as the long-time limit of the dynamical system (6). We formally derive asymptotes of solutions to (1) and state identities about the velocity $c_{\eta}$ in general cases. Also, we explain how to choose $c(t)$ in (6) to solve this equation in a comoving frame, namely, one which follows the translational motion of the core. An analytical solution to (1) that exists in a particular case is recalled. In Sec. III, we introduce our numerical representation for $\eta$ and discuss corresponding implementations of the diffusive operator $-|\partial_{x}|$ and the advective operator $\partial_{x}$, as well as methods to evaluate $c(t)$. Also, we make use of the asymptotic behavior when $|x|\rightarrow+\infty$ of the solution to (1) to circumvent the issue of the infinite domain of integration in (2). Once these fundamental elements have been introduced, we build in Sec. IV a Preconditioned Collocation Scheme (PCS) that applies to our problem, and justify this denomination. In Sec. V, we use this numerical approach on two test cases: one with a simple potential $F_{\sigma}$, for which the exact solution is known, and one with a more physically relevant potential $F_{\sigma}$. We also illustrate the robustness of our approach, concluding that the algorithm presented is unconditionally stable with respect to the time step $\Delta t$. We empirically derive error scalings with respect to the parameters involved in the discretization. A concluding discussion closes the article, underlining some limitations of our approach, and proposing a few possible extensions. An Appendix is devoted to examining further one such extension.

II Some properties of the Weertman equation

This section is devoted to an overview of some important properties of the Weertman equation (1) and of the companion dynamical equation (6).

II.1 Invariances

As the PN equation, Equation (1) is obviously invariant by translation. When invoking uniqueness of the solutions, we shall henceforth implicitly refer to ‘uniqueness up to arbitrary translations’. This invariance has consequences on the numerical solution, which usually undergoes an undesirable drift during calculations if no corrective action is undertaken. A special procedure is developed in Section II.6 below to eliminate this difficulty.

Moreover, Eq. (1) is invariant by reflection in the sense that if $(\eta(x),c_{\eta})$ is a solution for boundary conditions $\eta(\pm\infty)=\eta_{\rm l,r}$, then $(\eta(-x),-c_{\eta})$ is another solution for boundary conditions $\eta(\pm\infty)=\eta_{\rm r,l}$. Therefore, without loss of generality, we always assume throughout this article that $\eta_{\rm l}>\eta_{\rm r}$.

II.2 Existence and uniqueness of solutions to the Weertman equation.

There exists a unique solution to (1) when $F_{\sigma}$ is a bistable nonlinearity. More precisely it can be shown rigorously (the proof relies on a recent result Gui ) that for $F_{\sigma}$ sufficiently regular, if $\eta_{\rm l}>\eta_{\rm r}$ are such that: (i) $F_{\sigma}^{\prime}(\eta_{\rm l,r})=0$ and $F_{\sigma}^{\prime\prime}(\eta_{\rm l,r})>0$; (ii) any local minimizer $u$ of $F_{\sigma}$ between $\eta_{\rm r}$ and $\eta_{\rm l}$ satisfies $F_{\sigma}(u)>F_{\sigma}(\eta_{\rm r})$ and $F_{\sigma}(u)>F_{\sigma}(\eta_{\rm l})$, then there exists a unique velocity $c_{\eta}$ and a decreasing function $\eta$ satisfying (1), which is unique up to translation. Condition (i) means that $F_{\sigma}$ is a ‘bistable potential’. A typical example of such $F_{\sigma}$, to be used in Sec. V.5, is drawn in Fig 2. It corresponds to

[TABLE]

Possible non-decreasing solutions to (1) that might exist as in the PN equation (Nabarro46, , Equation (16)), e.g., for $\eta_{\rm l}=\eta_{\rm r}$, are disregarded in the present work. Hereafter, we always assume that $\eta_{\rm l,r}$ obey the above conditions, and we term such values of $\eta(x)$ at infinity consistent boundary conditions (CBCs).

II.3 Asymptotic behavior and characteristic lengths

By letting $|x|\to\infty$, we formally deduce the leading-order asymptotic expansions

[TABLE]

They are proved rigorously in MathWeert . The key ingredient of the proof is the following asymptotic behavior of integrals with Cauchy kernel (see (Musk, , p. 267)):

[TABLE]

It can be formally retrieved from the leading term of a ‘series expansion’ of (2) at $x$ large (note however that the integral involved in the next-to-leading term of such a formal expansion may not exist). Substituting expression (10) into (1), using the Taylor expansions at boundary values

[TABLE]

and noting that $\partial_{x}\eta(x)$ vanishes like $O\left(|x|^{-2}\right)$ Gui , we obtain (9).

Introducing the characteristic lengths

[TABLE]

Equation (9) can be rewritten as

[TABLE]

It will be argued in Sec. II.7 that $a_{\rm l,r}$ represent typical scales of variation of the dislocation density on both sides of the solution. For further purposes, it is useful to introduce a mean characteristic length as

[TABLE]

Finally, depending on the potential $F_{\sigma}(\eta)$ at hand, the solution can be either symmetric, in the sense that $\partial_{x}\eta(x)$ is an even function (an example will be given in Sec. II.7), or non-symmetric if the latter property does not hold.

II.4 Velocity determination

There are many ways to determine the velocity $c_{\eta}$ associated with the solution $\eta$ to Equation (1). One possibility is to multiply (1) by a function $g$ chosen such that the integrals involved make sense, and to integrate over $\mathbb{R}$. This yields

[TABLE]

In particular, the following choices of $g$ eliminate $|\partial_{x}|\eta(x)$. With $g(x)=1$ one gets

[TABLE]

where the improper integral is understood as a principal value at infinity (Gradshteyn, , p. 252). Taking instead $g(x)=\partial_{x}\eta(x)$ Gui , one arrives at an expression easier to implement than (16), namely,

[TABLE]

II.5 Convergence towards solutions to the Weertman equation

The first author (M.J.) proves in MathWeert that, under our working hypotheses, all the solutions to (6) converge towards the unique solution of (1) at exponential rate. Indeed, consider the following equation:

[TABLE]

where $u_{0}(x)$ is the initial condition of (6). The connection between (18) and (6) resides in that $\varphi(t,x)$ solves (18) if and only if

[TABLE]

solves (6). Now, under mild requirements similar to those of Sec. II.2, Equation (18) can be shown MathWeert to have a unique solution with the following property: if $(\eta,c_{\eta})$ is the solution to (1) with same boundary conditions at infinity as $u_{0}$, then there exist constants $\kappa>0$, $K>0$ and $\xi\in\mathbb{R}$ such that

[TABLE]

The proof relies on a comparison principle, which is a generic property of operators $\partial_{t}-D+F_{\sigma}^{\prime}\left[\cdot\right]$ where $D$ is a dissipative operator (see the classical references Fife ; XinfuChen ).

Combining Equations (19) and (20), one deduces that at large times and uniformly in $x$,

[TABLE]

Thus, given CBCs at infinity, the long-time limit of the solution to (6) is the solution to (1) with same CBCs, up to a time-dependent drift. In the next section, we show how a suitable choice of $c(t)$ eliminates this undesirable effect.

II.6 Centering and choice of $c(t)$

The present section essentially relies on formal arguments. To prescribe $c(t)$ in Equation (6), we need first to center at $x=0$ the core of the solution $\eta(x)$ of (1) by imposing the supplementary condition

[TABLE]

where we have introduced the quantity

[TABLE]

and where the constant $L>0$ represents the half-size of the computational box in the numerical calculations. Imposing condition (22) makes the solution of (1) unique (and not only unique up to translations), which is crucial for numerical purposes.

Next, on the basis of expression (17) for $c_{\eta}$ , we prescribe the function $c(t)$ as

[TABLE]

and $\kappa>0$ is a fixed parameter (the reciprocal of some characteristic time). Since by (21) $u(t,x)$ in (6) converges to $\eta(x)$ up to a translation, we formally have by comparing the following expression to (17):

[TABLE]

Substituting this limit into definition (24) one deduces that at large times

[TABLE]

If we can show that $I(t)\to 0$ when $t\to\infty$, which is equivalent to

[TABLE]

by definition of $I(t)$, then in the same limit we shall have by (26) $c(t)\to c_{\eta}$. The limit (27) indicates that the choice (24) forces the dynamical solution $u(t,x)$ to obey asymptotically the same centering as $\eta(x)$. Put differently, this amounts to computing $u(t,x)$ in a comoving frame. The rightmost term in (24)1 can thus be called a centering correction to the velocity.

To show that $I(t)\to 0$, we differentiate the large-time expression (21)1 of $u(t,x)$ with respect to time. Further invoking (26), we obtain at large times (the dot denotes a total time derivative)

[TABLE]

Integrating (28) over $x\in[-L,L]$ then yields the approximate first-order differential equation

[TABLE]

where the rightmost expression follows from having neglected the $1/x$ term in the asymptotic expansions (13) of $\eta(x)$. From (13) this is legitimate if $L\gg\max(a_{\rm l},a_{\rm r})/(2\pi^{2})+|\zeta(t)|$. The latter condition is compatible with $L$ being fixed in numerical computations if $\zeta(t)$ tends to a finite value at large times. The latter property follows from a simple self-consistent argument. Indeed, Equation (29) implies that $I(t)\simeq I(T){\rm e}^{-\kappa(t-T)}$ for $t>T$ where $T$ is some time above which (21) holds, so that effectively $I(t)\to 0$. Substituting this expression of $I(t)$ into (26), one deduces an approximate analytical expression for $c(t)$ that tends to $c_{\eta}$. Writing (21)2 in the form $\zeta(t)=\zeta(T)+c_{\eta}(T-t)+\smash{\int_{T}^{t}}c(s)\,{\rm d}s$, and substituting the obtained expression of $c(t)$ finally entails the desired saturation property in the form $\zeta(t\to\infty)\simeq\zeta(T)+I(T)/(\eta_{\rm l}-\eta_{\rm r})$. Unfortunately no quantitative estimate of the latter quantity is available. Yet, on the basis of numerical experiments, it can be expected to be of the order of the (unknown) dislocation core width.

As a final remark, we anticipate by indicating that, later on, $\kappa$ will be taken inversely proportional to the algorithmic time step. Thus, in practice $\kappa\to\infty$ in the limit of continuous times. Because of the $\exp(-\kappa t)$ dependence of $I(t)$, we observe that the centering correction in (24)1 remains well-behaved in this limit.

II.7 An analytical solution

For $|\sigma|<1$ and

[TABLE]

with CBCs

[TABLE]

the dimensionless Equation (1) admits the following analytical solution, easily deduced from the well-known solution Rosakis to the original Weertman equation:

[TABLE]

This solution is symmetric in the sense of Section II.3. We will use this test case in Sec. V as a benchmark. So, when $|\sigma|\to 1$ the core width $a$ and the velocity $c_{\eta}$ both blow up as

[TABLE]

Since $c_{\eta}$ and $a$ are not the physical velocity and core width, this behavior is not the hallmark of a physical pathology of the model. It however implies that the comptational box should be taken wider and wider to contain the core of $\eta$, and that the latter moves with nearly infinite velocity, which has numerical consequences to be examined in Sec. V.1.

In the form (13), the asymptotic behaviors deduced from a direct expansion of (32a)1 read

[TABLE]

where the length scales (12) are $a_{\rm l}=a_{\rm r}=a$. Thus, in this particular example where the solution is symmetric the asymptotic behaviors provide a connection between the core width and the next-to-leading terms in the expansion. This supports the interpretation put forward in Sec. II.3 that in generic asymmetric situations $a_{\rm l,r}$ represent characteristic scales of variation of the dislocation density.

III Building blocks

Our numerical scheme to solve the dynamical system (6) crucially rests on evaluating the action of the operator $|\partial_{x}|$ in the Fourier domain, by means of the Discrete Fourier Transform (DFT). This section explains the underlying spatial discretization procedure, and outlines the key features of the implementation.

III.1 Temporal and spatial discretization

At discrete times $t_{n}=n\Delta t$ with step $\Delta t>0$, we need a suitable representation $U^{n}(x)$ of $u(t_{n},x)$ over the whole $x$-axis. To this aim, we define a computational box $[-L,L]$, discretized into $2m$ elementary intervals of width $h=L/m$. We write the function $U^{n}(x)$ as

[TABLE]

where $V^{n}(x)$ is a time-evolving part, the support of which is contained within the box $[-L,L]$, and $\eta^{\rm{ref}}(x)$ is a fixed reference function that complies with the asymptotic behaviors (9). The latter function is prescribed once for all, and plays the role of boundary conditions for the operators $|\partial_{x}|$ and $\partial_{x}$. Decomposition (35) is motivated by the fact that the operator $|\partial_{x}|$ is non-local, and that the solution $\eta$ of (1) does not vanish at infinity. Thus, the tail contributions of $u(t,x)$ at infinity should be taken into account when computing $|\partial_{x}|u$. They are represented in (35) by those of $\eta^{\rm{ref}}(x)$. The need to properly account for tail contributions will be illustrated by means of numerical examples in Section V.4.

In this article, we take $\eta^{\rm{ref}}$ as the linear combination

[TABLE]

and where $a^{\rm ref}$ is an extra arbitrary scaling parameter. The four coefficients $A_{\alpha}$ are determined so as to satisfy the four constraints expressed by (9). Comparing the asymptotic expansions of $\eta^{\rm{ref}}(x)$ for $x\to\pm\infty$ directly deduced from (36) to the generic expressions (9) leads to

[TABLE]

where $\overline{a}$ and $\overline{\eta}$ and have been defined, respectively, in Equations (14) and (23). Thus, choosing $a^{\rm ref}=\overline{a}$ makes $A_{3}$ vanish. In the exact case of Sec. II.7, where furthermore $F^{\prime\prime}_{\sigma}(\eta_{\rm l})=F^{\prime\prime}_{\sigma}(\eta_{\rm r})$, only $A_{1}$ and $A_{2}$ are nonzero and the solution is already completely retrieved at the level of $\eta^{\rm{ref}}(x)$. Therefore we shall need to take $a^{\rm ref}$ different from $\overline{a}$ to be able to use the exact solution of Sec. II.7 as a non-trivial benchmark of the algorithm in Sec. V. We observe that the asymptotic expansion of $\eta^{\rm{ref}}(x)$ in (36) involves only odd inverse powers of $|x|$.

Finally, introducing discrete positions $x_{j}=j\,h$, the function $U^{n}(x)$ is represented inside the box by the vector $\mathbf{u}^{n}\in\mathbb{R}^{2m}$ of components $u^{n}_{j}=U^{n}(x_{j})$. It is decomposed according to (35) as

[TABLE]

where the vector $\mathbf{v}^{n}$ of components $v^{n}_{j}$ is now the unknown of the problem.

It must be emphasized that, however convenient, representation (36) is largely arbitrary. Indeed, any other smooth function $\eta^{\rm{ref}}(x)$ obeying the asymptotic conditions (9), mostly varying inside the computational box, and such that $|\partial_{x}|\eta^{\rm{ref}}(x)$ can be accurately computed once for all (either analytically, or even numerically), would equally well fit our purpose.

III.2 Numerical Fourier transform and discretization of the diffusion operator

This section addresses the discretization of the linear integro-differential operator $|\partial_{x}|$. In view of (2), the latter involves a convolution by a pseudofunction KANW04 , which can be done in the Fourier representation. Crucially, the present approach uses the continuous Fourier form of the operator because this representation is straightforwardly diagonal, and versatile in the sense that it could as well be employed to address the fractional Laplacian $|\partial_{x}|^{\alpha}$.

The FT will be implemented in DFT form, which allows one to benefit from Fast-Fourier-Transform (FFT) routines. As is well known, carrying out a convolution over the real axis by means of DFT and multiplications in the Fourier domain turns this convolution into a periodic one, which induces undesirable periodicity effects in the direct space near both extremities of the interval of interest Bracewell .

The usual workaround is to use zero-padding. Briefly, the technique consists in doubling the spatial extent of the interval of interest, in such a way that the undesirable effects that take place during convolution remain confined to the extra region. The vector to be convolved is continued by zero in the latter region. Once the convolution with the kernel has been carried out in the Fourier space by Fourier component multiplication, and the result has subsequently been transformed back to direct space by DFT inversion, the irrelevant contributions of the added region that have been spoiled by periodicity effects are dropped. This last step constitutes a projection of the result onto the initial interval of interest.

For details the reader is referred to the classical reference (NumRecipes, , p. 643). However, the latter reference only addresses situations where the initial vector and the convolution kernel are both known in the direct space. In contrast, since the present approach uses the convolution kernel in continuous Fourier form, it requires identifying the Fourier modes $k$ in the continuum to the discrete modes $k_{p}$ used by the DFT. For that reason, the labelling and ordering of the vector components hereafter somewhat differ from the ones in Ref. NumRecipes . However, as Equation (46) below shows, the final outcome will take the form of a standard convolution in the direct space.

Thus, zero-padding is carried out by means of the injection $\mathcal{I}$ that maps $\mathbb{R}^{2m}$ into $\mathbb{R}^{4m}$, and the inverse transform involves a projector $\mathcal{P}$ that maps $\mathbb{R}^{4m}$ onto $\mathbb{R}^{2m}$. These operators are defined by

[TABLE]

Note that $\mathcal{P}\mathcal{I}$ leaves $\mathbb{R}^{2m}$ invariant.

The DFT that operates on extended vectors $\mathbf{v}\in\mathbb{R}^{4m}$ is denoted hereafter by $\mathcal{F}_{d}\left\{\mathbf{v}\right\}$ (or, alternatively, by $\mathcal{F}_{d}\left\{v_{p}\right\}$). The corresponding enlarged spatial grid that extends the one in Equation (38) is $x_{j}$ for $j\in\mathcal{K}:=\{-2m,\cdots,2m-1\}$. It admits the dual wavemode grid

[TABLE]

Accordingly, DFTs are carried out for each wave mode $p$ as

[TABLE]

From (35) we can now compute the discretized form of $|\partial_{x}|U(x_{j})$ as

[TABLE]

in which $|\partial_{x}|\eta^{\rm{ref}}$ is evaluated at point $x_{j}$ by means of (36), and $|\partial_{x}|V$ is computed via the DFT as

[TABLE]

Combining (43) and (44) we shall introduce $\left|D_{x}\right|\mathbf{u}$, the discretization of $|\partial_{x}|U$, defined as

[TABLE]

A straightforward check shows that the rightmost term simplifies as

[TABLE]

Although there are $4m$ wave modes, the sum only involves indices $l\in\{-m,\cdots,m-1\}$, precisely because of the use of injection $\mathcal{I}$. As desired, the convolution in (46) is well-behaved and not spoiled by undesirable periodic effects, since it distinguishes at each point $x_{j}$ the contributions $l<j$ of the $u_{l}$ on the left of $x_{j}$, from those on the right for $l>j$. This would not have been the case, had fewer points been retained to compute DFTs.

III.3 Discretization of the advection operator

We discretize the advection operator $c\partial_{x}$ with the following upwind scheme Fletcher :

[TABLE]

where $c^{+}:=\max\left(c,0\right)$ and $c^{-}:=\min\left(c,0\right)$, and where the operators $D_{\pm}$ are implemented via expressions with ${\rm O}(h^{3})$ error (see (Fletcher, , p. 297) with $q=1/2$), namely,

[TABLE]

The choice a third-order implementation is motivated in Sec. V.4 by numerical considerations (see Ref. Fletcher for schemes of orders $1$ and $2$ and Ref. (Hildebrand, , p. 111) for order $4$). Since

[TABLE]

the discretization (47) of the advection operator actually reads

[TABLE]

where $D_{x}[c]$ only depends on the sign of $c$. In this expression, $\eta^{\rm{ref}}$ plays the role of boundary conditions when computing $\partial_{x}U$ at the extremities of the computational box, e.g., for $l\in\{-m,-m+1,m-2,m-1\}$ in (50).

Note that whereas, in principle, implementing $D_{\pm}$ should require adding only a few extra points at the exterior of the interval $[-L,L]$, our implementation (50) is carried out in practice on $4m$ points, due to the introduction of $\mathcal{I}$ and $\mathcal{P}$ for consistency with the implementation of the operator $|\partial_{x}|$. Also, remark that the operators $D_{\pm}$ are almost diagonal, in the sense that their periodizations are diagonal in the Fourier space. In this connection, we introduce for further use in Section IV.2 the vectors $\mathbf{D}_{\pm}\in\mathbb{R}^{4m}$ defined as

[TABLE]

with $\mathbf{e}_{0}=(1,0,\cdots,0)\in\mathbb{R}^{4m}$, and where the periodic convention $0=4m$ applies in the evaluation of the operator $D_{\pm}$ in (51).

III.4 Velocity computation

As seen above, Equation (6) can be solved in the comoving frame by using the velocity $c(t)$ given by (24). In this way, the dislocation core lies as remote as possible from the box boundaries to minimize the influence of the approximations made in handling the tails. To proceed, we substitute in (50), at each time $t_{n}$, the quantity $c$ by $c_{n}=c(t_{n})$ (with the convention that $c_{-1}=0$) computed from a discretized version over $2m-1$ points of expression (24), in which $L$ has been replaced by $L-h$ (since $x_{-(m-1)}=-L+h$ and $x_{m-1}=L-h$); namely,

[TABLE]

In these expressions $\eta_{{\rm l\,r},1}=(\eta_{l}-\eta_{r})a_{\rm l,r}/(2\pi^{2})$, and $\omega_{j}$ are the weights $1/3$, $4/3$, $2/3$, …, $2/3$, $4/3$, $1/3$ of the Simpson integration rule. Remark that $c_{n}$ depends only on $c_{n-1}$ via its sign, which is of little consequence except in calculations at $\sigma$ small where $c_{n}$ is close to [math] and may oscillate during iterations. The term within brackets results from a straightforward discretization of the integral in (24)1, in which the tail contributions have been evaluated analytically from the asymptotic expansions of $\eta^{\rm{ref}}$ (for simplicity, we have not evaluated the full tail contributions of $\eta^{\rm{ref}}$).

Moreover, a suitable value of $\kappa$ stems from the empirical consideration that if we discretize Equation (29) in explicit Euler form as $I_{n+1}=(1-\kappa\Delta t)I_{n}$, monotone convergence towards 0 of $I_{n}$ is ensured by taking $\kappa<1/\Delta t$. Correspondingly, the value of $\kappa$ used henceforth in (52a) is $\kappa=1/(2\Delta t)$.

IV Algorithm

This section describes the iterative numerical scheme used to compute $\bm{\eta}$ and $c_{\bm{\eta}}$. This algorithm is explicit in time, is consistent with the dynamical system (6), combines the above-discretized operators, and remains stable even with a large time step $\Delta t$.

IV.1 Procedure

Computations go as follows. First, for given local minimizers $\eta_{\rm l}$ and $\eta_{\rm r}$ of $F_{\sigma}$, the algorithm is typically initialized by choosing arbitrarily the values $u^{0}_{j}$ inside the box, preferentially not too far from the expected solution. This can be done in a number of ways; notably, by using as an initial condition the function $\smash{\eta^{\rm{ref}}}$ with $\smash{a^{\rm ref}}$ chosen large enough to encompass the expected overall width of the dislocation density (typically, a few times the characteristic scale $\overline{a}$), whence $u^{0}_{j}\equiv 0$. A few low-resolution runs may help adjusting $\smash{a^{\rm ref}}$. Obviously, to get a reasonable representation of the solution, the discretization step $h$ must necessarily be less than $\overline{a}$. Also, when carrying out incremental parametric studies, the solution computed from the previous value of the parameter under consideration can be used as an initial condition, to save CPU time. However, for studying stability issues, we shall purposely take initial data far from the expected shape of a dislocation.

Denoting by $\Phi$ the scheme presented below, we iterate

[TABLE]

until the difference between the results of two successive iterations is small, in the sense that

[TABLE]

where $\Delta_{0}$ is a user-defined stopping criterion. Upon completion at some $n$, the algorithm returns $\bm{\eta}:=\mathbf{u}_{n}$ and the associated $c_{\bm{\eta}}=c_{n}$, evaluated thanks to (52a). Unless otherwise stated, $\Delta_{0}=10^{-10}$ in the numerical calculations of Sec. V below.

IV.2 The Preconditioned Collocation Scheme

The scheme $\Phi$ described hereafter, which we call the Preconditioned Collocation Scheme (PCS), is based on the requirement that the long-time limit $\bm{\eta}$ of $\mathbf{u}^{n}$ should solve the following static equation:

[TABLE]

which is the discretized form of (1). This is a collocation method. A first naive way to proceed would be to attempt solving (6) by writing

[TABLE]

where $\Delta t$ should be adjusted in order to achieve convergence. At convergence, the solution $\bm{\eta}$ obeys (55) and depends on $\Delta t$ only through the evaluation of the numerical velocity $c_{n}$ defined by (52a). As will be justified in Section V.3 this dependence is of little relevance, which is an advantage of this approach.

However, system (56) is ill-conditioned, because it involves the stiff operators $-\left|D_{x}\right|$ and $D_{x}[c_{n}]$. As a consequence, (56) does not converge if $\Delta t$ is not small enough. This issue is dealt with by preconditioning (56) in the following way:

[TABLE]

With the PCS (57) the long-time limit of $\mathbf{u}^{n}$ also satisfies Equation (55). Therefore, the preconditioning achieved by introducing $M(\Delta t)$ is just a procedure to enable and speed up convergence.

Moreover, if we ignore the operators $\mathcal{I}$ and $\mathcal{P}$ in (57b), the matrix $M(\Delta t)$ is close to the effective inverse of $1+\Delta t\left|D_{x}\right|-\Delta tc_{n}D_{x}[c_{n}]$. Hence, (57a) preconditioned by (57b) is nothing but the following semi-implicit scheme (see (HairerWanner2, , p. 102) for a reference on semi-implicit schemes):

[TABLE]

As is well-known, treating stiff operators in an implicit way yields a stable scheme. Hence, this preconditioning naturally leads to stability (this assertion will be exemplified in Section V). We note that (58) is a consistent discretization of (6) —just as (56).

We now justify the preconditioned character of $\Phi$ by analogy with the task of solving iteratively a linear system. Indeed, ignoring nonlinearities by replacing $F_{\sigma}^{\prime}\left(\mathbf{u}^{n}\right)$ by a constant $\mathbf{b}$ and by setting $c_{n}=0$, (55) can be put in the form $\mathsf{A}\mathbf{u}=\mathbf{b}$, so that (56) reduces to a scheme of the type

[TABLE]

where $\mathsf{A}$ is a positive symmetric matrix. The latter scheme converges for general $\mathbf{b}$ if and only if all eigenvalues of $\left(1-\Delta t\,\mathsf{A}\right)$ belong to the interval $(-1,1)$; this can require $\Delta t$ to be very small. In a similar way, the scheme (57) can be abstracted as

[TABLE]

with $\mathsf{M}$ close to $\left(1+\Delta t\mathsf{A}\right)^{-1}$. Then, the eigenvalues of the latter scheme are close to those of $(1+\Delta t\,\mathsf{A})^{-1}$, which unconditionally belong to $(0,1)$ if $\Delta t>0$. Equation (60) amounts to solving $\mathsf{M}\mathsf{A}\mathbf{u}=\mathsf{M}\mathbf{b}$, which is the classical preconditioning method (in the spirit of a modified Richardson iteration, see (Kelley, , p. 6)).

V Numerical results

The above algorithm has been implemented in the form of a Matlab® code, and the following results have been obtained on a 2.3GHz standard laptop computer, with typical computation times of order one second to a few minutes per run, depending on the case at hand. Except in Sec. V.5, the calculations concern benchmark comparisons with the exact solution of Sec. II.7, and have been carried out with $F^{\prime}_{\sigma}$ as defined by (30), for $\sigma\in[0,1)$. Moreover, in Sections V.1 to V.4, a parameter value $a^{\rm ref}=\overline{a}/2$ (see Section III.1) is employed in $\smash{\eta^{\rm{ref}}(x)}$.

V.1 Convergence

As stated in Sec. II.7, the core width $a(\sigma)$ and the velocity $c_{\eta}(\sigma)$ blow up when $\sigma\rightarrow 1$. The first problem is easily solved by running one preliminary low-resolution run to provide a rough numerical estimate for the core width $\widetilde{a}(\sigma)$. The half computational box size $L$ is then adjusted to a value $L\gg\widetilde{a}(\sigma)$. Figure 3 displays the initial data $\mathbf{u}^{0}$, the converged result $\bm{\eta}$ and its discrete derivative. The figure illustrates the robustness of the PCS with respect to initial conditions in the sense that the initial data $\mathbf{u}^{0}$ can be non-monotone, irregular, and far from the solution $\eta$ to (1). In this calculation, the applied loading is $\sigma=1-1.973\times 10^{-3}$, which induces large values for the converged velocity and core width, $c_{\eta}(\sigma)\simeq a(\sigma)\simeq 15.9$, close to one another; see Equation (33). The half box size is $L=5\,a(\sigma)\simeq 80$. As seen in Fig. 3(c) the discretized derivative of $\bm{\eta}$ is regular inside the box $[-L,L]$, but the inset shows that some artifacts take place near the matching points between the solution inside the box and the tails of $\eta^{\rm{ref}}(x)$ outside it. Quite generally, it is observed that these artifacts diminish when either $L$ or $\sigma$ are increased (results not shown).

Under the same initial conditions, Figure 4 illustrates the convergence properties with time, via $\left\|\Delta\mathbf{u}^{n}\right\|$ defined in (54). The semi-logarithmic plot of Fig. 4(a) shows that, up to high-frequency oscillations, convergence is exponential with time, in agreement with the arguments of Section II.5. Moreover, the convergence rate (the slopes in Fig. 4(a)) depends on the applied loading $\sigma$, convergence being impeded when $\sigma$ approaches $1$. Not unexpectedly, this loss of performance coincides with $F_{\sigma}$ being ‘less and less bistable’ in the sense that the minima of $F_{\sigma}$ at $\eta=\eta_{\rm r}$ and $\eta=\eta_{\rm l}$ become less and less deep. Bistability is a crucial requirement for the existence of a solution to (1), and for proving convergence as expressed by (20). In this connection, it should be remarked that the tail expressions (9) involve the second derivatives $F^{\prime\prime}_{\sigma}(\eta)$, which presumably leads to pathologies when the latter are small. However, Fig. 4(b), which displays the data in parametric form of parameter $\sigma$, shows that the PCS copes well with high velocities, which is important from a physical standpoint. Thus, our preconditioning does a nice job of avoiding stability issues when $c_{\eta}$ is high.

V.2 Error indicators and overall accuracy

Fig. 5 compares the output of the PCS with the exact solution of Sec. II.7 in terms of the following indicators:

[TABLE]

The absolute error (61a) can as well be understood as a relative error since in the cases considered $\left\|\eta\right\|_{{\rm L}^{\infty}(\mathbb{R})}$ is of order $1$. In (61a) and (61b) the $\inf_{\xi}$ operation, motivated by the translation invariance of (1), takes care of the approximate character of the centering of the numerical solution, due to discretization errors. Formulated in this way, the error indicators are insensitive to small shifts in the position of the computed solution. The figure indicates that the outputs of the PCS accurately approximate the exact results, with errors of a similar order of magnitude for the three quantities represented. In addition, the small errors observed on Figure 5 depend only weakly on $c_{\eta}$.

V.3 Discretization parameters and error scaling

This section closely investigates the scalings of the error at convergence with respect to the time step $\Delta t$, the half box size $L$ and the space discretization step $h$. The same case as above is considered, with a variety of applied loadings.

Addressing first the influence of $\Delta t$, we observe that the PCS solution depends on $\Delta t$ only through the centering correction term in expression (52a) of $c_{n}$. However, as explained in Section III.4, this term vanishes in the limit of infinite times, and is therefore expected to be small at convergence. This is confirmed by the errors reported in Table 1 for applied loading $\sigma=\sigma_{1}$ (see caption) so that $c_{\eta}(\sigma_{1})\simeq 7.91$, and decreasing values of $\Delta t$. Errors are quasi-constant, which shows that the dependence on $\Delta t$ of the converged result is negligible. For definiteness, the rest of the calculations in the present paragraph is made with $\Delta t=0.1$.

Figure 6 illustrates how $\left\|\Delta\eta\right\|$ depends on $L$ and $h=L/m$, for the two contrasted velocities obtained with loadings $\sigma_{2}=0.951$ so that $c_{\eta}(\sigma_{2})\simeq 3.08$, and $\sigma_{3}=1-1.97\times 10^{-5}$ so that $c_{\eta}(\sigma_{3})\simeq 159$. Raw results are displayed in Figs. 6(a) and (b). Computations for $\sigma=\sigma_{3}$ did not converge with $h/a=0.31$, which is why this value is not considered in plots (b) and (d). It turns out that for $L$ large, $\left\|\Delta\eta\right\|$ scales as $h^{3}$ for both loadings. Besides, it is approximately proportional to $L^{-3}$ for $L$ small. This is demonstrated in the data-collapse plots (c) and (d) where the individual datasets of Fig. 6(a) and (b), respectively, are merged into one single master curve by means of appropriate rescalings of abscissas and ordinates. The data collapse of Fig. 6(b) for $\sigma=\sigma_{3}$ is only partially successful, as noticeable corrections to scaling arise for $Lh\ll a(\sigma_{3})^{2}$.

The scalings can be understood as follows. On the one hand, the scaling in $h$ is consistent with our choice of a third-order advection scheme in (48). On the other hand, since the error at any point is spread over the whole domain by the integro-differential operator $|\partial_{x}|$, the error $\eta(L)-\eta^{\rm{ref}}(L)$ at the boundary point $x=L$ (for instance) can be used to estimate the overall error. It behaves as $L^{-3}$ in the present case where $\eta(x)$ and $\eta^{\rm{ref}}(x)$ are symmetric, and where the next nonzero term in expansion (13) is $\propto x^{-3}$. This is consistent with the error scaling observed in the plots. Still, these elementary arguments do not explain the downwards bending of the high-velocity plots in Fig. 6(b), which indicates either that the $L^{-3}$ scaling regime has not been reached, or that it may not hold exactly. This bending causes deviations from ideal scaling, made conspicuous in Fig. 6(d).

Likewise, Figures 7 and 8 display the errors (61b) and (61c), respectively, for the loadings $\sigma_{2}$ and $\sigma_{3}$. Whereas Fig. 7(a) resembles Fig. 6(a), the dependence of the error with respect to $h$ is more involved. In particular for $\sigma=\sigma_{2}$ (Fig. 7(a)), there exists at fixed $L$ an optimal value of $h=L/m$ that minimizes the error. The reason for this behavior is unclear. However it should be realized that approximating $\partial_{x}\eta$ in the sense of the ${\rm L}^{\infty}$ norm is quite demanding. Figs. 7(b) and 8(b) display bendings similar as in Fig. 6(b). No data collapse is presented, as notable deviations from scaling take place in all figures.

For practical matters one needs, e.g., to determine the optimal value of $L$ at fixed number $2m$ of discretization points, which somehow corresponds to a fixed cost to go from $\bm{u}^{n}$ to $\bm{u}^{n+1}$ in (57a). Table 2 shows the minimal error $\left\|\Delta\eta\right\|$ deduced from the datasets of Fig. 6, together with the corresponding optimal value of $L$ via the ratio $(hL)/a^{2}=(L/a)^{2}/m$. One observes that the optimal value of the latter quantity does not depend on $m$ for $m$ sufficiently large. This can be understood from the above scaling arguments. Indeed, the data collapse in Fig. 6(c) indicates that for $Lh\ll C$, where $C$ is some constant, $\left\|\Delta\eta\right\|\propto L^{-3}$, and that for $Lh\gg C$, $\left\|\Delta\eta\right\|\propto h^{3}=L^{3}/m^{3}$. Thus the optimum takes place at the crossover between these regimes in which the error first decreases, then increases with $L$. Balancing the two regimes, it follows that $Lh=L^{2}/m\simeq C$ at the optimum. The analysis still qualitatively holds for higher $\sigma$ values. The optimum just discussed is not the one evoked above in connection with Fig. 7. However, in both cases, one sees that increasing the number of discretization points at fixed $L$ does not necessarily improve the accuracy: one also needs to increase $L$.

Same arguments as above suggests that, in the general case of a nonsymmetric solution where the next nonzero term in expansion (13) can be an ${\rm O}(x^{-2})$ (see Sec. II.3), the error on $\eta$ would behave as $L^{-2}$ instead of $L^{-3}$, and if the error scales as $h^{3}$ for $h$ small the optimum would take place for $Lh^{3/2}\simeq C$.

V.4 Influence of the tails and of the order of the advection scheme

At given $L$ the quality of the numerical solution also depends on how tail contributions are accounted for while handling the operator $|\partial_{x}|$. As the quantity $|\partial_{x}|\eta(x)$ enters the equation for $\eta$, errors on the former affect the latter via the nonlinear term. To illustrate this point, we focus on the point $x=0$ where $\eta(0)=\eta^{\rm{ref}}(0)=\overline{\eta}$. Since both the function $\eta^{\rm{ref}}$ and the solution $\eta$ to (1) satisfy the same asymptotic behavior (9), their difference behaves as $\eta(x)-\eta^{\rm{ref}}(x)=O\left(|x|^{-\tau}\right)$, where $\tau=2$ in the general case, and $\tau=3$ in the present, symmetric, benchmark case. Hence by (3), the long-range contribution to the error in $|\partial_{x}|\eta(0)$ due to the replacement of $\eta$ by $\eta^{\rm{ref}}$ outside the computational box is bounded by

[TABLE]

Thus, in the results of the previous section V.3 for which $\tau=3$ the observed scaling stems from local errors at boundary points, which scale as $L^{-3}$, and not from the present long-range contributions of errors in tails, which scale at most as $L^{-4}$. In the general non-symmetric case where $\tau=2$, the scalings of these errors are presumably changed into $L^{-2}$ (see end remark in the previous section) and $L^{-3}$, respectively. Therefore, using a function $\eta^{\rm{ref}}(x)$ with faithful tails (in the sense of Section III.1) should make local errors dominant over tail errors in all cases. In contrast, upon not using such a reference function, or upon periodizing the calculation (thus, disregarding tails in both cases), the error on $|\partial_{x}|\eta(0)$ would instead scale as

[TABLE]

We indeed obtained such a scaling while carrying out some preliminary studies (results not shown). To summarize, disregarding tails deteriorates the accuracy of the calculations. On the contrary, appealing to the zero-padding trick and using a reference function $\eta^{\rm{ref}}$ with faithful tails to handle properly the operator $|\partial_{x}|$ strongly improves it.

The overall error also depends on the order in $h$ of the discrete advection operator $D_{\pm}$ of Sec. III.3. Along with the calculations of the previous Section with an advection scheme of order 3, we also carried out similar calculations with schemes of order 2 and 4. In both latter cases, the power of the scaling in $h$ was found identical to the order of the scheme for $h$ small enough (results not shown). Thus, at least for orders 2 to 4, the scaling in $h$ is closely related to the order or the discretization scheme of the advection operator, the discretization errors involved by the DFT in the computation of the operator $|\partial_{x}|$ presumably being subdominant.

Figures 9 further illustrate these points by means of the indicators (61a) and (61c). The plots have been made with advection upwind schemes of order $1$, $2$, $3$ and $4$, and two different types of asymptotes; namely, ‘refined tails’ (r.t.) (9), and ‘constant tails’ (c.t.) in which the inverse power-law correction in (9) is dropped. Both figures show that using refined tails provides in most cases orders-of-magnitude gains of accuracy over using constant tails. The dependence of the error on the order of the upwind scheme is more difficult to interpret, although orders 3 and 4 lead to better accuracy. As the second-order scheme behaves irregularly in Fig. 9 (b), the third-order scheme is used in the rest of the article.

V.5 A generalized example: the camel-hump potential

We now evaluate our method with the ‘camel-hump’ type potential $F_{\sigma}$ defined by (8). Depending on $\sigma$, and on the parameter $r>0$, $F_{\sigma}$ can feature between its two main minima at $\eta_{\rm l}$ and $\eta_{\rm r}$, and its humps, an intermediate local minimum of depth controlled by $r$, leading thus to a more or less dissociated solution $\eta$ to (1). The potential (8) has been derived by means of Lejček’s method Lejcek75 so as to provide the following exact dissociated dislocation solution to the PN equation when $\sigma=0$:

[TABLE]

(note that for $r=0$, the camel-hump potential reduces to the simple sinusoidal form (30)).

In this section, calculations have been carried out with $a^{\rm ref}=5\overline{a}$, because the overall width of the dislocation is much larger than that of the individual peaks in the density. Figure 10 shows that the algorithm correctly recovers this solution. Due to the presence of a secondary hollow in the potential, the method takes longer to converge.

We finally present an application to a more physically relevant case, for which the exact solution is unknown. Indeed, when $\sigma>0$ and $r>0$ no analytical solution to (1) is available. However, as discussed in Sec. II.2, there exists a solution to (1) that can be computed with our algorithm for any $\sigma\in(0,\sigma_{\text{lim}})$, where $\sigma_{\text{lim}}=\max_{\eta}F^{\prime}_{\sigma}(\eta)$. Figure 11 displays the solution $\bm{\eta}$ (a) and its derivative (b), which is the dislocation density. In this example, the latter features two bumps that represent partial dislocations. The midpoint between the two partial dislocations corresponds to the local minimum of $F_{\sigma}$ in $\eta_{\rm m}\in\left(\eta_{\rm r},\eta_{\rm l}\right)$. The asymmetry of the dislocation density can be interpreted as a consequence of the nonzero driving force $\sigma$ coupled to the dissociation process induced by the camel-hump character of the potential. Figure 12 displays two quantities of physical interest, namely, the velocity $c_{\eta}(\sigma)$ and the effective core width $a(\sigma)$, for $\sigma$ varying between [math] and $\sigma_{\text{lim}}=\max\left(F_{0}^{\prime}\right)=0.5902$. Here, the quantity $a(\sigma)$ has been computed by minimizing the ${\rm{L}}^{2}$ norm of the difference between the numerical solution $\bm{\eta}$ and the function, parametrized by $a$ and $x_{0}$:

[TABLE]

As expected from an analogy with the simpler case of Section II.7, both $a(\sigma)$ and $c_{\eta}(\sigma)$ increase with $\sigma$ and blow up when $\sigma\rightarrow\sigma_{\text{lim}}$.

VI Concluding discussion

To summarize, we have proposed the Preconditioned Collocation Scheme (PCS), which is a numerical procedure to approximate solutions to (1), based on the dynamical system (6). The PCS uses the continuous FT of the operator $|\partial_{x}|$ and takes advantage of the FFT in its implementation, in spite of the strong constraint that the desired solution is not periodic, but has boundary conditions at infinity. We have taken advantage from the exact asymptotic expansion of the solution to improve the accuracy of the numerical approximation. Also, we have shown that an overall ${\rm O}(h^{3})$ error in the space discretization could be achieved by means of a third-order advection scheme.

The method employed remains stable when the (a priori unknown) advection part scaled by the velocity $c_{\eta}$ dominates over the diffusion part in (1), which allows to investigate the asymptotic behavior of (1) when $c_{\eta}$ is large. The PCS solves a discretized version of (1). Being preconditioned, it can be used with a large time step $\Delta t$, nonetheless delivering outputs that depend very weakly on $\Delta t$. Although this was not illustrated, we add that if $\eta_{\rm l}$ and $\eta_{\rm r}$ are not exactly computed as exact local minimizers of $F_{\sigma}$ (e.g., if they suffer from slight numerical errors), the algorithm converges as well. Our method has some limitations, however: the more advection dominates diffusion, the more iterations the method requires to converge.

Still, the time and space complexities of the algorithm give satisfactory accuracy at reasonable computational cost. Indeed, the PCS requires ${\rm O}(m)$ memory space. Moreover DFTs have been speeded up by means of Fast Fourier Transform routines. Therefore, each iteration step takes ${\rm O}(m\log m)$ CPU time. Given that the number of iterations to convergence obviously scales as $1/\Delta t$, the PCS therefore has an overall time complexity of order ${\rm O}(m\log m/\Delta t)$. Then, achieving $10^{-4}$ accuracy on $\eta$ and $c_{\eta}$ with 2048 discretization points takes about one second on a standard laptop with CPU running at 2.3GHz, except in difficult cases when $c_{\eta}\rightarrow+\infty$, where it is slower. This study involved runs with up to $10^{5}$ discretization points.

We point out that alternative schemes could be used to simulate the dynamical system (6). First, other time integrators can be appealed to. For instance, classical splitting schemes such as Strang or Lie splittings Hairer2010 aim at simulating a dynamical system that involves a sum of operators, and turn the latter into a composition of evolution operators. Although these investigations were not reported for conciseness, we have checked that such methods do indeed apply to our problem, and prove stable and robust. However, the solution they produce is inherently $\Delta t$-dependent. Hence, achieving high accuracy requires small $\Delta t$ values, which makes them expensive. Also, as mentioned in Section IV.2, an explicit Euler scheme could as well be employed. However, the latter approach proves unstable if $\Delta t/h$ is not small. In contrast, the alternative semi-implicit scheme sketched in (58) is stable whatever $\Delta t$. However, it is costly as it requires the inversion of the operator $1+\Delta t\left(\left|D_{x}\right|-c_{n}D_{x}[c_{n}]\right)$. To summarize, among all the schemes we have investigated, we deem the PCS the most robust and least expensive one. Second, we have deliberately chosen to implement the operator $|\partial_{x}|$ in continuous Fourier form. It would be equally possible to discretize the integral representation of the operator, e.g., in the convenient form (3). However, this may require adapted integration rules Movchan . In this respect, a method taking advantage of FTs proves more straightforward. It also proves more versatile because the analytical form of the kernel under consideration might be unavailable in cases involving a different integrodifferential operator. The main constraint for diagonalization by the FT is that the operator be translation-invariant.

As a perspective, although definite conclusions about the validity of the approach for more general equations are yet to be obtained, we have all reasons to believe that this method applies as well to equations of the type Gui

[TABLE]

where $F$ is a bistable nonlinearity, $\alpha\in(0,2]$, and $|\partial_{x}|^{\alpha}$ is the operator of Fourier symbol $|k|^{\alpha}$. However, we suspect that the preconditioning we use is insufficient for ensuring unconditional stability if $\alpha$ is larger than some value $\alpha_{0}\in(1,2)$. In fact, we have explored the classical advection-reaction-diffusion case (where $\alpha=2$, and $|\partial_{x}|^{2}=-\Delta$, see Appendix A), with

[TABLE]

where $c_{\eta}$ is the same as in the left-hand side of Equation (66). In this case, (66) admits the analytical solution $\eta(x)=-\tanh(x)$ (Griffiths, , p. 291). Employing a variant of the above-described method, we have recovered a numerical approximation of this analytical solution. Our method presumably applies as well to the modified Weertman equation with gradient term Rosakis , in which the operator $|\partial_{x}|$ in (1) is replaced by $|\partial_{x}|-\lambda\Delta$, where $\lambda>0$.

Acknowledgements.

M. Josien thanks for its hospitality the CEA-DAM Île-de-France where part of this work was carried out.

Appendix A Case of the Laplacian

We briefly justify here that, with little modification, the PCS presented here to address (6) can also be used for the classical reaction-advection-diffusion equation

[TABLE]

in which the operator $-\partial_{xx}$ replaces $|\partial_{x}|$. Then, the counterpart of (45) reads

[TABLE]

with

[TABLE]

Thus, in a first approach, one would only need to replace in (57a) the operator $\left|D_{x}\right|$ by $-D_{xx}$. The preconditioning is adapted as follows:

[TABLE]

This modified scheme solves a dynamical equation associated with (68), similar in spirit to (we do not write it down for conciseness). However, it turns out that this preconditioning does not ensure unconditional stability (namely, if $h$ is small, the method only converges for $\Delta t$ small). This instability is due to the truncation of the varying part $v$ of $u$ outside the box $[-L,L]$, which induces a numerical singularity in the discretization of the Laplacian near boundaries. In contrast, the preconditioning associated with $|\partial_{x}|$ in Section IV.2 suffices to damp potential oscillations, probably because the latter operator is not as stiff as the Laplacian. A simple albeit costly way to overcome this stability issue is to take $\Delta t$ small.

A cheaper way is to use the classical centered discretization of the Laplacian and set

[TABLE]

instead of (69). The associated $\kappa_{p}$ that replaces (70) in (71) reads

[TABLE]

In (72), the discretization at boundary points, e.g., at $j=-m$, can be explained by the following heuristic approximation:

[TABLE]

since $v(x)$ is supposed to vanish for $x<-L$. Actually, if we impose $\eta^{\rm{ref}}=0$, $c_{n}=0$ and $F=0$, a numerical study of the eigenvalues of the PCS induced by (74) indicates that the moduli of these eigenvalues lie between [math] and $1$ (although they can be very close to $1$), which implies that the scheme is stable. Then, with $F$ as in (67), the PCS delivers outputs close to the exact solution (results not shown), although no systematic error assessment has been made in this case.