On nonlinear viscoelastic deformations: a reappraisal of Fung's quasi-linear viscoelastic model

This paper offers a reappraisal of Fung's model for quasi-linear viscoelasticity. It is shown that a number of negative features exhibited in other works, commonly attributed to the Fung approach, are merely a consequence of the way it has been applied. The approach outlined herein is shown to yield improved behaviour and offers a straightforward scheme for solving a wide range of models. Results from the new model are contrasted with those in the literature for the case of uniaxial elongation of a bar: for an imposed stretch of an incompressible bar and for an imposed load. In the latter case, a numerical solution to a Volterra integral equation is required to obtain the results. This is achieved by a high-order discretization scheme. Finally, the stretch of a compressible viscoelastic bar is determined for two distinct materials: Horgan–Murphy and Gent.


Introduction
The study of nonlinear viscoelastic deformations of solid materials has a very long history, with a consequent proliferation of a diverse and extensive array of constitutive models. For a comprehensive overview of the topic, the reader is directed to the recent paper by Wineman [1] and related references therein. The subject is still active, with models continuing to be developed across the field, from highly mathematical approaches where implementation is not a concern to very applied studies where ease of application is essential. Each approach has its advantages: the models that more accurately describe the microphysics this relation must, at the minimum, satisfy objectivity, and hence the form chosen in (1.1). Nevertheless, a number of authors have erroneously expressed Fung's relation not for the second Piola-Kirchhoff stress, which guarantees objectivity, but in terms of other stress measures. For example, see Fung's original discussion on the subject in [2, page 253], in which the QLV relation is expressed in terms of the Kirchhoff stress! In the latest version [6.12] of ABAQUS, although the viscoelastic constitutive model is objective, the model is rather heuristic and as we shall see later, it assumes the same relaxation behaviour in both the shear and bulk parts of an isotropic material. Further technical aspects of objectivity are discussed in the paper by Liu [11].
Another approach to QLV sometimes employed in the literature (e.g. [12][13][14]) is, in an analogous fashion to nonlinear elasticity, to write the stress in terms of the instantaneous derivative (with respect to the principal stretch) of a strain energy function. Clearly, in this case the strain energy function must be dissipative and depend on the history of the strain, and so these authors express it as a fading-memory integral with integrand given as a hyperelastic strain energy function. This approach may be somewhat hard to justify and does not allow the user to consider the problem in terms of an auxiliary instantaneous measure of strain, i.e. the effective elastic stress term Π e given above.
Perhaps, a more subtle point to those mentioned above is the choice of Π e in (1.1). In this paper, the reasonable assertion is made that Π e must be zero whenever the deformation is zero. However, this point is often not recognized, especially for incompressible materials, as can be seen, for example, by inspection of the integrands (setting the stretch equal to unity) in the articles by [15] (see eqns (8) and (12) therein), [16] (eqns (9) and (13)) and [17] (eqn (9)). Note that this requirement is satisfied in the one-dimensional models in [8,9] but not in [10]. If Π e does not vanish for zero deformation then, in general, it can be expected that the actual stress field Π will be non-zero prior to the application of the imposed load or stretch, i.e. the solution will be noncausal! However, for incompressible materials, an arbitrary pressure term (Lagrange multiplier) can be adjusted so that the stress field is indeed causal, but then it must take a specific form at later times. That is, whenever the material is deformed, and then returns back to zero, the stress Π will instantaneously return to zero too. This is clearly an unrealistic consequence of the particular choice of model, for, as Fung states [2, page 228]: 'the tensile stress at any time t is equal to the instantaneous [elastic] stress response . . . decreased by an amount depending on the past history'. This issue is examined further in §4a.
There are other variants of QLV employed in the literature that offer slightly modified forms of the governing equation to that suggested by Fung. The reader is referred, for example, to articles [18][19][20]. In this paper, the method employed follows closely that proposed by Fung but does not exhibit any of the limitations just discussed.
The paper is organized as follows. In the following section, the usual (Boltzmann) linear viscoelastic model is reviewed by the way of introduction to a new interpretation of Fung's QLV theory, presented in §3. In §4, results from the new model are contrasted with those in the literature for the case of uniaxial elongation of a bar: in §4a for an imposed stretch of an incompressible bar, in §4b for an imposed load, and in §4c for stretch of a compressible bar. A numerical solution to a Volterra integral equation is required for the results in §4b, and the discretization (time-stepping) procedure used is described in appendix A. Finally, concluding remarks are offered in §5.

Boltzmann's linear viscoelastic law
It is helpful to commence analysis by recapping the theory of linear viscoelasticity. Under the assumption of isotropy, infinitesimal elastic deformations can be described by the constitutive law σ(t) = 2μ( (t) − 1 3 tr( (t))I) + κ tr( (t))I, (2.1) where σ and are the second-order stress and strain tensors, respectively, and μ and κ represent the infinitesimal shear modulus and modulus of compression (or bulk modulus), respectively. To incorporate viscoelastic behaviour, the most natural extension of (2.1) is to assume that the stress remembers the past history of the rate of strain with some fading memory, and then apply the superposition principle (or Boltzmann's principle); hence, where μ(t), κ(t) are now time-dependent relaxation functions, and the lower limit of the integral must be taken from −∞ in order to correctly consider the initial deformation. Integrating (2.2) by parts, and assuming that the deformation commences at t = 0, yields where the denotes differentiation with respect to the argument of the function. Note that this expression incorporates any jump discontinuity when the motion starts. Moreover, the first term in (2.2) (or equivalently the first two terms in (2.3)) is trace free, or deviatoric, and accounts for shear deformations and loss, while the second term accounts for the hydrostatic part representing compressive deformations and loss. Now, in the elastic case (2.1), incompressibility may be considered as the dual limit of κ/μ → ∞ and tr → 0, which results in the stress-strain relationship i.e. the second term in (2.1) yields a finite non-zero limit, where p(t) may be considered as a Lagrange multiplier of incompressibility. For the viscoelastic counterpart, the limit of incompressibility can be considered in a similar fashion, noting first that μ(∞) and κ(∞) are the long time shear and bulk moduli (the equivalent of μ, κ in the elastic case), respectively. Thus, taking κ(∞) → ∞ (which implies κ(t) → ∞ for all t owing to the fading memory assumption) and tr → 0 in (2.2) (or in (2.3)), it is found that where the limit is again assumed to have a non-zero finite value, −p(t). Thus, in the limit of incompressibility, equations (2.2) and (2.3) become and respectively.

Quasi-linear viscoelasticity
When the strain is not infinitesimal, linear theory becomes inappropriate to describe deformations; hence, a nonlinear constitutive law has to be considered. As discussed in the Introduction, Fung's hypothesis (1.1) is examined here as a means of describing the motion of viscous nonlinearly elastic materials. In his renowned work [2], Fung introduced this quasi-linear constitutive model in order to capture the nonlinear stress-strain relationship of living tissues; however, it also has applicability to elastomeric materials. Before deriving the quasi-linear theory, it is useful to introduce some standard definitions and notations. The deformation gradient tensor F is defined by with x(s) denoting the position of a generic particle P at time s ∈ [0, t], and X its position at the initial reference time. Note that the start time of the deformation, and any imposed tractions, will be taken as t = 0. The quantity J = det F, expressing the local volume change, is a constant J = 1 when the deformation is isochoric. Furthermore, from the deformation gradient tensor F the left Cauchy-Green tensor B = FF T is obtained, together with its principal invariants which, alternatively, can be expressed in terms of the principal stretches through Now, Fung [2] makes the assumption that the QLV stress depends linearly on the (superposed) time history of a related nonlinear elastic response (a nonlinear instantaneous measure of strain). This allows, for example, for incorporation of a finite hyperelastic theory in the analysis. In index notation, Fung's theory (1.1) can be written as 1 Fung refers to G ijkl as a reduced relaxation function tensor. The 'crucial' simplification of Fung's theory is that this term is independent of the strain. Moreover, if the material is isotropic then G, a tensor of rank four, has just two independent components, being therefore consistent with linear theory. 2 Following the analysis of the previous section, for isotropic materials it is convenient to split the equivalent (instantaneous) Cauchy stress into two parts, one which accounts for microscopic isochoric deformations of the material and the other that measures purely compressive deformations. These two components can be expected to have different instantaneous elastic behaviours as well as distinct relaxation rates, associated for example with the unwinding/unravelling of polymeric fibres as opposed to their stretching. It is assumed that this decomposition can be achieved by taking the deviatoric and hydrostatic components of the equivalent elastic Cauchy stress: which can be expressed as The split of equation (3.5) into shear and dilatational components is the nonlinear (hyperelastic) analogue of linear stress-strain law (2.1). However, it cannot be generalized to viscoelasticity, as in (2.2), because it would not preserve objectivity. Instead, the second Piola-Kirchhoff stress tensor associated with (3.5) has to be introduced, which is defined by with It must be emphasized that the subscripts D and H do not refer to the deviatoric and hydrostatic parts of the second Piola-Kirchhoff stress, but correspond to the second Piola-Kirchhoff stress of the deviatoric and hydrostatic Cauchy stress components, respectively. Assuming a superposition principle as for the linear case, it is now possible to introduce an objective viscoelastic law, relating the second Piola-Kirchhoff stress to the past history of the nonlinear rate of strain measure. This is taken as where now D(t) and H(t) are two scalar (independent)-reduced relaxation functions (with D(0) = H(0) = 1 without loss of generality). The latter relaxation functions relate to the inherent viscous processes involved with instantaneous shear and compressional (volumetric) deformations, respectively. Clearly, if the material was anisotropic, then a more complex tensorial relaxation function would be required. Furthermore, pre-multiplying by J −1 F and post-multiplying by F T yield the Cauchy viscoelastic stress (3.11) and integrating by parts, following (2.7), gives As mentioned earlier, to allow for large (nonlinear) deformations, a hyperelastic theory can be employed assuming that the measure of the effective elastic stress Π e is derived from an elastic potential. Let us then specialize equations (3.5)-(3.12), considering the existence of a strain energy function (SEF) W, say, which in the isotropic case is dependent on the principal invariants of the deformation I 1 , I 2 , I 3 , or of the principal stretches λ 1 , λ 2 , λ 3 : The general form of elastic Cauchy stress may be written (e.g. [22,23]) as where β j = β j (I 1 , I 2 , I 3 ) are the so-called material (or elastic) response functions. In terms of the strain energy function, they are given by It is straightforward to calculate the trace of T e , .
From this, the second Piola-Kirchhoff counterparts, (3.9), are given by and where it is recalled that C = F T F is the right Cauchy-Green deformation tensor. The viscoelastic stress is obtained from (3.19) to (3.20) via (3.10) or for the Cauchy stress from (3.12).
Note that (as required) objectivity is now preserved for the viscoelastic model (3.11), but the first part is not deviatoric, and so the second term is not purely hydrostatic. In fact, in general both the compressive and shear components of the stress history contribute to the deviatoric and hydrostatic parts of T. However, the main point to note is that when there is no deformation, i.e. the principal stretches are unity, λ j = 1, j = 1, 2, 3 and B ≡ I, then I 1 = I 2 = 3, I 3 = 1, which reveals that the effective deviatoric elastic stress T e D vanishes. Similarly, T e H vanishes as the strain energy function has always to satisfy the additional relation (e.g. [22]) Thus, both terms Π e D and Π e H vanish as λ i → 1, and hence, as there is no applied stress until the initial time t = 0, justifies taking the integration range for the pair of integrals in (3.12) as 0 to t.
If it proves convenient to express the strain energy function in terms of the principal stretches W =W(λ 1 , λ 2 , λ 3 ), then for diagonal F the quasi-linear viscoelastic stress relation (3.12) may be rewritten as where nowW i refers to the derivative ∂W/∂λ i . The final consideration of this section is the constraint of incompressibility for all possible deformations, i.e. J ≡ 1. In this limit, the speed of propagation of compressive disturbances tends to infinity, and hence the relaxation time for viscous dilatational motions can be assumed to tend to zero. Therefore, the second term in (3.11) (or (3.12)), in an analogous fashion to that for linear elasticity (2.7), reduces to  (3.25) where now from the first equation in (3.19), (3.26) As described above, the present theory has been developed for isotropic materials. The study of anisotropic nonlinear viscoelastic materials rapidly becomes rather demanding, although strain measures for anisotropic QLV were discussed in [24]. It is expected that, despite the complexity, the present approach is extendable to certain classes of anisotropy.

Uniaxial elongation
One of the simplest, useful and hence most popular experiments to measure the properties of a material is to subject a specimen to a simple elongation test. There is a huge literature of results on such a deformation, and it is therefore useful to examine this homogeneous solution, comparing the present result with extant published work. To appreciate the general QLV model developed herein, and the approach needed to obtain a solution, three specific cases will be examined. First, in §4a a uniaxial stretch is imposed on the specimen, with the stress determined as a relaxing function of the history of the stretch. For comparison, both Yeoh and Mooney-Rivlin incompressible strain energy functions are considered. In §4b, a tensile load is imposed, and the stretch history (creep) determined from this. In this case, the integral equation must be solved numerically, which is the focus of discussion in appendix A. Finally, in §4c, the procedure in §4a is repeated for a compressible material, with the effect of compressibility on the viscoelastic behaviour highlighted.

(a) Simple extension
In the case of simple extension, assumed homogeneous (in space) and incompressible, the principal stretches may be specified as where (X 1 , X 2 , X 3 ) and (x 1 , x 2 , x 3 ) are the Cartesian coordinates in the undeformed and deformed state, respectively, and the stresses are having assumed that the lateral surfaces are free of stress. The deformation gradient and the constraint of incompressibility, J = 1, gives a relationship between the stretches λ 1 and λ 2 , in particular λ 2 = λ −1/2 1 . Setting λ 1 = λ, the deformation gradient tensor F and the left Cauchy-Green tensor B = FF T become and The principal invariants are therefore and hence taking the diagonal terms of equation (3.25) yields the principal Cauchy stresses respectively. Thus, by subtraction, the Lagrange multiplier p can be eliminated, and so it is possible to rewrite the stress-strain relationship (4.6) and (4.7) as where from (3.26) and Π e D22 = 2 Note that equations (4.9) and (4.10) depend on the specific choice of strain energy function W, and so it is useful to examine two specific examples. Let us start by assuming the instantaneous response is modelled by a two-term Yeoh strain energy function [25] which yields the stress-strain relationship where α is a positive constant and μ is the shear modulus from infinitesimal theory. Thus, It is straightforward to obtain from (4.8) the relation where The second example is the instantaneous response modelled by a Mooney-Rivlin strain energy function which yields the stress-strain relationship where γ is a constant in the range −1/2 ≤ γ ≤ 1/2. Note that (4.17) reduces to the neo-Hookean model when γ = 1/2. The partial derivatives of W are and so (4.8) becomes, after simplification,  found for the Yeoh model (4.14) (dotted) and for the Mooney-Rivlin material (4.19) (dashed) where the solid curve is the neo-Hookean limit α = 0 (or γ = 1/2). (c) T/μ is plotted from the predictions of (4.24) (dotted), (4.23) (dashed) and (4.14) (solid), respectively. (d) The dimensionless stress, T/μ, is plotted against stretch, λ, from the predictions of (4.24) (dotted), (4.23) (dashed) and (4.14) (solid).
in which (4.20) These two viscoelastic models can be compared by considering an 'experiment' where the stretch is imposed and the stress is measured. To specify matters, the stress relaxation function is chosen to be the classical one-term Prony series where μ, μ ∞ are the infinitesimal shear modulus and the long-time infinitesimal shear modulus, respectively, and τ is the relaxation time. These are set as μ ∞ μ = 0.5 and τ = 1.0s. A dynamic imposed stretch history (shown in figure 1a) is applied. The time variation is assumed slow so that inertial terms in the balance equations can be neglected. Figure 1b shows the resultant stress predictions T/μ for the Yeoh hyperelastic model (4.14) (dotted curves with α = 1, 2) and for the Mooney-Rivlin hyperelastic model (4.19) (dashed curves with γ = 1/6, −1/3). It is interesting to observe that the two hyperelastic models depart from the solid curve, obtained when the strain energy function is of the neo-Hookean type, i.e. α = 0, γ = 1/2, in different ways. The Yeoh material is found to harden as the parameter α increases, whereas the effect of decreasing γ from 1/2 leads to a softening of the Mooney-Rivlin material's behaviour. As mentioned in the Introduction, it is useful to contrast the results presented here with those discussed, for example, by Ciambella et al. [15]. In that article, the authors employed a form of QLV, which, with the incompressible Yeoh SEF, yields the stretch-stress law (see eqn (16) in [15]): Similar equations can be derived from the constitutive models in [16] or [17]. In Ciambella et al.'s article, their equation (4.23) was compared with that offered by the formulation used in the ABAQUS finite-element analysis (FEA) package (Version 6.7, see [26]), namely The three alternative viscoelastic expressions (4.24), (4.23) and the present result (4.14) may be conveniently compared by again imposing a stretch and determining the resultant dimensionless stress T/μ. Employing the same stretch history as above (figure 1a), the relaxation function in respectively. An alternative way of representing this deformation is via a stretch-stress diagram (figure 1d). It is clear that the present QLV approach, illustrated by the solid curve in both figures, gives a result that lies remarkably close to that found from the ABAQUS model, whereas Ciambella et al.'s solution is quite distinct. The latter model was developed by the authors in order to address the deficiency they highlighted regarding v. 6.7 of ABAQUS. However, their result also exhibits a substantial limitation: it predicts instantaneous zero stress whenever λ recovers back to unity after some (arbitrary) deformation, thus showing no 'memory' of the history of the stress in this situation. This is at variance with physically observed behaviour, and in particular that found for linear viscoelasticity.

(b) Simple tensile load
Rather than an imposed stretch, perhaps a more important experiment for measuring viscoelastic properties of materials is application of a simple tensile load: a slowly varying uniaxial stress is imposed on the body for which the resultant stretch is measured over time. For linear viscoelasticity, equation (2.2) can be inverted to yield a straightforward creep relation to determine the strain. However, for nonlinear viscoelasticity, the solution procedure is somewhat more difficult, as inversion is not possible. A typical QLV relation is given in (4.19), which may be considered in the general form Although the integrand is of separable type, the presence of the various f j (λ(t)) terms means that an inversion operator cannot be introduced. Instead, a numerical scheme must be employed that can evaluate the stretch as a function of time, subjected to a prescribed stress. A suitable numerical discretization procedure is offered in appendix A, which has error O(δt 4 ), where δt is the step size and so gives rapid convergence. An example is illustrated in figure 2, with the imposed stress history shown on the left graph. The one-term Prony series relaxation function (4.21), with constant values as given in equation (4.22), is again employed. The resultant stretch is given in figure 2b for the Yeoh strain energy function (4.14) (dotted curves with α = 1, 2), and for the Mooney-Rivlin strain energy function (dashed curves with α = 1, 2). As before, the solid curve is the neo-Hookean prediction (α = 0,  , from the Yeoh model predictions (4.14) (dotted) and that from Mooney-Rivlin predictions (4.19) (dashed). The solid curve is the neo-Hookean limit α = 0 (or γ = 1/2). γ = 1/2) and it can be seen that increasing α in the Yeoh model leads to material hardening, while decreasing γ leads to softening of the Mooney-Rivlin material.

(c) Simple extension of a compressible bar
The purpose of this section is to repeat that analysis in §4a for a compressible material. As before, a bar, with zero tractions on its lateral faces, is assumed to undergo simple homogeneous extension, caused by an imposed stress or stretch in the X 1 direction, say. Then, the principal stretches may be given by x 1 (t) = λ 1 (t)X 1 , x 2 (t) = λ 2 (t)X 2 and x 3 (t) = λ 2 (t)X 3 , (4.26) where symmetry between the X 2 and X 3 directions is clear, and so the deformation gradient tensor is F(t) = diag(λ 1 (t), λ 2 (t), λ 2 (t)) (4.27) and the left Cauchy-Green tensor B becomes The principal invariants are therefore I 1 (t) = λ 2 1 (t) + 2λ 2 2 (t), I 2 (t) = 2λ 2 1 (t)λ 2 2 (t) + λ 4 2 (t) and I 3 = λ 2 1 (t)λ 4 2 (t). (4.29) The above can be substituted into the governing viscoelastic equation (3.12), with effective elastic stresses (3.19) and (3.20), and simplified. The X 1 component of the stress yields, after simplification, (W 1 (s) + W 2 (s)λ 2 2 (s)) ds + (2W 2 (s) + 3W 3 (s)λ 2 1 (s))λ 2 2 (s) ds where as before W j is the derivative of W with respect to the indicated invariant (3.16). Clearly, the particular choice of the strain energy function affects the results and so, as before, a couple of examples are presented. Assuming the Horgan-Murphy strain energy function [27], used for modelling a material with a small compressibility, the SEF is where γ is an arbitrary constant, μ is the infinitesimal shear modulus and κ is the infinitesimal bulk modulus. Then which are substituted in (4.30) and (4.31). Setting γ = 1/2, for ease of presentation, the stresses become respectively. By contrast, for the Gent model (e.g. [28]), where J m is a constant limiting value for I 1 − 3, taking into account the finite extensibility of polymeric chains within the material. Hence, which can then substituted into (4.30) and (4.31) to yield the equivalent results for the stresses. These are omitted for brevity.

Concluding remarks
This paper has focused on reappraising Fung's method for QLV. It has been shown that some of the negative features commonly associated with the approach are merely a consequence of the way it has been applied elsewhere. The approach outlined herein was shown in §4 to yield 'sensible' results and offers a straightforward approach in solving a wide range of models. The present method exhibits similarities with Simo's approach to nonlinear viscoelasticity [7], although the latter assumed a scalar relaxation function acting on the stress. Therefore, it is noted that Simo's method would not reduce in the linear limit to the usual relation (2.3) if the latter shear and bulk moduli have different temporal behaviour. (Note that ABAQUS v. 6.12 employs the same relaxation time constants in both the shear and bulk parts of the field.) Herein, the relaxation function is a tensor, which for isotropic materials reduces to two distinct scalar relaxation functions, one acting on the compressive part and the other on the shear component of the stress. This therefore is consistent with that found for linear theory. It was shown that for an imposed deformation (e.g. stretch) it is simple matter to solve the QLV equation directly. For imposed stress, the stretch has to be deduced by solving (inverting) the integral equation, and this is achieved here using a discretized scheme accurate to an order of the cube of the time-step size. Muliana et al. [29] have recently offered a QLV model where the strain is expressed as a function of the stress, which may be viewed as a dual model to (1.1). However, it is not clear as yet whether their approach is as effective at modelling viscoelasticity as Fung's scheme.
The authors are currently applying the present approach to a range of problems in viscoelasticity. The numerical procedure described in appendix A can be employed in a straightforward manner for any deformation that is incompressible and homogeneous, for example simple shear. It is presently being extended so that it can used for compressible materials undergoing simple extension or shear. The QLV method is also adaptable to inhomogeneous problems, for example those studied in [30], large amplitude radial deformations in two and three dimensions, and simple torsion. Such models lead to partial differential nonlinear Volterra integral equations and extending the present numerical scheme to these should offer improvement over current methods. It is also anticipated that the present approach will be more suited to satisfying the equations of equilibrium than, say, the more heuristic QLV relation employed in ABAQUS (v. 6.12). Finally, the authors are using the model for several biomechanics problems, to derive a perturbation theory for the viscoelastic evolution of small deformations on the top of large ones and to derive effective properties of ligaments and other tissues.
where T(t), D(t), g(X), f j (X) and h j (X) are all known functions of their respective arguments. Clearly, all the equations offered in §4a lie within this class but for more general problems, with inhomogeneous deformations, the analysis described below must be amended.
Let us now discretize the system and evaluate the equation at each time t n = nδt, with n = 0, 1, 2, 3, . . . , n max , where δt is a small time step. On writing the stretch at each time step as a Taylor series expansion can be employed to yield where denotes differentiation with respect to the argument, and Note that all functions are assumed sufficiently smooth between adjacent time steps so that the relations above and below hold. From (A 4), it is clear that and (λ n − λ n−1 ) 3 = (λ n−1 ) 3 δt 3 + O(δt 4 ).

(A 18)
As there is no stretch prior to an applied stress, an initial zero stress, T(0) = 0, would mean that λ 0 = 1. The derivatives of λ 0 are derived from (A 14) to (A 16) (setting n = 1 in these equations) and hence λ 1 is found from (A 4). Then n is incremented, the values of λ , λ , λ determined, and these are again used to determine the stretch at the next time step. The process is repeated at each following step. Note that size of the sums in A n−1 , B n−1 , C n−1 grows with n, and so for large times, or using small time steps, evaluation of the system can become rather slow. However, this difficulty can be alleviated, with the present system, by obtaining a relationship between A n and its previous increment A n−1 , and similarly for B n , C n . The precise forms of these relationships are determined from the given f j (λ) and the relaxation function D(t). In practice, the solution procedure works well, with an error of O(δt 4 ), i.e. for a time step of δt = 0.01 then the error is O(10 −8 ).