Finite element convergence for state-based peridynamic fracture modelsFunding: This material is based upon work supported by the U. S. Army Research Laboratory and the U. S. Army Research Office under contract/grant number W911NF1610456.

# Finite element convergence for state-based peridynamic fracture models††thanks: Funding: This material is based upon work supported by the U. S. Army Research Laboratory and the U. S. Army Research Office under contract/grant number W911NF1610456.

Prashant K. Jha Department of Mathematics, Louisiana State University, Baton Rouge, LA (prashant.j16o@gmail.com    Robert Lipton Department of Mathematics, Louisiana State University, Baton Rouge, LA (lipton@math.lsu.edu
###### Abstract

We establish the a-priori convergence rate for finite element approximations of a class of nonlocal nonlinear fracture models. We consider state based peridynamic models where the force at a material point is due to both the strain between two points and the change in volume inside the domain of nonlocal interaction. The pairwise interactions between points are mediated by a bond potential of multi-well type while multi point interactions are associated with volume change mediated by a hydrostatic strain potential. The hydrostatic potential can either be a quadratic function, delivering a linear force-strain relation, or a multi-well type that can be associated with material degradation and cavitation. We first show the well-posedness of the peridynamic formulation and that peridynamic evolutions exist in the Sobolev space . We show that the finite element approximations converge to the solutions uniformly as measured in the mean square norm. For linear continuous finite elements the convergence rate is shown to be , where is the size of horizon, is the mesh size, and is the size of time step. The constants and are independent of and and may depend on through the norm of the exact solution. We demonstrate the stability of the semi-discrete approximation. The stability of the fully discrete approximation is shown for the linearized peridynamic force. We present numerical simulations with dynamic crack propagation that support the theoretical convergence rate.

Key words. Nonlocal fracture models, peridynamic, state based peridynamic, numerical analysis, finite element approximation

AMS subject classifications. 34A34, 34B10, 74H55, 74S05

## 1 Introduction

In this work, we study the state-based peridynamic theory and obtain an a-priori error bound for the finite element approximation. The peridynamic theory is a reformulation of classical continuum mechanics carried out in the work of Silling in [Silling, 2000, Silling et al., 2007]. The strain inside the medium is expressed in terms of displacement differences as opposed to the displacement gradients. Acceleration of a point is now due to the sum of the forces acting on the point from near by points. The new kinematics bypasses the difficulty incurred by juxtaposing displacement gradients and discontinuities as in the case of classical fracture theories. The nonlocal fracture theory has been applied numerically to model the complex fracture phenomenon in materials, see [Weckner and Abeyaratne, 2005], [Silling and Bobaru, 2005], [Silling and Lehoucq, 2008],
[Silling et al., 2010], [Foster et al., 2011], [Ha and Bobaru, 2010], [Agwai et al., 2011],
[Bobaru and Hu, 2012], [Ghajari et al., 2014], [Du et al., 2013a], [Lipton et al., 2016]. Every point interacts with its neighbors inside a ball of fixed radius called the horizon. The size of horizon sets the length scale of nonlocal interaction. When the forces between points are linear and nonlocal length scale tends to zero, it is seen that peridynamic models converge to the classic model of linear elasticity, [Emmrich et al., 2013, Silling and Lehoucq, 2008, Aksoylu and Unlu, 2014, Mengesha and Du, 2015]. The work of [Tian and Du, 2014] provides an analytic framework for analyzing FEM for linear bond and state-based peridynamics. For nonlinear forces associated with double well potentials the peridynamic evolution converges in the small horizon limit to an evolution with a sharp evolving fracture set and the evolution is governed by the classic linear elastic wave equation away from the fracture set, see [Lipton, 2016, Lipton, 2014, Jha and Lipton, 2018a]. A recent review of the state of the art can be found in [Bobaru et al., 2016] and [Du, 2018a].

In this work we assume small deformation and work with linearized bond-strain. Let , for , be the material domain. For a displacement field , the bond-strain between two material points is given by

 S(y,x,t;u)=u(y,t)−u(x,t)|y−x|⋅y−x|y−x|. (1)

Let be the size of horizon and be the neighborhood of a material point . For pairwise interaction, we assume the following form of pairwise interaction potential

 Wϵ(S(y,x,t;u))=Jϵ(|y−x|)ϵ|y−x|f(√|y−x|S(y,x,t;u)), (2)

where is the influence function. We assume where for and for . The potential , see LABEL:fig:f_a, is assumed to be convex for small strains and becomes concave for larger strains. In the widely used prototypical micro-elastic brittle (PMB) peridynamic material, the strain vs force profile is linear up to some critical strain and is zero for any strain above . In contrast, the peridynamic force given by , is linear near zero strain and as the strain gets larger and reaches the critical strain, () for positive (negative) strain, the bond starts to soften, see LABEL:fig:f_b. For a given potential function , the critical strain is given by and where are the inflection points of the potential function as shown in LABEL:fig:f_a.

The spherical or hydrostatic strain at material point is given by

 θ(x,t;u)=1ϵdωd∫Hϵ(x)Jϵ(|y−x|)S(y,x,t;u)|y−x|dy, (3)

where is the volume of unit ball in dimension . The potential for hydrostatic interaction is of the form

 Vϵ(θ(x,t;u))=g(θ(x,t;u))ϵ2, (4)

where is the potential function associated to hydrostatic strain. Here can be of two types: 1) a quadratic function with only one well at zero strain, and 2) a convex-concave function with a wells at the origin and at , see LABEL:fig:g_a. If is assumed to be quadratic then the force due to spherical strain is linear. If is a multi-well potential, the material softens as the hydrostatic strains exceeds critical value. For the convex-concave type , the critical values are and beyond which the force begins to soften is related to the inflection point and of as follows

 θ+c=r+∗,θ−c=r−∗. (5)

The critical compressive hydrostatic strain where the force begins to soften for negative hydrostatic strain is chosen much larger in magnitude than , i.e. .

The finite element approximation has been applied to peridynamic fracture, however there remains a paucity of literature addressing the rigorous a-priori convergence rate of the finite element approximation to peridynamic problems in the presence of material failure. This aspect provides the motivation for the present work. In this paper we first prove existence of peridynamic evolutions taking values in that are twice differentiable in time, see LABEL:thm:existence. We note that as these evolutions will become more fracture like as the region of nonlocal interaction decreases. These evolutions can be thought of as inner approximations to fracture evolutions. On passing to subsequences it is possible to show that the evolutions converge in the limit of vanishing non-locality to a limit solution taking values in the space of special functions of bounded deformation SBD. Here the limit evolution has a well defined Griffith fracture energy bounded by the initial data, see [Lipton, 2016] and [Jha and Lipton, 2019]. We show here that higher temporal regularity can be established if the body force changes smoothly in time. Motivated by these considerations we develop finite element error estimates for solutions that take values in and for a bounded time interval.

In this paper we obtain an a-priori error bound for the finite element approximation of the displacement and velocity using a central in time discretization. Due to the nonlinear nature of the problem we get a convergence rate by using Lax Richtmyer stability together with consistency. Both stability and consistency are shown to follow from the Lipschitz continuity of the peridynamic force in , see section LABEL:consistent and section LABEL:stable. The bound on the error is uniform in time and is given by , where the constants and are independent of and mesh size , see LABEL:thm:convergence. A more elaborate discussion of the a-priori bound is presented in section LABEL:s:convergence. For the linearized model we obtain a stability condition on , LABEL:thm:cflcondition, that is of the same form as those given for linear local and nonlocal wave equations [Karaa, 2012, Guan and Gunzburger, 2015]. We demonstrate stability for the linearized model noting that for small strains the material behaves like a linear elastic material and that the stability of the linearized model is necessary for the stability of nonlinear model. We believe a more constructive CFL stability condition is possible for the linear case and will pursue this in future work.

Previous work [Jha and Lipton, 2018a] treated spatially Lipschitz continuous solutions and addressed the finite difference approximation and obtained bounds on the error for the displacement and velocity that are uniform in time and of the form , where constants the and are as before. For finite elements the convergence rate is seen to be slower than for the FEM model introduced here and is of order as opposed to . On the other hand the FEM method increases the computational work due to the inversion of the mass matrix.

We carry out numerical experiments for dynamic crack propagation and obtain convergence rates for Plexiglass that are in line with the theory, see LABEL:s:numerical. We also compare the Griffith’s fracture energy with the peridynamic energy of the material softening zone we show good agreement between the two energies, see LABEL:sec:_softzone. Finite difference methods are less expensive than finite element approximations for nonlocal problems, however the latter offers more control on the accuracy of solution see, [Macek and Silling, 2007, Littlewood, 2010, Du et al., 2013c, Gerstle et al., 2007, Du, 2018b].

Here the a-priori convergence rates for the FEM given by LABEL:thm:convergence include the effects of material degradation through the softening of material properties. The FEM simulations presented in this paper show that the material develops localized softening zones (region where bonds exceed critical tensile strain) as it deforms. This is in contrast to linear peridynamic models which are incapable of developing softening zones. For nonlinear peridynamic models with material failure the localization of zones of softening and damage is the hallmark of peridynamic modeling [Silling, 2000, Silling et al., 2007], [Ha and Bobaru, 2010], [Foster et al., 2011]. One notes that the a-priori error involves in the denominator and in many cases is chosen small. However typical dynamic fracture experiments last only hundreds of microseconds and the a-priori error is controlled by the product of simulation time multiplied by . So for material properties characteristic of Plexiglass and of size , the a-priori estimates predict a relative error of for simulations lasting around 100 microseconds. We point out that the a-priori error estimates assume the appearance of nonlinearity anywhere in the computational domain. On the other hand numerical simulation and independent theoretical estimates show that the nonlinearity concentrates along “fat” cracks of finite length and width equal to , see [Lipton, 2016, Lipton, 2014]. Moreover the remainder of the computational domain is seen to behave linearly and to leading order can be modeled as a linear elastic material up to an error proportional to , see [Proposition 6, [Jha and Lipton 2018b]]. Future work will use these observations to focus on adaptive implementation and a-posteriori estimates. A-posteriori convergence for FEM models of peridynamics with material degradation can be seen in the work [Macek and Silling, 2007], [Chen and Gunzburger, 2011],
[Ren et al., 2017]. For other nonlinear and nonlocal models adaptive mesh refinement within FE framework for nonlocal models has been explored in [Du et al., 2013c] and convergence of the adaptive FE approximation is rigorously shown. A-posteriori error analysis of linear nonlocal models is carried out in [Du et al., 2013b].

The paper is organized as follows. We introduce the equation of motion in LABEL:s:nonlocal_dynamics and present the Lipschitz continuity of the force, existence of peridynamic solution, and the higher temporal regularity necessary for the finite element error analysis. In LABEL:s:fem, we consider the finite element discretization. We prove the stability of a semi-discrete approximation in LABEL:ss:semifem. In LABEL:s:fullfem, we analyze the fully discrete approximation and obtain an a-priori bound on errors. The stability of the fully discrete approximation linearized peridynamic force is shown in LABEL:ss:fullfemstability. We present our numerical experiments in LABEL:s:numerical. Proofs of the Lipschitz bound on the peridynamic force and higher temporal regularity of solutions is provided in LABEL:s:proofs. In LABEL:s:conclusions we present our conclusions.

We conclude the introduction by listing the notation used throughout the paper. We denote material domain as where for . Points and vectors in are denoted as bold letters. Some of the key notations are as follows

[0,T] Time domain Size of horizon Density Horizon of x∈D, a ball of radius ϵ centered at x Volume of unit ball in dimension d=2,3 Boundary function defined on D taking value 1 in the interior and smoothly decaying to 0 as x approaches ∂D Displacement field defined over D×[0,T]. We may also use notation u to denote field defined over just D Initial condition on displacement Body force defined over D×[0,T] The unit vector pointing from a point y to the point x Bond strain S=u(y,t)−u(x,t)|y−x|⋅ey−x. We may also use S(y,x;u) if u is a filed defined over just D Spherical or hydrostatic strain. We may also use θ(x;u) if u is a filed defined over just D Critical bond strain Critical hydrostatic strain Influence function where J is integrable with J(r)=0 for r≥1 and 0≤J(r)≤M for r<1 Moment of function J over H1(0) with weight 1/(ωd|ξ|α) Potential functions for pairwise and state-based interaction Pairwise and state-based potential energy density Total peridynamic potential energy at time t Total dynamic energy at time t total peridynamic force, pairwise peridynamic force, and state-based peridynamic force respectively Nonlinear operator where u,v are vector fields over D Nonlinear pairwise and state-based operator L2 norm over D, L∞ norm over D, and Sobolev Hn norm over D (for n=1,2) respectively Size of mesh and size of time step Triangulation of D given by triangular/tetrahedral elements Continuous piecewise linear interpolation operator on Th Space of functions in H2(D;Rd) such that trace of function is zero on boundary ∂D, i.e. W=H2(D;Rd)∩H10(D;Rd) Space of continuous piecewise linear interpolations on Th Interpolation function of mesh node i Finite element projection of u onto Vh Total error in mean square norm at time step k Approximate displacement and velocity field at time step k Exact displacement and velocity field at time step k

## 2 Equation of motion, existence, uniqueness, and higher regularity

We assume to be an open set with boundary. To enforce zero displacement boundary conditions at and to insure a well posed evolution we introduce the boundary function . This function is introduced as a factor into the potentials and . Here the boundary function takes value in the interior of domain and is zero on the boundary. We assume and in our analysis. The hydrostatic strain is modified to include the boundary and is given by

 θ(x,t;u)=1ϵdωd∫Hϵ(x)ω(y)Jϵ(|y−x|)S(y,x,t;u)|y−x|dy. (6)

The peridynamic potentials LABEL:eq:tensilepot and LABEL:eq:hydropot are modified to see the boundary follows

 Wϵ(S(y,x,t;u)) =ω(x)ω(y)Jϵ(|y−x|)ϵ|y−x|f(√|y−x|S(y,x,t;u)), (7) Vϵ(θ(x,t;u)) =ω(x)g(θ(x,t;u))ϵ2. (8)

We assume that potential function is at least 4 times differentiable and satisfies following regularity condition:

 Cf0:=supr|f(r)|<∞,Cfi:=supr|f(i)(r)|<∞,∀i=1,2,3,4. (9)

If potential function is convex-concave type then we assume that satisfies same regularity condition as . We denote constants , for , similar to above.

The total potential energy at time is given by

 PDϵ(u(t))=1ϵdωd∫D∫Hϵ(x)|y−x|Wϵ(S(y,x,t;u))dydx (10) +∫DVϵ(θ(x,t;u))dx,

where potential and are described above. The material is assumed to be homogeneous and the density is given by . The applied body force is denoted by . We define the Lagrangian

 L(u,∂tu,t)=ρ2||˙u||2−PDϵ(u(t))+∫Db(t)⋅u(t)dx,

here is the velocity given by the time derivative of . Applying the principal of least action together with a straight forward calculation (see for example [Lipton et al., 2018a] for detailed derivation) gives the nonlocal dynamics

 ρ¨u(x,t)=Lϵ(u)(x,t)+b(x,t), % for x∈D, (11)

where

 (12)

is the peridynamic force due to bond-based interaction and is given by

 LϵT(u)(x,t) =2ϵdωd∫Hϵ(x)ω(x)ω(y)Jϵ(|y−x|)ϵ|y−x|∂Sf(√|y−x|S(y,x,t;u))ey−xdy, (13)

and is the peridynamic force due to state-based interaction and is given by

 LϵD(u)(x,t) =1ϵdωd∫Hϵ(x)ω(x)ω(y)Jϵ(|y−x|)ϵ2[∂θg(θ(y,t;u))+∂θg(θ(x,t;u))]ey−xdy. (14)

The dynamics is complemented with the initial data

 u(x,0)=u0(x),∂tu(x,0)=v0(x). (15)

We prescribe zero Dirichlet boundary condition on the boundary

 u(x)=0∀x∈∂D. (16)

We extend the zero boundary condition outside to whole . In our analysis we will assume mass density without loss of generality.

### 2.1 Existence of solutions and higher regularity in time

We recall that the space is the closure in the norm of the functions that are infinitely differentiable with compact support in . For suitable initial conditions and body force we show that solutions exist in

 W=H2(D;Rd)∩H10(D;Rd)={v∈H2(D;Rd):γv=0, on ∂D} (17)

where is the trace of the function on the boundary of . We will assume that is extended by zero outside . We first exhibit the Lipschitz continuity property and boundedness of the peridynamic force for displacements in . We will then apply [Theorem 3.2, [Jha and Lipton, 2017]] to conclude the existence of unique solutions.

We note the following Sobolev embedding properties of when is a domain.

• From Theorem 2.72 of [Demengel, 2012], there exists a constant independent of such that

 ||u||∞≤Ce1||u||2. (18)
• Further application of standard embedding theorems
(e.g., Theorem 2.72 of [Demengel, 2012]), shows there exists a constant independent of such that

 ||∇u||Lq(D;Rd×d)≤Ce2||∇u||1≤Ce2||u||2, (19)

for any such that when and when .

We have following result which show the Lipschitz continuity property of a peridynamic force .

###### Theorem 1.

Lipschitz continuity of peridynamic force
Let be a convex-concave function satisfying for and let either be a quadratic function, or be a convex-concave function with for . Also, let boundary function be such that and . Then, for any , we have

 ||Lϵ(u)−Lϵ(v)||2 ≤¯L1(1+||u||2+||v||2)2ϵ3||u−v||2, (20)

where constant does not depend on and . Also, for , we have

 ||Lϵ(u)||2 ≤¯L2(||u||2+||u||22)ϵ5/2, (21)

where constant does not depend on and .

Now let be any positive number, a straight-forward application of [Theorem 3.2, [Jha and Lipton, 2017]] gives:

###### Theorem 2.

Existence and uniqueness of solutions over finite time intervals
Let , , and satisfy the hypothesis of LABEL:thm:lipschitzproperty. For any initial condition , time interval , and right hand side continuous in time for such that satisfies , there is a unique solution of peridynamic equation LABEL:eq:equationofmotion. Also, and are Lipschitz continuous in time for .

We can also show higher regularity in time of evolutions under suitable assumptions on the body force:

###### Theorem 3.

Higher regularity
Suppose the initial data and righthand side satisfy the hypothesis of LABEL:thm:existence and suppose further that exists and is continuous in time for and  . Then and

 ||∂3tttu(x,t)||2 ≤C(1+sups∈I0||u(s)||2)2ϵ3sups∈I0||∂tu(s)||2+||˙b(x,t)||2, (22)

where is a positive constant independent of .

The proofs of Theorems LABEL:thm:lipschitzproperty, and LABEL:thm:higherregularity are given in LABEL:s:proofs. For future reference, we note that for any , we have

 ||Lϵ(u)−Lϵ(v)|| ≤Lϵ2||u−v||. (23)

where constant is given by

 L :={4(Cf2¯J1+Cg2¯J20)if g is a convex-concave type,4(Cf2¯J1+g′′(0)¯J20)if g % is a quadratic function, (24)

and .

### 2.2 Weak form

We multiply LABEL:eq:equationofmotion by a test function in and integrate over to get

 (¨u(t),~u)=(Lϵ(u(t)),~u)+(b(t),~u). (25)

We have the following integration by parts formula:

###### Lemma 4.

For any we have

 (Lϵ(u),v)=−aϵ(u,v), (26)

where

 aϵ(u,v)=aϵT(u,v)+aϵD(u,v) (27)

and

 aϵT(u,v) =1ϵd+1ωd∫D∫Dω(x)ω(y)Jϵ(|y−x|) ∂Sf(√|y−x|S(y,x;u))S(y,x;v)dydx, aϵD(u,v) =1ϵ2∫Dω(x)g′(θ(x;u))θ(x;v)dx. (28)

The proof of above lemma is identical to the proof of Lemma 4.2 in [Lipton et al., 2018a].

Using the above Lemma, the weak form of the peridynamic evolution is given by

 (¨u(t),~u)+aϵ(u(t),~u)=(b(t),~u). (29)

Total dynamic energy: We define the total dynamic energy as follows

 Eϵ(u)(t)=12||˙u(t)||2L2+PDϵ(u(t)), (30)

where is defined in LABEL:eq:totalenergy. Time derivative of total energy satisfies

 ddtEϵ(u)(t)=(¨u(t),˙u(t))+aϵ(u(t),˙u(t)). (31)

Remark. It is readily verified that the peridynamic force and energy are bounded for all functions in . Here the bound on the force follows from the Lipschitz property of the force in , see, LABEL:eq:lipschitzproperty_l2. The peridynamic force is also bounded for functions in . This again follows from the Lipschitz property of the force in using arguments established in LABEL:s:proofs. The boundedness of the energy in both and follows from the boundedness of the bond potential energy and used in the definition of , see LABEL:eq:tensilepotmod and LABEL:eq:hydropotmod. More generally this also shows that for .

We next discuss the spatial and the time discretization of peridynamic equation.

## 3 Finite element approximation

Let be given by linear continuous interpolations over tetrahedral or triangular elements where denotes the size of finite element mesh. Here we assume the elements are conforming and the finite element mesh is shape regular and .

For a continuous function on , is the continuous piecewise linear interpolant on . It is given by

 Ih(u)∣∣∣T=IT(u)∀T∈Th, (32)

where is the local interpolant defined over finite element and is given by

 IT(u)=n∑i=1u(xi)ϕi. (33)

Here is the number of vertices in an element , is the position of vertex , and is the linear interpolant associated to vertex .

Application of Theorem 4.4.20 and remark 4.4.27 in [Brenner and Scott, 2007] gives

 ||u−Ih(u)||≤ch2||u||2, ∀u∈W. (34)

Let denote the projection of on . For the norm it is defined as

 ||u−rh(u)|| =inf~u∈Vh||u−~u||. (35)

and satisfies

 (rh(u),~u)=(u,~u),∀~u∈Vh. (36)

Since , and LABEL:eq:interpolationerror we see that

 ||u−rh(u)||≤ch2||u||2, ∀u∈W (37)

### 3.1 Semi-discrete approximation

Let be the approximation of satisfying following for all

 (¨uh,~u)+aϵ(uh(t),~u) =(b(t),~u),∀~u∈Vh. (38)

We have following result:

###### Theorem 5.

Energy stability of semi-discrete approximation
The semi-discrete scheme is stable and the energy , defined in LABEL:eq:totaldynamicenergy, satisfies the following bound

 Eϵ(uh)(t) ≤[√Eϵ(uh)(0)+∫t0||b(τ)||dτ]2.

We note that while proving the stability of semi-discrete scheme corresponding to nonlinear peridynamics, we do not require any assumption on the strain . The proof is similar to [Section 6.2, [Lipton, 2016]].

###### Proof.

Letting in LABEL:eq:feweak and noting the identity LABEL:eq:energyidentity, we get

 ddtEϵ(uh)(t)=(b(t),˙uh(t))≤||b(t)||||˙uh(t)||. (39)

We also have

 ||˙uh(t)||≤2√12||˙uh||2+PDϵ(uh(t))=2√Eϵ(uh)(t)

where we used the fact that is nonnegative. We substitute above inequality in LABEL:eq:ineq_stab_1 to get

 ddtEϵ(uh)(t) ≤2√Eϵ(uh)(t)||b(t)||.

We fix and define as . Then, from the above equation we easily have

 ddtA(t) ≤2√A(t)||b(t)||⇒12ddtA(t)√A(t)≤||b(t)||.

Noting that , integrating from to and relabeling as , we get

 √A(t) ≤√A(0)+∫t0||b(s)||ds.

Proof is complete once we let and take the square of both sides.

## 4 Central difference time discretization

In LABEL:s:convergence we calculate the convergence rate for the central difference time discretization of the fully nonlinear problem. We then present a CFL like condition on the time step for the linearized peridynamic equation in LABEL:ss:fullfemstability.

At time step , the exact solution is given by where , and their projection onto is given by . The solution of fully discrete problem at time step is given by .

We approximate the initial data on displacement and velocity by their projections and . Let and . For , satisfies, for all ,

 (uk+1h−ukhΔt,~u) =(vk+1h,~u), (vk+1h−vkhΔt,~u) =(Lϵ(ukh),~u)+(bkh,~u), (40)

where we have denoted the projection of , i.e. , as . Combining the two equations delivers central difference equation for . We have

 (uk+1h−2ukh+uk−1hΔt2,~u) =(Lϵ(ukh),~u)+(bkh,~u),∀~u∈Vh. (41)

For , we have

 (u1h−u0hΔt2,~u) =12(Lϵ(u0h),~u)+1Δt(v0h,~u)+12(b0h,~u). (42)

### 4.1 Implementation details

For completeness we describe the implementation of the time stepping method using FEM interpolants. Let be the shape tensor then are given by

 ukh=NUk,~u=N~U, (43)

where and are dimensional vectors, where is the number of nodal points in the mesh and is the dimension.

From LABEL:eq:central, for all with elements of zero on the boundary, then the following holds for

 [MUk+1−2Uk+Uk−1Δt2]⋅~U=Fk⋅~U. (44)

Here the mass matrix and force vector are given by

 M :=∫DNTNdx, Fk :=Fkpd+∫DNTb(x,tk)dx, (45)

where is defined by

 Fkpd :=∫DNT(Lϵ(ukh)(x))dx. (46)

We remark that a similar equation holds for .

At the time step we must invert to solve for using

 Uk+1=Δt2M−1Fk+2Uk−Uk−1. (47)

As is well known this inversion amounts to an increase of computational complexity associated with discrete approximation of the weak formulation of the evolution. Further, the matrix-vector multiplication needs to be carried out at each time step. On the other hand the quadrature error in the computation of the force vector is reduced when using the weak form.

We next show the convergence of approximation.

### 4.2 Convergence of approximation

In this section, we prove the uniform bound on the error and show that the approximate solution converges to the exact solution with rate given by . Here horizon is assumed to be fixed. We first compare the exact solution with its projection in and then compare the projection with the approximate solution. We further divide the calculation of error between the projection and the approximate solution in two parts, namely consistency analysis and error analysis.

Error is given by

 Ek:=||ukh−u(tk)||+||vkh−v(tk)||.

The error is split into two parts as follows

 Ek ≤(||uk−rh(uk)||+||vk−rh(vk)||)+(||rh(uk)−ukh||+||rh(vk)−vk