To simulate the dynamics of the fish model, the differential equations of motion () must be integrated over time. The system is intrinsically stiff and can become unstable in a variety of situations that may occur in the simulated dynamic world (e.g. executing a left turn to avoid a sudden dynamic obstacle before completing a right turn). Therefore, the stability of the numerical solution is an important consideration. We employ a simple, numerically stable, semi-implicit Euler method.

The system of equations of motion Eq. () may be written in the standard matrix form for nonlinear, damped elasto-dynamic systems [Bathe and Wilson1976] as follows:

where , and are the mass, damping and
stiffness matrices, respectively. Matrix is a constant
diagonal matrix with on the main diagonal, while and
are time-dependent nondiagonal matrices. These three
matrices are of dimension where *n* is the number of
nodes in the biomechanical model. Matrices , ,
consist of nodal position, velocity and acceleration
vectors, respectively, and is the matrix of external force
vectors. These four matrices are of dimension , such that
each row contains the three components of each vector.

Because of the special properties of the viscoelastic units in our
dynamic fish system, the damping term
and the stiffness term in the above equation
are most conveniently collected into the form of an *effective
stiffness matrix* :

and hence Eq. () can be re-written as follows:

where is time-dependent since the viscoelastic forces depend nonlinearly on the nodal position and velocity variables and .

Our method computes the left-hand side of Eq. () implicitly and the right-hand side explicitly. In particular, we discretize continuous time into time steps, and solve for in

where *t* represents simulation time and is the simulation
time step. We do not solve for the position vectors directly because, first, the velocity vectors is
readily available since it has to be stored at each time step for
computing (more details can be found in the next section)
and, second, using saves us from having to store at each time step. We use the following finite
difference approximations for and :

and

Substituting the above equations into Eq. () and collecting all known vectors on the right-hand side, we obtain a linear system for :

where

is the effective system matrix and

is the effective load matrix. As we will explain in the next section, the structure the system matrix is such that an efficient solution to Eq. () can be obtained by assembling and factorizing it as follows:

where is diagonal and is lower triangular. The new nodal velocities are then computed by first solving the lower triangular system

for by forward-substitution and finally solving the upper triangular system

for by back-substitution. Hence, the numerical simulation of the fish biomechanical model requires the assembly and factorization solution of a sequence of linear systems of the form of Eq. for .

Xiaoyuan Tu | January 1996 |