Iterated Crank–Nicolson Method for Peridynamic Models

: In this paper, we explore the iterated Crank–Nicolson (ICN) algorithm for the one-dimensional peridynamic model. The peridynamic equation of motion is an integro-differential equation that governs structural deformations such as fractures. The ICN method was originally developed for hyperbolic advection equations. In peridynamics, we apply the ICN algorithm for temporal discretization and the midpoint quadrature method for spatial integration. Several numerical tests are carried out to evaluate the performance of the ICN method. In general, the ICN method demonstrates second-order accuracy, consistent with the Störmer–Verlet (SV) method. When the weight is 1/3, the ICN method behaves as a third-order Runge–Kutta method and maintains strong stability-preserving (SSP) properties for linear problems. Regarding energy conservation, the ICN algorithm maintains at least second-order accuracy, making it superior to the SV method, which converges linearly. Furthermore, selecting a weight of 0.25 results in fourth-order superconvergent energy variation for the ICN method. In this case, the ICN method exhibits energy variation similar to that of the fourth-order Runge–Kutta method but operates approximately 20% faster. Higher-order convergence for energy can also be achieved by increasing the number of iterations in the ICN method.


Introduction
Peridynamics, introduced by Stewart Silling in 2000 [1], is a nonlocal theory designed to model structural deformations such as fractures [2][3][4].Its applications extend beyond mechanics and material science to encompass fields like biology and multi-physics problems [5].The peridynamic equation of motion is an integro-differential equation governing nonlocal wave interaction.Unlike classical continuum mechanics, which relies on partial differential equations and faces challenges in handling discontinuities, the peridynamic model employs summation of bond forces over a finite region known as the horizon.This unique feature enables the effective handling of discontinuous solutions.Numerical methods for solving nonlocal wave equations involve temporal discretization of derivatives and numerical integration in the spatial domain.Various numerical algorithms have been adopted into peridynamics, including the meshfree method [6], quadrature methods [7,8], the finite element method [9], the discontinuous Galerkin method [10], and spectral methods [11][12][13][14].
The Störmer-Verlet (SV) algorithm [15] is a popular temporal discretization method for peridynamics, implemented in software packages like PDMATLAB2D [16].The SV method offers advantages such as time reversal and symplectic properties, leading to total energy conservation in a Hamiltonian system.However, as noted by [17], the magnitude of the energy variation being greater than zero is influenced by the choice of the numerical quadrature method.
The Crank-Nicolson (CN) algorithm [18] is another widely used numerical method capable of preserving the total energy of a Hamiltonian system.The CN method is secondorder accurate, implicit, and unconditionally stable.Because it is implicit, the discretized equations can be solved using iterations, leading to the explicit iterated Crank-Nicolson (ICN) algorithm [19].Originally developed for hyperbolic advection equations with applications in numerical relativity [20][21][22], the ICN method has been extended to other fields, including beam propagation equations [23,24] and Maxwell's equations for electromagnetic waves [25,26].
In this paper, we explore the ICN method for peridynamics.We employ the weighted ICN algorithm for temporal discretization and the midpoint quadrature method to evaluate the integral in the spatial domain.Teukolsky [20] discovered that the number of iterations does not affect the convergence rate of the ICN method for hyperbolic advection equations.As long as the number of iterations is greater than one, the method maintains second-order accuracy.Nevertheless, employing a large number of iterations can enhance its energy conservation property, as it gradually converges to the energy-conserving solution of the CN method.To enhance stability, the ICN method is generalized by introducing a weight θ, leading to the θ-ICN method [27].When the weight equals 0.5, the method reverts to the original ICN method.However, the θ-ICN method degrades to first-order accuracy when θ is not equal to 0.5.In another study [28], the θ-ICN method is refined to achieve secondorder accuracy by employing two different weights in consecutive iterations, with the geometric mean of the two weights set to 0.5.Our investigation of the ICN method for peridynamics aims to identify an efficient method with enhanced performance in terms of energy conservation.
The paper is organized as follows: In Section 2, we review the peridynamic equation and its corresponding numerical algorithms.In Section 3, we present the ICN method.Numerical examples are presented in Section 4, followed by the conclusion.

Peridynamic Model
In this section, we briefly review the peridynamic equation and the corresponding numerical methods.
Consider a one-dimensional homogeneous and infinitely long bar.The peridynamic equation of motion can be written as where u, ρ, b, and f represent the displacement field, the density, the external force, and the pairwise force function, respectively.The limits of the integral in (1) are usually replaced with finite numbers by ignoring the interaction between particles with a distance greater than the horizon δ, or equivalently, f = 0 when | x − x| > δ.Equation ( 1) can be written as and a system of first-order equations where is linear, then it can be written as A(x, u) = Āu, where Ā is a matrix that does not depend on u.For the linear case, the total energy of the system is given by [17] By viewing the system (3) and ( 4) as a Hamiltonian system, we see that E is the Hamiltonian (energy), so it is a conserved quantity.
The system of peridynamic Equations ( 3) and ( 4) can be written as a vector equation: where U = (u, v) T and F(U) = (v, A(x, u) + B) T .In the linear case, A(x, u) = Āu, and we obtain where A numerical algorithm for the system (3) and ( 4) involves the evaluation of the integral A(x, u) in the spatial domain and the temporal discretization of the derivatives u t and v t .Using the Störmer-Verlet (SV) algorithm, we obtain the following discretized equations: Another popular method is the Runge-Kutta (RK) algorithm.A general s-stage RK method for Equation ( 6) can be written in the following standard form: where a ij and b i are the coefficients in the Butcher tableau [29].In particular, when s = 4, a 21 = a 32 = 1/2, a 43 = 1, a 31 = a 41 = a 42 = 0, b 1 = b 4 = 1/6, and b 2 = b 3 = 1/3, we obtain the classical fourth-order RK (RK4) method [30] To evaluate the integral A(x, u), a commonly used method is the quadrature algorithm.For example, using the midpoint rule, we have where f j = f (x j − x, u(x j , t) − u(x, t)).In the linear case, using the midpoint method, the matrix Ā is symmetric and negative definite [17].Other Newton-Cotes quadrature rules include the composite trapezoidal rule and the composite Simpson's rule where N ≥ 3 must be an odd number for Simpson's rule.
In general, the SV method ( 9)-( 11) is second-order in time, and the midpoint quadrature is second-order in space.Consequently, the resulting midpoint SV method attains second-order accuracy in both space and time.The SV method exhibits time reversal and symplectic properties.Therefore, it preserves the total energy of the Hamiltonian system.
In contrast, although commonly used RK methods have orders larger than two, they do not preserve the total energy.Let E 0 represent the energy at the initial time t = 0 and E n the energy at time step t n .For an energy-conserving method, we should have |E n − E 0 | = 0 (or of the order of the machine round-off error), where |E n − E 0 | is the total energy variation.However, it has been reported in [17] that |E n − E 0 | > 0 for the SV method.

The Iterated Crank-Nicolson Method
In [19,20,27,28], the iterated Crank-Nicolson (ICN) and the weighted ICN algorithms are developed as predictor-corrector methods for hyperbolic partial differential equations.The peridynamic equation is a second-order integro-differential Equation ( 1), which can be expressed as a system of first-order differential Equations ( 3) and ( 4) and in its vector form (6).
In this section, we start with the Crank-Nicolson (CN) method for system of first-order differential equations and derive the ICN method as an iterative solver for the discretized equations from the CN algorithm.
There are two versions of the CN algorithms for the system (6): and where U n and U n+1 represent the solution at time t n = n∆t and t n+1 = (n + 1)∆t, respectively.If the function F is linear in terms of the variable U, then Equations ( 22) and ( 23) are equivalent.
The implicit nature of the method requires the solution of a system of linear or nonlinear equations.When solved by iterations, it leads to the ICN method.In [28], the ICN method is developed using a predictor-corrector strategy, and it can also be derived using the CN Equation (23).Here, we use the CN Equation ( 22) to demonstrate the ICN method, and the other one can be derived following a similar procedure.The iterative solver is a fixed-point iteration.Replacing U n+1 on the right-hand side of Equation ( 22) with the k-th approximation U k and letting the initial guess U 0 = U n , we obtain the following ICN solver: where m is the number of iterations.
In the ICN Equations ( 24)-( 26), if we use a weight θ when averaging F(U k ) and F(U n ) in Equation ( 25), then we obtain the θ-ICN method: For hyperbolic equations, the stability region of the θ-ICN method depends on the weight θ.As a result, by changing the weight θ, we can obtain more stable solutions.However, the θ-ICN method is only first-order accurate [27].In [28], the θ-ICN method is improved to second-order accuracy by employing different weights at consecutive steps where the number of iterations is fixed at two.The modified θ-ICN method can be written as and it is second-order accurate when θ = 0.25/θ.This method is referred as the geometric- averaging (GA) θ-ICN method as the geometric mean of θ and θ is 0.5.When θ = 0.5, it becomes the original ICN method.For simplicity, we call it the ICN method in the proceedings part of the paper.
The ICN algorithm ( 30)-( 32) is a three-stage RK method and can be written in standard form as: A special case is when θ = 1/3: and it is a third-order RK method as it satisfies the criteria given in Gottlieb and Shu's paper [31].Furthermore, if F is linear, this method can further be written as which makes it a strong stability-preserving (SSP) method [32].

Numerical Examples
In this section, we evaluate the performance of the ICN method for peridynamics through a series of numerical simulations.Our computer codes were implemented in MATLAB and executed on a 64-bit computer equipped with a 12th generation Intel(R) Core(TM) i7-12700 CPU @ 2.10GHz.
For simplicity, we set b(x, t) = 0 and ρ(x) = 1.The initial condition is u 0 (x) = exp(−x 2 ) and v 0 = 0 unless otherwise specified.We consider the following pairwise force function: where the micromodulus function C is given by: The problem is linear when r = 1 and nonlinear otherwise.We use a finite horizon δ, so that The computational spatial domain [L x , U x ] is divided into N grid cells so ∆x = (U x − L x )/N.We let the solution be defined at the center of each cell, so U n j = U(x j , n∆t) where x j = L x + (j − 1 2 )∆x, j = 1, 2, . . ., N. Similar to [17], a periodic boundary condition is implemented to approximate the original problem of an infinitely long bar.The simulation stops before the solution reaches the boundaries.
To test the rate of convergence, we compute the L 2 error norm and the where u e is the exact solution.
In our simulations, we primarily test the ICN methods with three weights: θ = 0.5, 1/3, and 0.25.The original ICN method corresponds to the case when θ = 0.5.The number of iterations for the ICN method is fixed at two unless specified.

Linear Peridynamic Equation
In the linear case where r = 1, the computational domain is chosen to be [−20, 20], and the mesh size is N.We let δ = 5 and ∆t = ∆x.The exact solution is given by [33] The solutions of the ICN method (θ = 0.5) at four time instances are shown in Figure 1.The results are very similar for other values of θ and the SV method.The grid size is N = 800.We can see that the numerical solutions agree with the exact solutions very well.Table 1 shows the rates of convergence of the solutions measured in the L 2 and L ∞ norms, together with the energy variation |E n − E 0 |.In terms of L 2 and L ∞ norms, the SV method converges quadratically, the same as the ICN methods with θ = 0.5 and 0.25.When θ = 1/3, the ICN method is third-order.
A symplectic method preserves the total energy, meaning the energy variation |E n − E 0 | = 0 is ideally zero or is on the order of the machine round-off error.However, as shown in Figure 2, the energy variation of the SV method is approximately 10 −3 , and as illustrated in Table 1, besides being zero, it converges linearly.In the case of the ICN method with θ = 0.5, the magnitude of energy variation grows and reaches approximately 10 −2 when t = 10.On the other hand, when θ = 0.25, the energy variation grows much more slowly, with its magnitude around 10 −6 , several orders of magnitude lower than the other two methods.From Table 1, we have the following observations regarding the convergence of the energy variation.Firstly, the ICN method is at least second-order accurate, surpassing the linear convergence of the SV method.Secondly, when θ = 0.5, both the energy and the solution exhibit second-order convergence.Thirdly, when θ = 1/3, the energy variation converges quadratically, despite its third-order solution.Lastly, when θ = 0.25, the ICN method achieves fourth-order superconvergence.Table 2 compares the SV, the ICN with θ = 0.25, and the classical RK4 methods.For RK4, both the solution and the energy variation have fourth-order convergence, making it one of the optimal choices for temporal discretization.The energy variation of the SV method is four or five orders of magnitude larger than the other two methods.The ICN method operates approximately 30% slower than the SV method due to the additional evaluation of integrals in the updating equations.The ICN and RK4 methods yield similar energy variation for both resolutions N = 1600 and N = 3200, while the ICN solver runs approximately 20% faster than the RK4 method at both resolutions.

Energy Conservation for the ICN Method with Different Weights
In this test, we investigate the relationship between the energy variation and the weight θ in the ICN method.The simulation setup remains consistent with the previous section, where the grid size is fixed at N = 800.The results presented in Table 3 demonstrate that as the weight θ decreases from 0.7 to 0.2, the energy variation E n − E 0 increases from −1.7 × 10 −2 to 1.9 × 10 −3 .Notably, when θ = 0.25, the energy variation reaches its minimum amplitude of 3.7 × 10 −6 .

The ICN Method with Multiple Iterations
In this example, we study the original ICN method (θ = 0.5) with a varying number of iterations.Firstly, we examine the relationship between the energy variation and the number of iterations in the ICN method.The simulation setup remains consistent with Section 4.1.We consider two grid sizes: N = 800 and 1600.As shown in Table 4, as m increases, the energy variation decreases.For N = 800 and N = 1600, the energy variations reach the order of the machine round-off error of 10 −15 when m is greater than 10 and 8, respectively.This example illustrates that the ICN method can preserve energy up to machine precision by employing a large number of iterations.Furthermore, when the grid size increases, the energy variation reaches machine precision at a faster rate.For N = 1600, achieving the order of 10 −7 requires 3 iterations, with a CPU run time of approximately 37 s.This runtime is not substantially greater than that required for the ICN method with θ = 0.25 (26 s) and the RK4 method (33 s) in the previous test, as shown in Table 2. From the previous multi-iteration test, we observe that when the number of iterations is three, the original ICN method exhibits an energy variation at the same order of magnitude as the RK4 method.This motivates us to conduct a convergence test for the three-iteration ICN method.The results are presented in Table 5.We find that the three-iteration ICN method demonstrates second-order convergence in terms of the L 2 and L ∞ error norms, and fourth-order convergence in terms of energy conservation.Since both the three-iteration ICN and the RK4 methods are four-stage methods, they have comparable CPU run times.
Furthermore, we conduct a series of simulations to examine the rate of convergence of the ICN method with varying numbers of iterations, as shown in Table 6.We observe that as the number of iterations increases, the rate of convergence also increases.When the number of iterations is an even number, the rate of convergence remains the same as the number of iterations.However, for an odd number of iterations m, the rate of convergence is m + 1.These results indicate that a higher order of convergence (in terms of energy conservation) can be achieved by increasing the number of iterations in the ICN method.As shown in [17], the choice of quadrature formula for integration impacts the accuracy of energy conservation.In this example, we examine three different Newton-Cotes quadrature formulas: the midpoint rule (Equation ( 19)), the composite trapezoidal rule (Equation ( 20)), and the composite Simpson's rule (Equation ( 21)).For Simpson's rule, we adjust the grid size to N + 1, as N is an even number in our simulations.Results are presented in Tables 7 and 8.We do not observe significant differences among the three quadrature formulas.All three methods exhibit the same order of convergence.The L 2 , L ∞ , and energy variation are of comparable magnitudes across all three methods, with the differences being negligible.This example suggests that the degree of the Newton-Cotes quadrature formulas does not significantly affect the convergence rate.

Table 7.
Comparison of three quadrature formulas for the peridynamic simulation using the SV method.ϵ L 2 and ϵ L ∞ represent the numerical errors measured in the L 2 and L ∞ norms, respectively.|E n − E 0 | denotes the energy variation from t = 0 to t = 10.

Linear Peridynamic Model with Discontinuous Initial Condition
In this example, we explore our method's capability for handling discontinuities, focusing on discontinuous initial conditions [34,35].The simulation setup mirrors the previous linear example in Section 4.1, except for the initial condition.
First, we utilize the following piecewise constant function as the initial displacement field u 0 : The initial velocity field v 0 is set to zero.The solutions of the ICN method (θ = 0.25) at four time instances are shown in Figure 3, and the convergence rates are listed in Table 9.
For the second experiment in this section, we apply a piecewise constant function as the initial velocity field v 0 : The initial displacement field is a continuous function u 0 (x) = exp(−x 2 ).The results are presented in Figure 4 and in Table 10.When computing the convergence rates, we utilize the solution obtained on a fine mesh (N = 3200) as the reference solution.In both scenarios, the convergence rates are akin to those observed with continuous initial conditions.Specifically, with regard to the L 2 and L ∞ error norms, both the SV method and the ICN method demonstrate quadratic convergence.In terms of energy conservation, the SV method displays linear convergence, whereas the ICN method with θ = 0.25 achieves fourth-order accuracy.

Nonlinear Peridynamic Equation
For the nonlinear case, we set r = 3 in Equation (43).The computational domain spans [−10, 10], with ∆t = ∆x.The solution obtained through the ICN method is presented in Figure 5, and Table 11 illustrates the convergence rates of these solutions.We use the solution on a fine mesh (N = 3200) as the exact solution when computing the convergence rates.From the tabulated data, it is evident that when θ = 1/3, the ICN method demon-strates third-order convergence.For the other two values of θ, the ICN method converges quadratically, similar to the SV method.

Figure 2 .
Figure 2. Time history of the total energy variation E(t) − E(0).The grid size is N = 800.

Table 1 .
Comparison of the SV method and ICN methods with various weights for the linear peridynamic equation.ϵ L 2 and ϵ L ∞ represent the numerical errors measured in the L 2 and L ∞ norms, respectively.|E n − E 0 | denotes the energy variation from t = 0 to t = 10.

Table 2 .
Comparison of energy variation and run time for three methods (SV, ICN, and RK4).

Table 3 .
Energy variation of the ICN method with different weights.

Table 4 .
Energy variation of the ICN method when the number of iterations changes.E n is the total energy at t = 10.

Table 6 .
Energy convergence rate as a function of the number of iterations in the ICN method.
4.4.Comparison of Three Quadrature Formulas for Linear Peridynamics

Table 9 .
Numerical errors and convergence rates of the SV and the ICN methods for the linear peridynamic equation with piecewise constant initial displacement field (48).ϵ L 2 and ϵ L ∞ represent the numerical errors measured in the L 2 and L ∞ norms, respectively.|E n − E 0 | denotes the energy variation from t = 0 to t = 10.

Table 10 .
Numerical errors and convergence rates of the SV and the ICN methods for the linear peridynamic equation with piecewise constant initial velocity field (49).ϵ L 2 and ϵ L ∞ represent the numerical errors measured in the L 2 and L ∞ norms, respectively.|E n − E 0 | denotes the energy variation from t = 0 to t = 10.

Table 11 .
Comparison of the SV method and the ICN methods with various weights for the nonlinear peridynamic equation.ϵ L 2 and ϵ L ∞ represent the numerical errors measured in the L 2 and L ∞ norms, respectively.