Numerical analysis of nonlocal fracture models in H\"older space

TL;DR

This paper analyzes the convergence rate of finite difference methods for nonlocal fracture models, demonstrating stability and convergence to brittle fracture evolution in H"older spaces.

Contribution

It provides a detailed convergence analysis for finite difference schemes applied to nonlocal fracture models within H"older spaces, including stability results.

Findings

Convergence rate depends on $C_s h^b3/b5^2$ factor.

Both explicit and implicit schemes achieve the convergence rate.

H"older continuous evolutions converge to brittle fracture as nonlocality vanishes.

Abstract

In this work, we calculate the convergence rate of the finite difference approximation for a class of nonlocal fracture models. We consider two point force interactions characterized by a double well potential. We show the existence of a evolving displacement field in H\"{o}lder space with H\"{o}lder exponent $\gamma \in (0,1]$. The rate of convergence of the finite difference approximation depends on the factor $C_s h^\gamma/\epsilon^2$ where $\epsilon$ gives the length scale of nonlocal interaction, $h$ is the discretization length and $C_s$ is the maximum of H\"older norm of the solution and its second derivatives during the evolution. It is shown that the rate of convergence holds for both the forward Euler scheme as well as general single step implicit schemes. A stability result is established for the semi-discrete approximation. The H\"older continuous evolutions are seen to converge to a brittle fracture evolution in the limit of vanishing nonlocality.

Equations

$$S=S(\mathbf{y},\mathbf{x};\mathbf{u})$$

$${\mathbf{W}}^{ϵ}(S,\mathbf{x})=\frac{1}{{ϵ}^d{\omega }_d}{\int }_{H_{ϵ}(\mathbf{x})}|\mathbf{y}-\mathbf{x}|W^{ϵ}(S,\mathbf{y}-\mathbf{x})d\mathbf{y},$$

$$P{D}^{ϵ}(\mathbf{u})$$

$$\rho {\partial }_{tt}^2\mathbf{u}(t,\mathbf{x})$$

$$-\nabla P{D}^{ϵ}(\mathbf{u})(\mathbf{x})=\frac{2}{{ϵ}^d{\omega }_d}{\int }_{H_{ϵ}(\mathbf{x})}{\partial }_S{W}^{ϵ}(S,\mathbf{y}-\mathbf{x})\frac{\mathbf{y}-\mathbf{x}}{\left|\mathbf{y}-\mathbf{x}\right|}d\mathbf{y}$$

$$\mathbf{u}(\mathbf{x})=0\forall \mathbf{x}\in \partial D,$$

$$\mathbf{u}(0,\mathbf{x})={\mathbf{u}}_0(\mathbf{x}){\partial }_t\mathbf{u}(0,\mathbf{x})={\mathbf{v}}_0(\mathbf{x})$$

$${\mathcal{E}}^{ϵ}(\mathbf{u})(t)$$

$$P{D}^{ϵ}(\mathbf{u})$$

$$\frac{d}{dt}{\mathcal{E}}^{ϵ}(\mathbf{u})(t)=(\ddot{\mathbf{u}}(t),\dot{\mathbf{u}}(t))-(-\nabla P{D}^{ϵ}(\mathbf{u}(t)),\dot{\mathbf{u}}(t)),$$

$$W^{ϵ}(S,\mathbf{y}-\mathbf{x})$$

$$\underset{r\to 0^{+}}{\lim }\frac{f(r)}{r}=f^{\prime }(0),\underset{r\to \infty }{\lim }f(r)=f_{\infty }<\infty$$

$$\left|S(\mathbf{y},\mathbf{x};\mathbf{u})\right|>S_c.$$

$$f^{\prime }(r^2)+2r^2{f}^{\prime \prime }(r^2)=0.$$

$${\Vert \mathbf{u}\Vert }_{C^{0,\gamma }(D;{\mathbb{R}}^d)}$$

$${\left[\mathbf{u}\right]}_{C^{0,\gamma }(D;{\mathbb{R}}^d)}$$

$$F_1^{ϵ}(y,t)$$

$$F_2^{ϵ}(y,t)$$

$$\frac{d}{dt}y=F^{ϵ}(y,t),$$

$${\Vert F^{ϵ}(y,t)-F^{ϵ}(z,t)\Vert }_X$$

$$\le \frac{\left(L_1+L_2\left(\Vert \omega {\Vert }_{C^{0,\gamma }(D)}+{\Vert y\Vert }_X+{\Vert z\Vert }_X\right)\right)}{{ϵ}^{2+\alpha (\gamma )}}{\Vert y-z\Vert }_X$$

$$\alpha (\gamma )=\{\begin{array}{ll}0& \text{if }\gamma \ge 1/2\\ 1/2-\gamma & \text{if }\gamma <1/2.\end{array}$$

$${\Vert F^{ϵ}(y,t)\Vert }_X$$

$${\Vert F^{ϵ}(y,t)-F^{ϵ}(z,t)\Vert }_X$$

$$y(t)=x_0+{\int }_0^t{F}^{ϵ}(y(\tau ),\tau )d\tau ,$$

$$y^{\prime }(t)=F^{ϵ}(y(t),t),withy(0)=x_0,$$

$$\frac{{\hat{\mathbf{u}}}_i^{k+1}-{\hat{\mathbf{u}}}_i^k}{\Delta t}$$

$$\frac{{\hat{\mathbf{v}}}_i^{k+1}-{\hat{\mathbf{v}}}_i^k}{\Delta t}$$

$${\hat{\mathbf{u}}}^k(\mathbf{x})$$

$${\hat{\mathbf{v}}}^k(\mathbf{x})$$

$$E^k$$

$$\underset{0\le k\le T/\Delta t}{\sup }{E}^k\le O\left(C_t\Delta t+C_s\frac{h^{\gamma }}{{ϵ}^2}\right),$$

Full text

{AMS}

34A34, 34B10, 74H55, 74S20

1 Introduction

Nonlocal formulations have been proposed to describe the evolution of deformations which exhibit loss of differentiability and continuity, see [28] and [31]. These models are commonly referred to as peridynamic models. The main idea is to define the strain in terms of displacement differences and allow nonlocal interactions between material points. This generalization of strain allows for the participation of a larger class of deformations in the dynamics. Numerical simulations based on peridynamic modeling exhibit formation and evolution of sharp interfaces associated with phase transformation and fracture [9], [32], [27], [16], [1], [11], [24], [6], [20], [30], [34], [17]. A recent summary of the state of the art can be found in [14].

In this work, we provide a numerical analysis for the class of nonlocal models introduced in [22] and [23]. These models are defined by a double well two point potential. Here one potential well is centered at zero and associated with elastic response while the other well is at infinity and associated with surface energy. The rational for studying these models is that they are shown to be well posed over the class of square integrable non-smooth displacements and, in the limit of vanishing non-locality, their dynamics recover features associated with sharp fracture propagation see, [22] and [23]. The numerical simulation of prototypical fracture problems using this model is carried out in [24]. In order to develop an $L^{2}$ approximation theory, we show the nonlocal evolution is well posed over a more regular space of functions. To include displacement fields which have no well-defined derivatives, we consider displacement fields in the Hölder space $C^{0,\gamma}$ with Hölder exponent $\gamma$ taking any value in $(0,1]$. We show that a unique evolution exists in $C^{0,\gamma}$ for $C^{0,\gamma}$ initial data and body force. The semi-discrete approximation to the Hölder continuous evolution is considered and it is shown that at any time its energy is bounded by the initial energy and the work done by the body force. We develop an approximation theory for the forward Euler scheme and show that these ideas can be easily extended to the backward Euler scheme as well other implicit one step time discretization schemes. It is found that the discrete approximation converges to the exact solution in the $L^{2}$ norm uniformly over finite time intervals with the rate of convergence proportional to $(C_{t}\Delta t+C_{s}h^{\gamma}/\epsilon^{2})$, where $\Delta t$ is the size of time step, $h$ is the size of spatial mesh discretization, and $\epsilon$ is the length scale of nonlocal interaction relative to the size of the domain. The constant $C_{t}$ depends on the $L^{2}$ norm of the time derivatives of the solution, $C_{s}$ depends on the Hölder norm of the solution and the Lipschitz constant of peridynamic force. We point out that the constants appearing in the convergence estimates with respect to $h$ can be dependent on the horizon and be large when $\epsilon$ is small. This is discussed in section 3 and an example is provided in section 6. These results show that while errors can grow with each time step they can be controlled over finite times $t$ by suitable spatial temporal mesh refinement. We then apply the methods developed in [22] and [23], to show that in the limit $\epsilon\rightarrow 0$, the Hölder continuous evolutions converge to a limiting sharp fracture evolution with bounded Griffiths fracture energy. Here the limit evolution is differentiable off the crack set and satisfies the linear elastic wave equation.

In the language of nonlocal operators, the integral kernel associated with the nonlocal model studied here is Lipschitz continuous guaranteeing global stability of the finite difference approximation. This is in contrast to PDE based evolutions where stability can be conditional. In addition we examine local stability. Unfortunately the problem is nonlinear so we don’t establish CFL conditions but instead identify a mode of dynamic instability that can arise during the evolution. This type of instability is due to a radial perturbation of the solution and causes error to grow with each time step for the Euler scheme. For implicit schemes this perturbation can become unstable in parts of the computational domain where there is material softening, see subsection 3.4. Of course stability conditions like the CFL conditions for linear nonlocal equations are of importance for guidance in implementations. In the case of $d=1$, a CFL type condition is obtained for the finite difference and finite element approximation of the linear peridynamic equation, see [18]. Recent work develops a new simple CFL condition for one dimensional linearized peridynamics in the absence of body forces [15]. Related analysis for the linear peridynamic equation in one dimension is taken up in [35] and [7]. The recent and related work [12] and [19] addresses numerical approximation for problems of nonlocal diffusion.

There is now a large body of contemporary work addressing the numerical approximation of singular kernels with application to nonlocal diffusion, advection, and mechanics. Numerical formulations and convergence theory for nonlocal $p$-Laplacian formulations are developed in [13], [26]. Numerical analysis of nonlocal steady state diffusion is presented in [33] and [25], and [8]. The use of fractional Sobolev spaces for nonlocal problems is investigated and developed in [12]. Quadrature approximations and stability conditions for linear peridynamics are analyzed in [35] and [29]. The interplay between nonlocal interaction length and grid refinement for linear peridynamic models is presented in [7]. Analysis of adaptive refinement and domain decomposition for linearized peridynamics are provided in [3], [21], and [2]. This list is by no means complete and the literature on numerical methods and analysis continues to grow.

The paper is organized as follows. In section 2, we describe the nonlocal model. In subsection 2.2, we state theorems which show Lipschitz continuity of the nonlocal force (Proposition 2.1) and the existence and uniqueness of an evolution over any finite time interval (Theorem 2.2). In section 3, we compute the convergence rate of the forward Euler scheme as well as implicit one step methods. We identify stability of the semi-discrete approximation with respect to the energy in subsection 3.3. In subsection 3.4, we identify local instabilities in the fully discrete evolution caused by suitable radial perturbations of the solution. In section 4, we give the proof of Proposition 2.1, Theorem 4.1, and Theorem 2.2. The convergence of Hölder continuous evolutions to sharp fracture evolutions as $\epsilon\rightarrow 0$ is shown in section 5. In section 6 we present an example showing the effect of the constants $C_{t}$ and $C_{s}$ on the convergence rate and summarize our results.

2 Double well potential and existence of a solution

In this section, we present the nonlinear nonlocal model. Let $D\subset\mathbb{R}^{d}$, $d=2,3$ be the material domain with characteristic length-scale of unity. Let $\epsilon\in(0,1]$ be the size of horizon across which nonlocal interaction between points takes place. The material point $\boldsymbol{x}\in D$ interacts nonlocally with all material points inside a horizon of length $\epsilon$. Let $H_{\epsilon}(\boldsymbol{x})$ be the ball of radius $\epsilon$ centered at $\boldsymbol{x}$ containing all points $\boldsymbol{y}$ that interact with $\boldsymbol{x}$. After deformation the material point $\boldsymbol{x}$ assumes position $\boldsymbol{z}=\boldsymbol{x}+\boldsymbol{u}(\boldsymbol{x})$. In this treatment we assume infinitesimal displacements and the strain is written in terms of the displacement $\boldsymbol{u}$ as

[TABLE]

Let $W^{\epsilon}(S,\boldsymbol{y}-\boldsymbol{x})$ be the nonlocal potential density per unit length between material point $\boldsymbol{y}$ and $\boldsymbol{x}$. The energy density at $\boldsymbol{x}$ is given by

[TABLE]

where $\omega_{d}$ is the volume of a unit ball in $d$-dimension and $\epsilon^{d}\omega_{d}$ is the volume of the ball of radius $\epsilon$. The potential energy is written as

[TABLE]

and the displacement field satisfies following equation of motion

[TABLE]

for all $\boldsymbol{x}\in D$. Here we have

[TABLE]

where $\boldsymbol{b}(t,\boldsymbol{x})$ is the body force, $\rho$ is the density and $\partial_{S}W^{\epsilon}$ is the derivative of potential with respect to the strain.

We prescribe the zero Dirichlet condition on the boundary of $D$

[TABLE]

where we have denoted the boundary by $\partial D$. We extend the zero boundary condition outside $D$ to $\mathbb{R}^{3}$.

The peridynamic equation, boundary conditions, and initial conditions

[TABLE]

determine the peridynamic evolution $\boldsymbol{u}(t,\boldsymbol{x})$.

Peridynamics energy

The total energy $\mathcal{E}^{\epsilon}(\boldsymbol{u})(t)$ is given by the sum of kinetic and potential energy given by

[TABLE]

where potential energy $PD^{\epsilon}$ is given by

[TABLE]

Differentiation of Equation 4 gives the identity

[TABLE]

where $(\cdot,\cdot)$ is the inner product on $L^{2}(\mathbb{R}^{d},D)$ and $\|\cdot\|_{L^{2}(\mathbb{R}^{d},D)}$ is the associated norm.

2.1 Nonlocal potential

We consider the nonlocal two point interaction potential density $W^{\epsilon}$ of the form

[TABLE]

where $f:\mathbb{R}^{+}\to\mathbb{R}$ is assumed to be positive, smooth and concave with following properties

[TABLE]

The potential $W^{\epsilon}(S,\boldsymbol{y}-\boldsymbol{x})$ is of double well type and convex near the origin where it has one well and concave and bounded at infinity where it has the second well. $J^{\epsilon}(\left\lvert\boldsymbol{y}-\boldsymbol{x}\right\rvert)$ models the influence of separation between points $\boldsymbol{y}$ and $\boldsymbol{x}$. We define $J^{\epsilon}$ by rescaling $J(\left\lvert\boldsymbol{\xi}\right\rvert)$, i.e. $J^{\epsilon}(\left\lvert\boldsymbol{\xi}\right\rvert)=J(\left\lvert\boldsymbol{\xi}\right\rvert/\epsilon)$. Here $J$ is zero outside the ball $H_{1}(\mathbf{0})$ and satisfies $0\leq J(\left\lvert\boldsymbol{\xi}\right\rvert)\leq M$ for all $\boldsymbol{\xi}\in H_{1}(\mathbf{0})$. The domain function $\omega$ enforces boundary conditions on $\partial_{S}W^{\epsilon}$ at the boundary of the body $D$. Here the boundary is denoted by $\partial D$ and $\omega$ is a nonnegative differentiable function $0\leq\omega\leq 1$. On the boundary $\omega=0$ and $\omega=1$ for points $\boldsymbol{x}$ inside $D$ with distance greater than $\epsilon$ away from the boundary. We continue $\omega$ by zero for all points outside $D$.

The potential described in Equation 6 gives the convex-concave dependence of $W(S,\boldsymbol{y}-\boldsymbol{x})$ on the strain $S$ for fixed $\boldsymbol{y}-\boldsymbol{x}$, see Figure 1. Initially the force is elastic for small strains and then softens as the strain becomes larger. The critical strain where the force between $\boldsymbol{x}$ and $\boldsymbol{y}$ begins to soften is given by $S_{c}(\boldsymbol{y},\boldsymbol{x}):=\bar{r}/\sqrt{\left\lvert\boldsymbol{y}-\boldsymbol{x}\right\rvert}$ and the force decreases monotonically for

[TABLE]

Here $\bar{r}$ is the inflection point of $r:\to f(r^{2})$ and is the root of following equation

[TABLE]

2.2 Existence of solution

Let $C^{0,\gamma}(D;\mathbb{R}^{d})$ be the Hölder space with exponent $\gamma\in(0,1]$. The closure of continuous functions with compact support on $D$ in the supremum norm is denoted by $C_{0}(D)$. We identify functions in $C_{0}(D)$ with their unique continuous extensions to $\overline{D}$. It is easily seen that functions belonging to this space take the value zero on the boundary of $D$, see e.g. [10]. We introduce $C_{0}^{0,\gamma}(D)=C^{0,\gamma}(D)\cap C_{0}(D)$. In this paper we extend all functions in $C_{0}^{0,\gamma}(D)$ by zero outside $D$. The norm of $\boldsymbol{u}\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$ is taken to be

[TABLE]

where $\left[\boldsymbol{u}\right]_{C^{0,\gamma}(D;\mathbb{R}^{d})}$ is the Hölder semi norm and given by

[TABLE]

and $C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$ is a Banach space with this norm. Here we make the hypothesis that the domain function $\omega$ belongs to $C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$.

We write the evolution Equation 1 as an equivalent first order system with $y_{1}(t)=\boldsymbol{u}(t)$ and $y_{2}(t)=\boldsymbol{v}(t)$ with $\boldsymbol{v}(t)=\partial_{t}\boldsymbol{u}(t)$. Let $y=(y_{1},y_{2})^{T}$ where $y_{1},y_{2}\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$ and let $F^{\epsilon}(y,t)=(F^{\epsilon}_{1}(y,t),F^{\epsilon}_{2}(y,t))^{T}$ such that

[TABLE]

The initial boundary value associated with the evolution Equation 1 is equivalent to the initial boundary value problem for the first order system given by

[TABLE]

with initial condition given by $y(0)=(\boldsymbol{u}_{0},\boldsymbol{v}_{0})^{T}\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})\times C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$.

The function $F^{\epsilon}(y,t)$ satisfies the Lipschitz continuity given by the following theorem.

Proposition 2.1**.**

**Lipschitz continuity and bound

Let $X=C^{0,\gamma}_{0}(D;\mathbb{R}^{d})\times C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$. The function $F^{\epsilon}(y,t)=(F^{\epsilon}_{1},F^{\epsilon}_{2})^{T}$, as defined in Equation 8 and Equation 9, is Lipschitz continuous in any bounded subset of $X$. We have, for any $y,z\in X$ and $t>0$,**

[TABLE]

where $L_{1},L_{2}$ are independent of $\boldsymbol{u},\boldsymbol{v}$ and depend on peridynamic potential function $f$ and influence function $J$ and the exponent $\alpha(\gamma)$ is given by

[TABLE]

Furthermore for any $y\in X$ and any $t\in[0,T]$, we have the bound

[TABLE]

*where $b=\sup_{t}{\left\lVert\boldsymbol{b}(t)\right\rVert}_{C^{0,\gamma}(D;\mathbb{R}^{d})}$ and $L_{3}$ is independent of $y$. *

We easily see that on choosing $z=0$ in Proposition 2.1 that $-\nabla PD^{\epsilon}(\boldsymbol{u})(\boldsymbol{x})$ is in $C^{0,\gamma}(D;\mathbb{R}^{3})$ provided that $\boldsymbol{u}$ belongs to $C^{0,\gamma}(D;\mathbb{R}^{3})$. Since $-\nabla PD^{\epsilon}(\boldsymbol{u})(\boldsymbol{x})$ takes the value [math] on $\partial D$ we conclude that $-\nabla PD^{\epsilon}(\boldsymbol{u})(\boldsymbol{x})$ belongs to $C^{0,\gamma}_{0}(D;\mathbb{R}^{3})$.

In Theorem 6.1 of [23], the Lipschitz property of a peridynamic force is shown in $X=L^{2}(D;\mathbb{R}^{d})\times L^{2}(D;\mathbb{R}^{d})$. It is given by

[TABLE]

for all $y,z\in L^{2}_{0}(D;\mathbb{R}^{d})^{2}$. For this case $L$ does not depend on $\boldsymbol{u},\boldsymbol{v}$. We now state the existence theorem.

The following theorem gives the existence and uniqueness of solution in any given time domain $I_{0}=(-T,T)$.

Theorem 2.2**.**

**Existence and uniqueness of Hölder solutions of cohesive dynamics over finite time intervals

For any initial condition $x_{0}\in X={C^{0,\gamma}_{0}(D;\mathbb{R}^{d})\times C^{0,\gamma}_{0}(D;\mathbb{R}^{d})}$, time interval $I_{0}=(-T,T)$, and right hand side $\boldsymbol{b}(t)$ continuous in time for $t\in I_{0}$ such that $\boldsymbol{b}(t)$ satisfies $\sup_{t\in I_{0}}{||\boldsymbol{b}(t)||_{C^{0,\gamma}}}<\infty$, there is a unique solution $y(t)\in C^{1}(I_{0};X)$ of**

[TABLE]

or equivalently

[TABLE]

*where $y(t)$ and $y^{\prime}(t)$ are Lipschitz continuous in time for $t\in I_{0}$. *

The proof of this theorem is given in section 4. We now describe the finite difference scheme and analyze its convergence to Hölder continuous solutions of cohesive dynamics.

3 Finite difference approximation

In this section, we present the finite difference scheme and compute the rate of convergence. We first consider the semi-discrete approximation and prove the bound on energy of semi-discrete evolution in terms of initial energy and the work done by body forces.

Let $h$ be the size of a mesh and $\Delta t$ be the size of time step. We will keep $\epsilon$ fixed and assume that $h<\epsilon<1$. Let $D_{h}=D\cap(h\mathbb{Z})^{d}$ be the discretization of material domain. Let $i\in\mathbb{Z}^{d}$ be the index such that $\boldsymbol{x}_{i}=hi\in D$. Let $U_{i}$ is the unit cell of volume $h^{d}$ corresponding to the grid point $\boldsymbol{x}_{i}$. The exact solution evaluated at grid points is denoted by $(\boldsymbol{u}_{i}(t),\boldsymbol{v}_{i}(t))$.

3.1 Time discretization

Let $[0,T]\cap(\Delta t\mathbb{Z})$ be the discretization of time domain where $\Delta t$ is the size of time step. Denote fully discrete solution at $(t^{k}=k\Delta t,\boldsymbol{x}_{i}=ih)$ as $(\hat{\boldsymbol{u}}^{k}_{i},\hat{\boldsymbol{v}}^{k}_{i})$. Similarly, the exact solution evaluated at grid points is denoted by $(\boldsymbol{u}^{k}_{i},\boldsymbol{v}^{k}_{i})$. We enforce boundary condition $\hat{\boldsymbol{u}}^{k}_{i}=\mathbf{0}$ for all $\boldsymbol{x}_{i}\notin D$ and for all $k$.

We begin with the forward Euler time discretization, with respect to velocity, and the finite difference scheme for $(\hat{\boldsymbol{u}}^{k}_{i},\hat{\boldsymbol{v}}^{k}_{i})$ is written

[TABLE]

The scheme is complemented with the discretized initial conditions $\hat{\boldsymbol{u}}^{0}_{i}=(\hat{\boldsymbol{u}}_{0})_{i}$ and $\hat{\boldsymbol{v}}^{0}_{i}=(\hat{\boldsymbol{v}}_{0})_{i}$. If we substitute Equation 14 into Equation 15, we get standard Central difference scheme in time for second order in time differential equation. Here we have assumed, without loss of generality, $\rho=1$.

The piecewise constant extensions of the discrete sets $\{\hat{\boldsymbol{u}}^{k}_{i}\}_{i\in\mathbb{Z}^{d}}$ and $\{\hat{\boldsymbol{v}}^{k}_{i}\}_{i\in\mathbb{Z}^{d}}$ are given by

[TABLE]

In this way we represent the finite difference solution as a piecewise constant function. We will show this function provides an $L^{2}$ approximation of the exact solution.

3.1.1 Convergence results

In this section we provide upper bounds on the rate of convergence of the discrete approximation to the solution of the peridynamic evolution. The $L^{2}$ approximation error $E^{k}$ at time $t^{k}$, for $0<t^{k}\leq T$ is defined as

[TABLE]

The upper bound on the convergence rate of the approximation error is given by the following theorem.

Theorem 3.1**.**

**Convergence of finite difference approximation (forward Euler time discretization)

Let $\epsilon>0$ be fixed. Let $(\boldsymbol{u},\boldsymbol{v})$ be the solution of peridynamic equation Equation 10. We assume $\boldsymbol{u},\boldsymbol{v}\in C^{2}([0,T];C^{0,\gamma}_{0}(D;\mathbb{R}^{d}))$. Then the finite difference scheme given by Equation 14 and Equation 15 is consistent in both time and spatial discretization and converges to the exact solution uniformly in time with respect to the $L^{2}(D;\mathbb{R}^{d})$ norm. If we assume the error at the initial step is zero then the error $E^{k}$ at time $t^{k}$ is bounded and to leading order in the time step $\Delta t$ satisfies**

[TABLE]

*where constant $C_{s}$ and $C_{t}$ are independent of $h$ and $\Delta t$ and $C_{s}$ depends on the Hölder norm of the solution and $C_{t}$ depends on the $L^{2}$ norms of time derivatives of the solution. *

Here we have assumed the initial error to be zero for ease of exposition only.

We remark that the explicit constants leading to Equation 16 can be large. The inequality that delivers Equation 16 is given to leading order by

[TABLE]

where the constants $\bar{C}$, $C_{t}$ and $C_{s}$ are given by Equation 35, Equation 37, and Equation 38. The explicit constant $C_{t}$ depends on the spatial $L^{2}$ norm of the time derivatives of the solution and $C_{s}$ depends on the spatial Hölder continuity of the solution and the constant $\bar{C}$. This constant is bounded independently of horizon $\epsilon$. Although the constants are necessarily pessimistic they deliver a-priori error estimates and an example is discussed in section 6.

An identical convergence rate can be established for the general one step scheme and we state it below.

Theorem 3.2**.**

**Convergence of finite difference approximation (General single step time discretization)

Let us assume that the hypothesis of Theorem 3.1 holds. Fix $\theta\in[0,1]$, and let $(\hat{\boldsymbol{u}}^{k},\hat{\boldsymbol{v}}^{k})^{T}$ be the solution of following finite difference equation**

[TABLE]

Then, for any fixed $\theta\in[0,1]$, there exists a constant $K>0$ independent of $(\hat{\boldsymbol{u}}^{k},\hat{\boldsymbol{v}}^{k})^{T}$ and $(\boldsymbol{u}^{k},\boldsymbol{v}^{k})^{T}$, such that for $\Delta t<K{\epsilon^{2}}$ the finite difference scheme given by Equation 18 and Equation 19 is consistent in time and spatial discretization. if we assume the error at the initial step is zero then the error $E^{k}$ at time $t^{k}$ is bounded and satisfies

[TABLE]

The constant $K$ is given by the explicit formula $K=1/\bar{C}$ where $\bar{C}$ is described by equation Equation 35. Furthermore for the Crank Nicholson scheme, $\theta=1/2$, if we assume the solutions $\boldsymbol{u},\boldsymbol{v}$ belong to $C^{3}([0,T];C^{0,\gamma}_{0}(D;\mathbb{R}^{d}))$, then the approximation error $E^{k}$ satisfies

[TABLE]

*where $\bar{C}_{t}$ is independent of $\Delta t$ and $h$ and is given by Equation 49. *

As before we assume that the error in the initial data is zero for ease of exposition. The proofs of Theorem 3.1 and Theorem 3.2 are given in the following sections.

Remark. In Theorem 3.2, we have stated a condition on $\Delta t$ for which the convergence estimate holds. This condition naturally occurs in the analysis and is related to the Lipschitz continuity of the peridnamic force with respect to the $L^{2}$ norm, see Equation 41.

3.1.2 Error analysis

Theorem 3.1 and Theorem 3.2 are proved along similar lines. In both cases we define the $L^{2}$-projections of the actual solutions onto the space of piecewise constant functions defined over the cells $U_{i}$. These are given as follows. Let $(\tilde{\boldsymbol{u}}^{k}_{i},\tilde{\boldsymbol{v}}^{k}_{i})$ be the average of the exact solution $(\boldsymbol{u}^{k},\boldsymbol{v}^{k})$ in the unit cell $U_{i}$ given by

[TABLE]

and the $L^{2}$ projection of the solution onto piecewise constant functions are $(\tilde{\boldsymbol{u}}^{k},\tilde{\boldsymbol{v}}^{k})$ given by

[TABLE]

The error between $(\hat{\boldsymbol{u}}^{k},\hat{\boldsymbol{v}}^{k})^{T}$ with $(\boldsymbol{u}(t^{k}),\boldsymbol{v}(t^{k}))^{T}$ is now split into two parts. From the triangle inequality, we have

[TABLE]

In subsubsection 3.1.3 and subsection 3.2 we will show that the error between the $L^{2}$ projections of the actual solution and the discrete approximation for both forward Euler and implicit one step methods decay according to

[TABLE]

In what follows we can estimate the terms

[TABLE]

and show they go to zero at a rate of $h^{\gamma}$ uniformly in time. The estimates given by Equation 22 together with the $O(h^{\gamma})$ estimates for Equation 23 establish Theorem 3.1 and Theorem 3.2. We now establish the $L^{2}$ estimates for the differences $\tilde{\boldsymbol{u}}^{k}-\boldsymbol{u}(t^{k})$ and $\tilde{\boldsymbol{v}}^{k}-\boldsymbol{v}(t^{k})$.

We write

[TABLE]

where we used Cauchy’s inequality and Jensen’s inequality. For $\boldsymbol{x},\boldsymbol{y}\in U_{i}$, $\left\lvert\boldsymbol{x}-\boldsymbol{y}\right\rvert\leq ch$, where $c=\sqrt{2}$ for $d=2$ and $c=\sqrt{3}$ for $d=3$. Since $\boldsymbol{u}\in C^{0,\gamma}_{0}$ we have

[TABLE]

and substitution in Equation 24 gives

[TABLE]

A similar estimate can be derived for $||\tilde{\boldsymbol{v}}^{k}-\boldsymbol{v}^{k}||_{L^{2}}$ and substitution of the estimates into Equation 23 gives

[TABLE]

In the next section we establish the error estimate Equation 22 for both forward Euler and general one step schemes in subsubsection 3.1.3 and subsection 3.2.

3.1.3 Error analysis for approximation of $L^{2}$ projection of the the exact solution

In this sub-section, we estimate the difference between approximate solution $(\hat{\boldsymbol{u}}^{k},\hat{\boldsymbol{v}}^{k})$ and the $L^{2}$ projection of the exact solution onto piece wise constant functions given by $(\tilde{\boldsymbol{u}}^{k},\tilde{\boldsymbol{v}}^{k})$, see Equation 20 and Equation 21 . Let the differences be denoted by $\boldsymbol{e}^{k}(u):=\hat{\boldsymbol{u}}^{k}-\tilde{\boldsymbol{u}}^{k}$ and $\boldsymbol{e}^{k}(v):=\hat{\boldsymbol{v}}^{k}-\tilde{\boldsymbol{v}}^{k}$ and their evaluation at grid points are $\boldsymbol{e}^{k}_{i}(u):=\hat{\boldsymbol{u}}^{k}_{i}-\tilde{\boldsymbol{u}}^{k}_{i}$ and $\boldsymbol{e}^{k}_{i}(v):=\hat{\boldsymbol{v}}^{k}_{i}-\tilde{\boldsymbol{v}}^{k}_{i}$. Subtracting $(\tilde{\boldsymbol{u}}^{k+1}_{i}-\tilde{\boldsymbol{u}}^{k}_{i})/\Delta t$ from Equation 14 gives

[TABLE]

Taking the average over unit cell $U_{i}$ of the exact peridynamic equation Equation 10 at time $t^{k}$, we will get $\tilde{\boldsymbol{v}}^{k+1}_{i}-\dfrac{\partial\tilde{\boldsymbol{u}}^{k+1}_{i}}{\partial t}=0$. Therefore, the equation for $\boldsymbol{e}^{k}_{i}(u)$ is given by

[TABLE]

where we identify the discretization error as

[TABLE]

Similarly, we subtract $(\tilde{\boldsymbol{v}}^{k+1}_{i}-\tilde{\boldsymbol{v}}^{k}_{i})/\Delta t$ from Equation 15 and add and subtract terms to get

[TABLE]

where we identify $\tau^{k}_{i}(v)$ as follows

[TABLE]

Note that in $\tau^{k}(u)$ we have $\dfrac{\partial\tilde{\boldsymbol{u}}^{k+1}_{i}}{\partial t}$ and from the exact peridynamic equation, we have

[TABLE]

Combining subsubsection 3.1.3, Equation 29, and Equation 30, to get

[TABLE]

The spatial discretization error $\sigma^{k}_{i}(u)$ and $\sigma^{k}_{i}(v)$ is given by

[TABLE]

We finally have

[TABLE]

We now show the consistency and stability properties of the numerical scheme.

3.1.4 Consistency

We deal with the error in time discretization and the error in spatial discretization error separately. The time discretization error follows easily using the Taylor’s series while spatial the discretization error uses properties of the nonlinear peridynamic force.

**Time discretization: **We first estimate the time discretization error. A Taylor series expansion is used to estimate $\tau^{k}_{i}(u)$ as follows

[TABLE]

Computing the $L^{2}$ norm of $\tau^{k}_{i}(u)$ and using Jensen’s inequality gives

[TABLE]

Similarly, we have

[TABLE]

**Spatial discretization: **We now estimate the spatial discretization error. Substituting the definition of $\tilde{\boldsymbol{v}}^{k}$ and following the similar steps employed in subsubsection 3.1.2, gives

[TABLE]

Taking the $L^{2}$ norm of error $\sigma^{k}_{i}(v)$ and substituting the estimate above delivers

[TABLE]

Now we estimate $\left\lvert\sigma^{k}_{i}(u)\right\rvert$. We use the notation $\bar{\boldsymbol{u}}^{k}(\boldsymbol{x}):=\boldsymbol{u}^{k}(\boldsymbol{x}+\epsilon\boldsymbol{\xi})-\boldsymbol{u}^{k}(\boldsymbol{x})$ and $\overline{\tilde{\boldsymbol{u}}}^{k}(\boldsymbol{x}):=\tilde{\boldsymbol{u}}(\boldsymbol{x}+\epsilon\boldsymbol{\xi})-\tilde{\boldsymbol{u}}^{k}(\boldsymbol{x})$ and choose $\boldsymbol{u}=\boldsymbol{u}^{k}$ and $\boldsymbol{v}=\tilde{\boldsymbol{u}}^{k}$ in Equation 76 to find that

[TABLE]

Here $C_{2}$ is the maximum of the second derivative of the profile describing the potential given by Equation 71. Following the earlier analysis, see subsubsection 3.1.2, we find that

[TABLE]

For reference, we define the constant

[TABLE]

We now focus on subsubsection 3.1.4. We substitute the above two inequalities to get

[TABLE]

Therefore, we have

[TABLE]

This completes the proof of consistency of numerical approximation.

3.1.5 Stability

Let $e^{k}$ be the total error at the $k^{\text{th}}$ time step. It is defined as

[TABLE]

To simplify the calculations, we define new term $\tau$ as

[TABLE]

From our consistency analysis, we know that to leading order

[TABLE]

where,

[TABLE]

We take $L^{2}$ norm of Equation 26 and subsubsection 3.1.3 and add them. Noting the definition of $\tau$ as above, we get

[TABLE]

We only need to estimate the last term in above equation. Similar to the subsubsection 3.1.4, we have

[TABLE]

By $\boldsymbol{e}^{k}(u)(\boldsymbol{x})$ we mean evaluation of piecewise extension of set $\{\boldsymbol{e}^{k}_{i}(u)\}_{i}$ at $\boldsymbol{x}$. We proceed further as follows

[TABLE]

Using inequality $|ab|\leq(\left\lvert a\right\rvert^{2}+\left\lvert b\right\rvert^{2})/2$, we get

[TABLE]

and

[TABLE]

Since $\boldsymbol{e}^{k}(u)(\boldsymbol{x})=\sum_{\boldsymbol{x}_{i}\in D}\boldsymbol{e}^{k}_{i}(u)\chi_{U_{i}}(\boldsymbol{x})$, we have

[TABLE]

where $\bar{C}$ is given by Equation 35. In summary Equation 40 shows the Lipschitz continuity of the peridynamic force with respect to the $L^{2}$ norm, see Equation 13, expressed in this context as

[TABLE]

Finally, we substitute above inequality in subsubsection 3.1.5 to get

[TABLE]

We add positive quantity $\Delta t||e^{k+1}(u)||_{L^{2}(D;\mathbb{R}^{d})}+\Delta t6\bar{C}/\epsilon^{2}||\boldsymbol{e}^{k}(v)||_{L^{2}(D;\mathbb{R}^{d})}$ to the right side of above equation, to get

[TABLE]

We recursively substitute $e^{j}$ on above as follows

[TABLE]

Since $1/(1-\Delta t)=1+\Delta t+\Delta t^{2}+O(\Delta t^{3})$, we have

[TABLE]

Now, for any $k\leq T/\Delta t$, using identity $(1+a)^{k}\leq\exp[ka]$ for $a\leq 0$, we have

[TABLE]

We write above equation in more compact form as follows

[TABLE]

We use above estimate in Equation 42 and get following inequality for $e^{k}$

[TABLE]

where we used the fact that $1/(1-\Delta t)=1+\Delta t+O(\Delta t^{2})$.

Assuming the error in initial data is zero, i.e. $e^{0}=0$, and noting the estimate of $\tau$ in Equation 36, we have

[TABLE]

and we conclude to leading order that

[TABLE]

Here the constants $C_{t}$ and $C_{s}$ are given by Equation 37 and Equation 38.This shows the stability of the numerical scheme. We now address the general one step time discretization.

3.2 Extension to the implicit schemes

Let $\theta\in[0,1]$ be the parameter which controls the contribution of the implicit and explicit scheme. Let $(\hat{\boldsymbol{u}}^{k},\hat{\boldsymbol{v}}^{k})$ be the solution of Equation 18 and Equation 19 for given fixed $\theta$.

The forward Euler scheme, backward Euler scheme, and Crank Nicholson scheme correspond to the choices $\theta=0$, $\theta=1$, and $\theta=1/2$ respectively.

To simplify the equations, we define $\Theta$ acting on discrete set $\{f^{k}\}_{k}$ as $\Theta f^{k}:=(1-\theta)f^{k}+\theta f^{k+1}$. By $\Theta\left\lVert f^{k}\right\rVert$, we mean $(1-\theta)\left\lVert f^{k}\right\rVert+\theta\left\lVert f^{k+1}\right\rVert$. Following the same steps as in the case of forward Euler, we write down the equation for $\boldsymbol{e}^{k}_{i}(u):=\hat{\boldsymbol{u}}^{k}_{i}-\tilde{\boldsymbol{u}}^{k}_{i}$ and $\boldsymbol{e}^{k}_{i}(v):=\hat{\boldsymbol{v}}^{k}_{i}-\tilde{\boldsymbol{v}}^{k}_{i}$ as follows

[TABLE]

where $\tau^{k}_{i}(v),\sigma^{k}_{i}(u),\sigma^{k}_{i}(v)$ are defined in Equation 29, Equation 31, and Equation 32 respectively. In this section $\tau^{k}(u)$ is defined as follows

[TABLE]

We take the $L^{2}$ norm of $\boldsymbol{e}^{k}(u)(\boldsymbol{x})$ and $\boldsymbol{e}^{k}(v)(\boldsymbol{x})$ and for brevity we denote the $L^{2}$ norm by $||\cdot||$. Recall that $\boldsymbol{e}^{k}(u)$ and $\boldsymbol{e}^{k}(v)$ are the piecewise constant extension of $\{\boldsymbol{e}^{k}_{i}(u)\}_{i}$ and $\{\boldsymbol{e}^{k}_{i}(v)\}_{i}$ and we get

[TABLE]

From our consistency analysis, we have

[TABLE]

where $C_{t}$ and $C_{s}$ are given by Equation 37 and Equation 38. Since $0\leq 1-\theta\leq 1$ and $0\leq\theta\leq 1$ for all $\theta\in[0,1]$, we have

[TABLE]

**Crank Nicholson scheme: **If $\theta=1/2$, and if $\boldsymbol{u},\boldsymbol{v}\in C^{3}([0,T];C^{0,\gamma}(D;\mathbb{R}^{d}))$, then we can show that

[TABLE]

A similar result holds for $1/2\tau^{k}_{i}(v)+1/2\tau^{k+1}_{i}(v)$. Therefore, the consistency error will be bounded by $\bar{C}_{t}\Delta t^{2}+C_{s}h^{\gamma}/\epsilon^{2}$ with

[TABLE]

and $C_{s}$ is given by Equation 38.

We now estimate Equation 47. Similar to Equation 40, we have

[TABLE]

where $\bar{C}$ is the constant given by Equation 35. Let $e^{k}:=\left\lVert\boldsymbol{e}^{k}(u)\right\rVert+\left\lVert\boldsymbol{e}^{k}(v)\right\rVert$. Adding Equation 46 and Equation 47 and noting Equation 50, Equation 51, and Equation 48, we get

[TABLE]

where we assumed $\bar{C}/\epsilon^{2}\geq 1$. We further simplify the equation and write

[TABLE]

where we have assumed that $1-\Delta t\theta\bar{C}/\epsilon^{2}>0$, i.e.

[TABLE]

Thus, for fixed $\epsilon>0$, the error calculation in this section applies when the time step $\Delta t$ satisfies Equation 53. We now define $a$ and $b$ by

[TABLE]

We use the fact that, for $\Delta t$ small, $(1-\alpha\Delta t)^{-1}=1+\alpha\Delta t+\alpha^{2}(\Delta t)^{2}+O((\Delta t)^{3})$, to get

[TABLE]

Now since $\Delta t<\epsilon^{2}/\bar{C}$, we have

[TABLE]

We have the estimates for $a$ given by

[TABLE]

Therefore, for any $k\leq T/\Delta t$, we have

[TABLE]

where we simplified the bound by incorporating Equation 53. Then, from Equation 52, we get

[TABLE]

From the estimates on $a^{k}$, we have

[TABLE]

Combining Equation 54 and Equation 55, to get

[TABLE]

Since $\tau=O(C_{t}\Delta t+C_{s}h^{\gamma}/\epsilon^{2})$, we conclude that, for any $\epsilon>0$ fixed,

[TABLE]

Where we assumed $e^{0}=0$. Similarly, for $\theta=1/2$, we have $\sup_{k}e^{k+1}\leq O(\bar{C}_{t}(\Delta t)^{2}+C_{s}h^{\gamma}/\epsilon^{2})$. Therefore, the scheme is stable and consistent for any $\theta\in[0,1]$.

3.3 Stability of the energy for the semi-discrete approximation

We first spatially discretize the peridynamics equation Equation 1. Let $\{\hat{\boldsymbol{u}}_{i}(t)\}_{i,\boldsymbol{x}_{i}\in D}$ denote the semi-discrete approximate solution which satisfies following, for all $t\in[0,T]$ and $i$ such that $\boldsymbol{x}_{i}\in D$,

[TABLE]

where $\hat{\boldsymbol{u}}(t)$ is the piecewise constant extension of discrete set $\{\hat{\boldsymbol{u}}_{i}(t)\}_{i}$ and is defined as

[TABLE]

The scheme is complemented with the discretized initial conditions $\hat{\boldsymbol{u}}_{i}(0)=\boldsymbol{u}_{0}(\boldsymbol{x}_{i})$ and $\hat{\boldsymbol{v}}_{i}(0)=\boldsymbol{v}_{0}(\boldsymbol{x}_{i})$. We apply boundary condition by setting $\hat{\boldsymbol{u}}_{i}(t)=\mathbf{0}$ for all $t$ and for all $\boldsymbol{x}_{i}\notin D$.

We have the stability of semi-discrete evolution.

Theorem 3.3**.**

**Energy stability of the semi-discrete approximation

Let $\{\hat{\boldsymbol{u}}_{i}(t)\}_{i}$ satisfy Equation 56 and $\hat{\boldsymbol{u}}(t)$ is its piecewise constant extension. Similarly let $\hat{\boldsymbol{b}}(t,\boldsymbol{x})$ denote the piecewise constant extension of $\{\boldsymbol{b}(t,\boldsymbol{x}_{i})\}_{i,\boldsymbol{x}_{i}\in D}$. Then the peridynamic energy $\mathcal{E}^{\epsilon}$ as defined in Equation 4 satisfies, $\forall t\in[0,T]$,**

[TABLE]

*The constant $C$, defined in Equation 63, is independent of $\epsilon$ and $h$. *

Proof 3.4**.**

We multiply Equation 56 by $\chi_{U_{i}}(\boldsymbol{x})$ and sum over $i$ and use definition of piecewise constant extension in Equation 57 to get

[TABLE]

where $-\boldsymbol{\nabla}\hat{PD^{\epsilon}}(\hat{\boldsymbol{u}}(t))(\boldsymbol{x})$ and $\hat{\boldsymbol{b}}(t,\boldsymbol{x})$ are given by

[TABLE]

We define set as follows

[TABLE]

We use following result which we will show after few steps

[TABLE]

We then have

[TABLE]

Multiplying above with $\dot{\hat{\boldsymbol{u}}}(t)$ and integrating over $D$ to get

[TABLE]

Consider energy $\mathcal{E}^{\epsilon}(\hat{\boldsymbol{u}})(t)$ given by Equation 4 and note the identity Equation 5, to have

[TABLE]

where we used Hölder inequality in last step. Since $PD^{\epsilon}(\boldsymbol{u})$ is positive for any $\boldsymbol{u}$, we have

[TABLE]

Using above, we get

[TABLE]

Let $\delta>0$ is some arbitrary but fixed real number and let $A(t)=\delta+\mathcal{E}^{\epsilon}(\hat{\boldsymbol{u}})(t)$. Then

[TABLE]

Using the fact that $\frac{1}{\sqrt{A(t)}}\frac{d}{dt}A(t)=2\frac{d}{dt}\sqrt{A(t)}$, we have

[TABLE]

where we used bound on $||\sigma(s)||_{L^{2}(D;\mathbb{R}^{d})}$ from Equation 60. Noting that $\delta>0$ is arbitrary, we send it to zero to get

[TABLE]

and Equation 58 follows by taking square of above equation.

It remains to show Equation 60. To simplify the calculations, we use following notations: let $\boldsymbol{\xi}\in H_{1}(\mathbf{0})$ and let

[TABLE]

With above notations and using expression of $-\boldsymbol{\nabla}PD^{\epsilon}$ from subsection 4.1, we have for $\boldsymbol{x}\in U_{i}$

[TABLE]

Using the fact that $0\leq\omega(\boldsymbol{x})\leq 1$ and $|F_{1}^{\prime}(r)|\leq C_{1}$, where $C_{1}$ is $\sup_{r}|F_{1}^{\prime}(r)|$, we get

[TABLE]

where $\bar{J}_{1/2}=(1/\omega_{d})\int_{H_{1}(\mathbf{0})}J(|\boldsymbol{\xi}|)|\boldsymbol{\xi}|^{-1/2}d\boldsymbol{\xi}$.

Taking the $L^{2}$ norm of $\sigma(t,\boldsymbol{x})$, we get

[TABLE]

thus

[TABLE]

where

[TABLE]

*This completes the proof. *

3.4 Local instability under radial perturbations

We observe that both explicit and implicit schemes treated in previous sections show that any increase in local truncation error is controlled at each time step. From the proofs above (and the general approximation theory for ODE), this control is adequate to establish convergence rates as $\Delta t\rightarrow 0$. We now comment on a source of error that can grow with time steps in regions where the strain is large and peridynamic bonds between material points begin to soften.

We examine the Jacobian matrix of the peridynamic system associated by perturbing about a displacement field and seek to understand the stability of the perturbation. Suppose the solution is near the displacement field $\overline{\boldsymbol{u}}(\boldsymbol{x})$ and let $\boldsymbol{s}(t,\boldsymbol{x})=\boldsymbol{u}(t,\boldsymbol{x})-\overline{\boldsymbol{u}}(\boldsymbol{x})$ be the perturbation. We write the associated strain as $S(\boldsymbol{y},\boldsymbol{x};\overline{\boldsymbol{u}})$ and $S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{s})$. Expanding the peridynamic force in Taylor series about $\overline{\boldsymbol{u}}$ assuming $\boldsymbol{s}$ is small gives

[TABLE]

where $\mathcal{W}^{\epsilon}(S,\boldsymbol{y}-\boldsymbol{x})=W^{\epsilon}(S,\boldsymbol{y}-\boldsymbol{x})/|\boldsymbol{y}-\boldsymbol{x}|$ and $W^{\epsilon}$ is given by Equation 6.

To recover a local stability formula in terms of a spectral radius we consider local radial perturbations $\boldsymbol{s}$ with spatially constant strain $S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{s})$ of the form $S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{s})=-\delta(t){\boldsymbol{\mu}}\cdot\boldsymbol{e}$ where $\boldsymbol{\mu}$ is in $\mathbb{R}^{d}$ and $\boldsymbol{s}$ has radial variation about $\boldsymbol{x}$ with $\boldsymbol{s}(\boldsymbol{y})=\delta(t){\boldsymbol{\mu}}(1-|\boldsymbol{y}-\boldsymbol{x}|)$. This delivers the local ODE

[TABLE]

where the stability matrix $A$ is selfadjoint and given by

[TABLE]

and

[TABLE]

A stability criterion for the perturbation is obtained on analyzing the linear system $\delta^{\prime\prime}(t)\boldsymbol{\mu}=A\delta(t){\boldsymbol{\mu}}$. Writing it as a 1st order system gives

[TABLE]

where $\boldsymbol{\mu}$ is a vector in $\mathbb{R}^{d}$. The eigenvalues of $A$ are real and denoted by $\lambda_{i}$, $i=1,\ldots,d$ and the associated eigenvectors are denoted by $\boldsymbol{v}^{i}$. Choosing $\boldsymbol{\mu}=\boldsymbol{v}^{i}$ gives

[TABLE]

Applying the Forward Euler method to this system gives the discrete iterative system

[TABLE]

The spectral radius of the matrix associated with this iteration is

[TABLE]

It is easy to see that the spectral radius is larger than $1$ for any choice of $\lambda_{i}$ and we conclude local instability for the forward Euler scheme under radial perturbation.

For the implicit scheme given by backward Euler we get the discrete iterative system

[TABLE]

and

[TABLE]

The spectral radius for the iteration matrix is

[TABLE]

where $\theta=1-\lambda_{i}(\Delta t)^{2}$. If we suppose that the stability matrix $A$ is not negative definite and there is a $\lambda_{j}>0$, then the spectral radius is larger than one, i.e,

[TABLE]

Thus it follows from Equation 65 that we can have local instability of the backward Euler scheme for radial perturbations. Inspection of Equation 64 shows the sign of the eigenvalues of the matrix $A$ depend explicitly on the sign of $\partial^{2}_{{S}}\mathcal{W}^{\epsilon}({S}(\boldsymbol{y},\boldsymbol{x};\overline{\boldsymbol{u}})$. It is shown in [23] that

[TABLE]

From the model we see that bonds are loosing stiffness when $|{S}(\boldsymbol{y},\boldsymbol{x};\overline{\boldsymbol{u}})|>S_{c}$ and the points for which $A$ is non negative definite correspond to points where Equation 67 holds for a preponderance of bonds inside the horizon. We conclude noting that both explicit and implicit schemes treated in previous sections have demonstrated convergence rates $O((C_{t}\Delta t+C_{s}h^{\gamma}/\epsilon^{2}))$ as $\Delta t\rightarrow 0$. However, the results of this section show that the error can grow with time for this type of radial perturbation.

4 Lipschitz continuity in Hölder norm and existence of a solution

In this section, we prove Proposition 2.1, Theorem 4.1, and Theorem 2.2.

4.1 Proof of Proposition 1

Let $I=[0,T]$ be the time domain and $X=C^{0,\gamma}_{0}(D;\mathbb{R}^{d})\times C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$. Recall that $F^{\epsilon}(y,t)=(F^{\epsilon}_{1}(y,t),F^{\epsilon}_{2}(y,t))$, where $F^{\epsilon}_{1}(y,t)=y^{2}$ and $F^{\epsilon}_{2}(y,t)=-\boldsymbol{\nabla}PD^{\epsilon}(y^{1})+\boldsymbol{b}(t)$. Given $t\in I$ and $y=(y^{1},y^{2}),z=(z^{1},z^{2})\in X$, we have

[TABLE]

Therefore, to prove the Proposition 2.1, we only need to analyze the second term in above inequality. Let $\boldsymbol{u},\boldsymbol{v}\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$, then we have

[TABLE]

Note that the force $-\boldsymbol{\nabla}PD^{\epsilon}(\boldsymbol{u})(\boldsymbol{x})$ can be written as follows

[TABLE]

where we substituted $\partial_{S}W^{\epsilon}$ using Equation 6. In second step, we introduced the change in variable $\boldsymbol{y}=\boldsymbol{x}+\epsilon\boldsymbol{\xi}$.

Let $F_{1}:\mathbb{R}\to\mathbb{R}$ be defined as $F_{1}(S)=f(S^{2})$. Then $F^{\prime}_{1}(S)=f^{\prime}(S^{2})2S$. Using the definition of $F_{1}$, we have

[TABLE]

Because $f$ is assumed to be positive, smooth, and concave, and is bounded far away, we have following bound on derivatives of $F_{1}$

[TABLE]

where $\bar{r}$ is the inflection point of $f(r^{2})$, i.e. $F^{\prime\prime}_{1}(\bar{r})=0$. $\{0,\hat{u}\}$ are the maxima of $F^{\prime\prime}_{1}(r)$. $\{\bar{u},\tilde{u}\}$ are the maxima of $F^{\prime\prime\prime}_{1}(r)$. By chain rule and by considering the assumption on $f$, we can show that $\bar{r},\hat{u},\bar{u}_{2},\tilde{u}_{2}$ exists and the $C_{1},C_{2},C_{3}$ are bounded. Figure 4, Figure 5, and Figure 6 shows the generic graphs of $F^{\prime}_{1}(r)$, $F^{\prime\prime}_{1}(r)$, and $F^{\prime\prime\prime}_{1}(r)$ respectively.

The nonlocal force $-\boldsymbol{\nabla}PD^{\epsilon}$ can be written as

[TABLE]

To simplify the calculations, we use following notation

[TABLE]

and $\overline{(\boldsymbol{u}-\boldsymbol{v})}(\boldsymbol{x})$ is defined similar to $\bar{\boldsymbol{u}}(\boldsymbol{x})$. Also, let

[TABLE]

In what follows, we will come across the integral of type $\int_{H_{1}(\mathbf{0})}J(\left\lvert\boldsymbol{\xi}\right\rvert)\left\lvert\boldsymbol{\xi}\right\rvert^{-\alpha}d\boldsymbol{\xi}$. Recall that $0\leq J(\left\lvert\boldsymbol{\xi}\right\rvert)\leq M$ for all $\boldsymbol{\xi}\in H_{1}(\mathbf{0})$ and $J(\left\lvert\boldsymbol{\xi}\right\rvert)=0$ for $\boldsymbol{\xi}\notin H_{1}(\mathbf{0})$. Therefore, let

[TABLE]

With notations above, we note that $S(\boldsymbol{x}+\epsilon\boldsymbol{\xi},\boldsymbol{x};\boldsymbol{u})=\bar{\boldsymbol{u}}(\boldsymbol{x})\cdot\boldsymbol{e}/s$. $-\boldsymbol{\nabla}PD^{\epsilon}$ can be written as

[TABLE]

We first estimate the term $\left\lvert-\boldsymbol{\nabla}PD^{\epsilon}(\boldsymbol{u})(\boldsymbol{x})-(-\boldsymbol{\nabla}PD^{\epsilon}(\boldsymbol{v})(\boldsymbol{x}))\right\rvert$ in subsection 4.1.

[TABLE]

Here we have used the fact that $|\omega(\boldsymbol{x})|\leq 1$ and for a vector $\boldsymbol{e}$ such that $\left\lvert\boldsymbol{e}\right\rvert=1$, $\left\lvert\boldsymbol{a}\cdot\boldsymbol{e}\right\rvert\leq\left\lvert\boldsymbol{a}\right\rvert$ holds and $\left\lvert\alpha\boldsymbol{e}\right\rvert\leq\left\lvert\alpha\right\rvert$ holds for all $\boldsymbol{a}\in\mathbb{R}^{d},\alpha\in\mathbb{R}$. Using the fact that $\boldsymbol{u},\boldsymbol{v}\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$, we have

[TABLE]

Substituting the estimate given above, we get

[TABLE]

where $C_{2}$ is given by Equation 71 and $\bar{J}_{1-\gamma}$ is given by Equation 74.

We now estimate the second term in subsection 4.1. To simplify notation we write $\tilde{\omega}(\boldsymbol{x},\boldsymbol{\xi})=\omega(\boldsymbol{x})\omega(\boldsymbol{x}+\epsilon\boldsymbol{\xi})$ and with the help of Equation 75, we get

[TABLE]

We analyze the integrand in above equation. We let $H$ be defined by

[TABLE]

Let $\boldsymbol{r}:[0,1]\times D\to\mathbb{R}^{d}$ be defined as

[TABLE]

Note $\partial\boldsymbol{r}(l,\boldsymbol{x})/\partial l=\bar{\boldsymbol{u}}(\boldsymbol{x})-\bar{\boldsymbol{v}}(\boldsymbol{x})$. Using $\boldsymbol{r}(l,\boldsymbol{x})$, we have

[TABLE]

Similarly, we have.

[TABLE]

Note that

[TABLE]

Combining Equation 80, Equation 81, and Equation 82, gives

[TABLE]

Adding and subtracting $\tilde{\omega}(\boldsymbol{x},\boldsymbol{\xi})F^{\prime\prime}_{1}(\boldsymbol{r}(l,\boldsymbol{x})\cdot\boldsymbol{e}/\sqrt{s})(\bar{\boldsymbol{u}}(\boldsymbol{y})-\bar{\boldsymbol{v}}(\boldsymbol{y}))$, and noting $0\leq\tilde{\omega}(\boldsymbol{x},\boldsymbol{\xi})\leq 1$ gives

[TABLE]

Estimating $H_{1}$ first. Note that $|F^{\prime\prime}_{1}(r)|\leq C_{2}$. Since $\boldsymbol{u},\boldsymbol{v}\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$, it is easily seen that

[TABLE]

Therefore, we have

[TABLE]

We now estimate $H_{2}$. We add and subtract $\tilde{\omega}(\boldsymbol{x},\boldsymbol{\xi})F^{\prime\prime}_{1}(\boldsymbol{r}(l,\boldsymbol{y})\cdot\boldsymbol{e}/\sqrt{s}))$ in $H_{2}$ to get

[TABLE]

where

[TABLE]

and

[TABLE]

Now we estimate $H_{3}$. Since $|F^{\prime\prime\prime}_{1}(r)|\leq C_{3}$, see Equation 72, we have

[TABLE]

Where we have used the fact that $\left\lvert 1-l\right\rvert\leq 1,\left\lvert l\right\rvert\leq 1$, as $l\in[0,1]$. Also, note that

[TABLE]

We combine above estimates with Equation 84, to get

[TABLE]

Next we estimate $H_{4}$. Here we add and subtract $\omega(\boldsymbol{y})\omega(\boldsymbol{x}+\epsilon\boldsymbol{\xi})$ to get

[TABLE]

Recalling that $\omega$ belongs to $C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$ and in view of the previous estimates a straight forward calculation gives

[TABLE]

Combining Equation 83, Equation 85, and Equation 86 gives

[TABLE]

Substituting $H$ in LABEL:eq:H_def, gives

[TABLE]

We combine subsection 4.1, Equation 77, and Equation 87, and get

[TABLE]

where we introduce new constants $\bar{C}_{1},\bar{C}_{2},\bar{C}_{3}$. We let $\alpha(\gamma)=0$, if $\gamma\geq 1/2$, and $\alpha(\gamma)=1/2-\gamma$, if $\gamma\leq 1/2$. One can easily verify that, for all $\gamma\in(0,1]$ and $0<\epsilon\leq 1$,

[TABLE]

To complete the proof, we combine Equation 88 and subsection 4.1, and get

[TABLE]

This proves the Lipschitz continuity of $F^{\epsilon}(y,t)$ on any bounded subset of $X$. The bound on $F^{\epsilon}(y,t)$, see Equation 12, follows easily from Equation 75. This completes the proof of Proposition 2.1.

4.2 Existence of solution in Hölder space

In this section, we prove Theorem 2.2. We begin by proving a local existence theorem. We then show that the local solution can be continued uniquely in time to recover Theorem 2.2.

The existence and uniqueness of local solutions is stated in the following theorem.

Theorem 4.1**.**

**Local existence and uniqueness

Given $X=C^{0,\gamma}_{0}(D;\mathbb{R}^{d})\times C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$, $\boldsymbol{b}(t)\in C^{0,\gamma}_{0}(D;\mathbb{R}^{d})$, and initial data $x_{0}=(\boldsymbol{u}_{0},\boldsymbol{v}_{0})\in X$. We suppose that $\boldsymbol{b}(t)$ is continuous in time over some time interval $I_{0}=(-T,T)$ and satisfies $\sup_{t\in I_{0}}{\left\lVert\boldsymbol{b}(t)\right\rVert}_{C^{0,\gamma}(D;\mathbb{R}^{d})}<\infty$. Then, there exists a time interval $I^{\prime}=(-T^{\prime},T^{\prime})\subset I_{0}$ and unique solution $y=(y^{1},y^{2})$ such that $y\in C^{1}(I^{\prime};X)$ and**

[TABLE]

or equivalently

[TABLE]

where $y(t)$ and $y^{\prime}(t)$ are Lipschitz continuous in time for $t\in I^{\prime}\subset I_{0}$.

To prove Theorem 4.1, we proceed as follows. We write $y(t)=(y^{1}(t),y^{2}(t))$ and $||y||_{X}=||y^{1}(t)||_{C^{0,\gamma}}+||y^{2}(t)||_{C^{0,\gamma}}$. Define the ball $B(0,R)=\{y\in X:\,||y||_{X}<R\}$ and choose $R>||x_{0}||_{X}$. Let $r=R-{\left\lVert x_{0}\right\rVert}_{X}$ and we consider the ball $B(x_{0},r)$ defined by

[TABLE]

see figure Figure 7.

To recover the existence and uniqueness we introduce the transformation

[TABLE]

Introduce $0<T^{\prime}<T$ and the associated set $Y(T^{\prime})$ of Hölder continuous functions taking values in $B(x_{0},r)$ for $I^{\prime}=(-T^{\prime},T^{\prime})\subset I_{0}=(-T,T)$. The goal is to find appropriate interval $I^{\prime}=(-T^{\prime},T^{\prime})$ for which $S_{x_{0}}$ maps into the corresponding set $Y(T^{\prime})$. Writing out the transformation with $y(t)\in Y(T^{\prime})$ gives

[TABLE]

and there is a positive constant $K=C/\epsilon^{2+\alpha(\gamma)}$, see Equation 12, independent of $y^{1}(t)$, for $-T^{\prime}<t<T^{\prime}$, such that estimation in Equation 92 gives

[TABLE]

and from Equation 91

[TABLE]

We write $b=\sup_{t\in I_{0}}||\boldsymbol{b}(t)||_{C^{0,\gamma}}$ and adding Equation 93 and Equation 94 gives the upper bound

[TABLE]

Since $B(x_{0},r)\subset B(0,R)$, see Equation 90, we make the choice $T^{\prime}$ so that

[TABLE]

For this choice we see that

[TABLE]

Now it is easily seen that $\theta(R)$ is increasing with $R>0$ and

[TABLE]

So given $R$ and $||x_{0}||_{X}$ we choose $T^{\prime}$ according to

[TABLE]

and set $I^{\prime}=(-T^{\prime},T^{\prime})$. We have found the appropriate time domain $I^{\prime}$ such that the transformation $S_{x_{0}}(y)(t)$ as defined in subsection 4.2 maps $Y(T^{\prime})$ into itself. We now proceed using standard arguments, see e.g. [[10], Theorem 6.10], to complete the proof of existence and uniqueness of solution for given initial data $x_{0}$ over the interval $I^{\prime}=(-T^{\prime},T^{\prime})$.

We now prove Theorem 2.2. From the proof of Theorem 4.1 above, we see that a unique local solution exists over a time domain $(-T^{\prime},T^{\prime})$ with $\frac{\theta(R)}{2}<T^{\prime}$. Since $\theta(R)\nearrow 1/K$ as $R\nearrow\infty$ we can fix a tolerance $\eta>0$ so that $[(1/2K)-\eta]>0$. Then given any initial condition with bounded Hölder norm and $b=\sup_{t\in[-T,T)}||\boldsymbol{b}(t)||_{C^{0,\gamma}}$ we can choose $R$ sufficiently large so that $||x_{0}||_{X}<R$ and $0<(1/2K))-\eta<T^{\prime}$. Thus we can always find local solutions for time intervals $(-T^{\prime},T^{\prime})$ for $T^{\prime}$ larger than $[(1/2K)-\eta]>0$. Therefore we apply the local existence and uniqueness result to uniquely continue local solutions up to an arbitrary time interval $(-T,T)$.

5 Limit behavior of Hölder solutions in the limit of vanishing nonlocality

In this section, we consider the behavior of bounded Hölder continuous solutions as the peridynamic horizon tends to zero. We find that the solutions converge to a limiting sharp fracture evolution with bounded Griffiths fracture energy and satisfy the linear elastic wave equation away from the fracture set. We look at a subset of Hölder solutions that are differentiable in the spatial variables to show that sharp fracture evolutions can be approached by spatially smooth evolutions in the limit of vanishing non locality. As $\epsilon$ approaches zero derivatives can become large but must localize to surfaces across which the limiting evolution jumps.

We consider a sequence of peridynamic horizons $\epsilon_{k}=1/k$, $k=1,\ldots$ and the associated Hölder continuous solutions $\boldsymbol{u}^{\epsilon_{k}}(t,\boldsymbol{x})$ of the peridynamic initial value problem Equation 1, Equation 2, and Equation 3. We assume that the initial conditions $\boldsymbol{u}_{0}^{\epsilon_{k}},\boldsymbol{v}_{0}^{\epsilon_{k}}$ have uniformly bounded peridynamic energy and mean square initial velocity given by

[TABLE]

Moreover we suppose that $\boldsymbol{u}_{0}^{\epsilon_{k}},\boldsymbol{v}_{0}^{\epsilon_{k}}$ are differentiable on $D$ and that they converge in $L^{2}(D;\mathbb{R})$ to $\boldsymbol{u}_{0}^{0},\boldsymbol{v}_{0}^{0}$ with bounded Griffith free energy given by

[TABLE]

where $J_{\boldsymbol{u}_{0}^{0}}$ denotes an initial fracture surface given by the jumps in the initial deformation $\boldsymbol{u}_{0}^{0}$ and $\mathcal{H}^{2}(J_{\boldsymbol{u}^{0}(t)})$ is its $2$ dimensional Hausdorff measure of the jump set. Here $\mathcal{E}\boldsymbol{u}^{0}_{0}$ is the elastic strain and ${\rm div}\,\boldsymbol{u}^{0}_{0}=Tr(\mathcal{E}\boldsymbol{u}^{0}_{0})$. The constants $\mu$, $\lambda$ are given by the explicit formulas

[TABLE]

and

[TABLE]

where $f^{\prime}(0)$ and $f_{\infty}$ are defined by Equation 7. Here $\mu=\lambda$ and is a consequence of the central force model used in cohesive dynamics. Last we suppose as in [23] that the solutions are uniformly bounded, i.e.,

[TABLE]

The Hölder solutions $\boldsymbol{u}^{\epsilon_{k}}(t,\boldsymbol{x})$ naturally belong to $L^{2}(D;\mathbb{R}^{d})$ for all $t\in[0,T]$ and we can directly apply the Grönwall inequality (equation (6.9) of [23]) together with Theorems 6.2 and 6.4 of [23] to conclude similar to Theorems 5.1 and 5.2 of [23] that there is at least one “cluster point” $\boldsymbol{u}^{0}(t,\boldsymbol{x})$ belonging to $C([0,T];L^{2}(D;\mathbb{R}^{d}))$ and subsequence, also denoted by $\boldsymbol{u}^{\epsilon_{k}}(t,\boldsymbol{x})$ for which

[TABLE]

Moreover it follows from [23] that the limit evolution $\boldsymbol{u}^{0}(t,\boldsymbol{x})$ has a weak derivative $\boldsymbol{u}_{t}^{0}(t,\boldsymbol{x})$ belonging to $L^{2}([0,T]\times D;\mathbb{R}^{d})$. For each time $t\in[0,T]$ we can apply methods outlined in [23] to find that the cluster point $\boldsymbol{u}^{0}(t,\boldsymbol{x})$ is a special function of bounded deformation (see, [4], [5]) and has bounded linear elastic fracture energy given by

[TABLE]

for $0\leq t\leq T$ where $J_{\boldsymbol{u}^{0}(t)}$ denotes the evolving fracture surface The deformation - crack set pair $(\boldsymbol{u}^{0}(t),J_{\boldsymbol{u}^{0}(t)})$ records the brittle fracture evolution of the limit dynamics.

Arguments identical to [23] show that away from sets where $|S(\boldsymbol{y},\boldsymbol{x};\boldsymbol{u}^{\epsilon_{k}})|>S_{c}$ the limit $\boldsymbol{u}^{0}$ satisfies the linear elastic wave equation. This is stated as follows: Fix $\delta>0$ and for $\epsilon_{k}<\delta$ and $0\leq t\leq T$ consider the open set $D^{\prime}\subset D$ for which points $\boldsymbol{x}$ in $D^{\prime}$ and $\boldsymbol{y}$ for which $|\boldsymbol{y}-\boldsymbol{x}|<\epsilon_{k}$ satisfy,

[TABLE]

Then the limit evolution $\boldsymbol{u}^{0}(t,\boldsymbol{x})$ evolves elastodynamically on $D^{\prime}$ and is governed by the balance of linear momentum expressed by the Navier Lamé equations on the domain $[0,T]\times D^{\prime}$ given by

[TABLE]

where the stress tensor $\boldsymbol{\sigma}$ is given by,

[TABLE]

where $I_{d}$ is the identity on $\mathbb{R}^{d}$ and $Tr(\mathcal{E}\,\boldsymbol{u}^{0})$ is the trace of the strain. Here the second derivative $\boldsymbol{u}_{tt}^{0}$ is the time derivative in the sense of distributions of $\boldsymbol{u}^{0}_{t}$ and ${\rm div}\boldsymbol{\sigma}$ is the divergence of the stress tensor $\boldsymbol{\sigma}$ in the distributional sense. This shows that sharp fracture evolutions can be approached by spatially smooth evolutions in the limit of vanishing non locality.

6 Conclusions

In this article, we have presented a numerical analysis for class of nonlinear nonlocal peridynamic models. We have shown that the convergence rate applies, even when the fields do not have well-defined spatial derivatives. We treat both the forward Euler scheme as well as the general implicit single step method. The convergence rate is found to be the same for both schemes and is given by $C(\Delta t+h^{\gamma}/\epsilon^{2})$. Here the constant $C$ depends on $\epsilon$ and Hölder and $L^{2}$ norm of the solution and its time derivatives. The Lipschitz property of the nonlocal, nonlinear force together with boundedness of the nonlocal kernel plays an important role. It ensures that the error in the nonlocal force remains bounded when replacing the exact solution with its approximation. This, in turn, implies that even in the presence of mechanical instabilities the global approximation error remains controlled by the local truncation error in space and time.

Taking $\gamma=1$, a straight forward estimate of Equation 38 using Equation 12 gives to leading order

[TABLE]

where $C_{t}$ is independent of $\epsilon$ and depends explicitly on the $L^{2}$ norms of time derivatives of the solution see Equation 37 and

[TABLE]

[TABLE]

It is evident that the exponential factor could be large. However we can choose times $T$ for which the effects of the exponential factor can be diminished and $C_{1}$ and $C_{2}$ are not too large. To fix ideas consider a $1$ cubic meter sample and a corresponding 1400 meter per second shear wave speed. This wave speed is characteristic of plexiglass. Then the time for a shear wave to traverse the sample is $718$ $\mu$-seconds. This is the characteristic time $T^{\ast}$ and a fracture experiment can last a few hundred $\mu$-seconds. The actual time in $\mu$-seconds of a fracture simulation is given by $TT^{*}$ where $T$ is the non-dimensional simulation time. The dimensionless constant $\bar{C}$ is $1.19$ and we take $\epsilon=1/10$ and dimensionless body force unity. For a simulation cycle of length $TT^{\ast}=1.5\mu$-seconds the constants $C_{1}$ and $C_{2}$ in Equation 100 are $0.0193$ and $7.976$ respectively. The solution after cycle time $T$ can be used as initial conditions for a subsequent run and the process can be iterated. Unfortunately these estimates predict a total simulation time of $15\mu$-second before the relative error becomes greater than one even for a large number of spatial degrees of freedom. We point out that because the constants in the a-priori bound are necessarily pessimistic the predicted simulation time is an order of magnitude below what is seen in experiment. Future work will focus on a-posteriori estimates for simulations and adaptive implementations of the finite difference scheme.

In conclusion the analysis shows that the method is stable and one can control the error by choosing the time step and spatial discretization sufficiently small. However errors do accumulate with time steps and this limits the time interval for simulation. We have identified local perturbations for which the error accumulates with time step for the implicit Euler method. These unstable local perturbations correspond to regions for which a preponderance of bonds are in the softening regime.

Acknowledgements

RL would like to acknowledge the support and kind hospitality of the Hausdorff Institute for Mathematics in Bonn during the trimester program on multiscale problems.