Page 329 - Advances in Biomechanics and Tissue Regeneration

16.5 IMPLEMENTATION 327

numerical integrators have been implemented, including forward Euler, backward Euler, Midpoint, Adams-
Bashforth, or adaptive Runge-Kutta solvers, such as Dormand-Prince scheme [43] or Bogacki-Shampine scheme [44].

16.5.1.3.1 FORWARD EULER METHOD
Forward Euler is the simplest numerical integrator. Using forward Euler integration, Eq. (16.78) writes
x t +1 ¼ x t + f t Δt (16.79)

where x t ¼ x(t), x t+1 ¼ x(t + Δt), and f t ¼ f(x t , t). Eq. (16.79) is just an evaluation because it has an explicit nature. How-
ever, it is known that the forward Euler method can also be numerically unstable, especially for stiff equations, requir-
ing very small time steps for obtaining accurate results.

16.5.1.3.2 BACKWARD EULER COMBINED WITH THE BROYDEN METHOD
The backward Euler method is a numerical integrator that may work for greater time steps than forward Euler, due
to its implicit nature. However, because of this, at each time-step, a multidimensional nonlinear equation must be
solved. Eq. (16.78) discretized by means of the backward Euler method writes
x t +1 ¼ x t + f t +1 Δt (16.80)
where x t ¼ x(t), x t+1 ¼ x(t + Δt), and f t+1 ¼ f(x t+1 , t + Δt). Eq. (16.80) is nonlinear in general and has to be solved iter-
atively. Computation of the tangent operator ∂f is computationally expensive and hard to derive from the expression of
∂x
the function f. Instead of using the classical Newton method, based on the tangent operator, the multidimensional
generalization of the secant method is going to be used, the Broyden method [45]. At each iteration, a secant operator,
1
J n , is going to be computed and an update of variables is performed as x t +1,n +1 ¼ x t +1,n J r t +1,n , where we have
n
defined the residual at iteration n r t+1, n ¼x t+1, n x t, n f t+1, n Δt and the secant operator is defined such as satisfying
(16.81)
J Δx t +1,n ¼ Δf t +1,n
n
where Δx t+1, n ¼x t+1, n+1 x t+1, n and Δf t+1, n ¼f t+1, n+1 f t+1, n .
Of course, if the dimension of x and f are greater than 1, Eq. (16.81) is undetermined and further conditions shall be
supplied. A possibility is to use the current estimate of the secant operator J n 1 and improving upon it by taking the
solution to the secant equation that is a minimal modification to J n 1 in terms of the Frobenius norm, that is kJ n J n 1 k F
is minimal. Thus
Δf n J n 1 Δx n T
J ¼ J n 1 + 2 Δx n (16.82)
n
k Δx n k
Using the Sherman-Morrison formula [46], the inverse of the secant operator may be updated, using Eq. (16.82)as
follows
Δx n J 1 Δf n T 1
n 1
1
1
J ¼ J + Δx J (16.83)
n n 1 T 1 n n 1
Δx J Δf n
n n 1
Finally, at each iteration, the residual f t+1, n is evaluated, where
r t +1,n ¼ x t +1,n x t,n f t +1,n Δt (16.84)
The iteration step stops when kr t+1, n k < TOL. Note that with this approach, at each iteration, we have an evaluation of
function f n , that is, a construction of operators M, K, and f.

16.5.2 1D Finite Element Implementation

16.5.2.1 Unidimensional Equations
Differential equations (16.1), (16.4) with boundary conditions (16.2), (16.3), (16.5), (16.6) may be expressed easily for
unidimensional problems.
Let us define
¼ C i , i ¼ 1,…,n
u i
(16.85)
u n + i ¼ S i , i ¼ 1,…,m
T
and u ¼ (u 0 , …, u n , u n+1 , …, u n+m ) .

II. MECHANOBIOLOGY AND TISSUE REGENERATION