The Java Applets

The following are links to the two solution graphs.

• Applet 1 - Linear and Non-linear applet
• Applet 2 - Non-linear animated applet with profile

Contents

• Introduction
• The Non-linear Model
• Linearisation
• The Numerical Investigation
• Explanation of Programme
• Observations and Results
• References
• Bibliography

Introduction

This project has been based on the paper Anguilliform body dynamics: modelling the interaction between muscle activation and body curvature by Graham Bowtell & Thelma L. Williams (1991). (Here after referred to as the paper.)

The paper explores the model of a lamprey as a simple two-dimensional rod and pivot model for it's mechanical structure, each pivot being controlled by a muscle segment attached via perpendicular extensions to the two rods. The muscle segments contain terms to model the damping and stiffness, and also produce a time dependent forcing term, these combine to produce the thrust.

The equations of motion are derived using a standard Langrangian approach, with the adapted co-ordinates of the angles each rod makes with the horizontal and the co-ordinates of the head position at the end of the first rod.

The system examined in the paper is concerned with the general case of N rods. This project examines the case where N=3.

The non-linear system is explored, unlike the paper, and linearised, assuming low curvature dynamics. The two models are compared.

Both the linear and non-linear models are solved numerically and programmed in Java. Applets are produced to allow parameter changes to the model and display the linear and non-linear solutions on a web page.

The non-linear model is also animated to infinite time. And the corresponding profile of the lamprey is calculated and displayed.

The possible parameters for the model will be discussed. The results will be interpreted, and their relevance to a real life animal examined.

The Non-linear Model

The Mechanics

Figure 1. Representation of the body. 3 light pivoted rods of length l with a mass m at each pivot and at both ends which are controlled by the muscle segments (M_s .) The fixed ends of each M_s are a distance w from the (s)th and (s+1)th rod. The (s)th rod makes an angle theta_s with the positive x-axis.

Figure 2. Each muscle M_s consists of a spring (stiffness mu), a dashpot (damping coefficient gamma) and a power source producing contractile tension F_s (t), taken as positive in the direction indicated by the arrows. Total thrust is given by the sum of the contributions of the three parallel components, and is positive in the direction shown by the arrows along AB.

From figures 1 and 2, it is clear that:

Figure 3.

The total thrust tending to increase the length AB is given by

where x_s / l is the extension per unit length which, from (1) is given by

And therefore:

Equations of Motion

To derive the equations of motion we use a standard Lagrangian approach (Woodhouse 1987) for a non-conservative system, using a kinetic energy term and generalized forces. We use the adapted coordinates, i.e. the minimum set of co-ordinates needed to specify the position of the system, of theta₁ , theta₂ , theta₃ , representing the angles between the 3 rods and the positive x-axis, and the position of the head defined as (x₀ ,y₀ ), see figure 1. The positions of the masses at each of the pivots and the tail are given, for segments of equal length l, by:

and the velocity of the kth mass is, therefore, given by the derivative:

So, the total kinetic energy is the energy at the head and the other three masses and is given by:

where

and

Thus after substituting, and applying the normal trigonometric identities, we have:

The equations of motion with respect to the angular coordinates are given by:

Where G-theta_s is the generalized force corresponding to the generalized coordinate theta_s .

Hence, we have:

Since no work is done by the system for small variations of the head position, the corresponding generalized forces are zero, hence we have:

and:

Substituting in (5), we have:

Again using the normal identities, this reduces to:

Similarly, we have:

Generalized Forces

From the paper we have:

And we choose the forcing term:

And introduce a coefficient to the forcing A, to control the amplitude, giving

Linearisation

If we assume low curvature dynamics, we can linearise the system by using the approximations:

then we have, from (3)

And substituting in (2):

and so from (11):

so we have:

And from (8,9,10):

Then solving for theta double dot

The Numerical Investigation

Consider the scheme

by forward difference approximations, where h is the step length.

The Non-Linear Model

Re-arranging and putting into matrix form, we have:

This can be solved by Gaussian Elimination, with from (11):

and

from equation (4).

The Linear Model

From equations (15,16,17):

Obtaining the Solution

So theta_i can be found by the backward difference

Explanation of Programme

The numerical schemes for both the linear and non-linear models are programmed in Java (as an applet.) They allow the input of the parameters mu, gamma, kappa and omega, as well as the coefficient to the forcing term, A. The step length and the number of steps can also be varied, subject to a maximum number of steps of 1200.

The output is in the form of three graphs, representing theta₁ ,theta₂ and theta₃ plotted against time steps. The solution to the non-linear scheme is shown coloured and the linear in white.

The vertical axis can be re-scaled to suit a particular parameter set, and the plots redrawn (but not recalculated.)

There is a second applet with the same parameter input, but solving only the non-linear model, which runs to infinite time. The vertical axis value must be selected before running it, but any paramter can be changed and the graphs `restarted'.

In addition, the animated profile of the lamprey is also shown synchronized to the animated graphs. It assumes the lamprey to be initially at rest lying along the positive x-axis with the head point (x₀ ,y₀ ) at the origin. With these initial conditions, the evolution is given by (as in the paper with N=3):

Observations and Results

Choice of Parameters

Both the total mass and length are taken to be 1. i.e. m=0.25 and l=0.333. The width of a lamprey is approximately 4\% of it's body length, so w=0.02.

After some experimentation the stiffness, mu, was set to 1000, and the damping, gamma, to 100. This project does not endeavor to explore these parameters.

Omega and kappa are kept equal so we consider a constant speed of the traveling wave. Rearranging (12):

and hence speed omega / kappa. We explore omega = kappa = pi and omega = kappa = 2pi.

The step length, h, for the numerical investigation is set to a default of 0.02, values greater than this introduce rapid rounding errors. And both models fail for values greater than 0.04.

Comparison of Linear and Non-linear Models

Firstly considering omega = kappa = pi. For a value of the forcing coefficient, A=3, the agreement is fair between the two models, with theta_1,3 only differing by approximately 0.008. There was excellent agreement for theta₂. If A=5 is used the agreement is poorer with a difference of approximately 0.03.

A step length of 0.02 was used over 600 time steps. Halving the step length and doubling the number of steps yielded the same result. Decreasing the step length further didn't improve the agreement either.

In the above, both models stayed constant, which is not true for the non-linear model if A is increased. The non-linear model becomes more negative as the time increases and the graph `falls away' from the linear model, giving a large disagreement. The non-linear model fails for A less than or equal to 7.

Secondly considering omega = kappa = 2pi. Here, F₁ = F₂ , equations (13) and (14), as the period of sin is 2pi, and it follows that G_s = G_s-1 so the forcing terms cancel in equation (16) for the linear model, and we expect similar behaviour in the non-linear model, so we expect theta₂ to be zero.

Here A=3 gave good agreement between the models initially, but the non-linear model falls off as it did for larger forcing above. For A=5, the two solutions parted much quicker, and theta₂ strays from the axis after around 100 steps. The non-linear model failed for A greater than 5.5.

We put the failure of the numerical method for the non-linear model down to rounding errors.

Stability of the Linear Model

Valid Angles

The methods used to obtain the three angles does not take into consideration the physical constraints of the model. From figure 4, it can be seen that if phi_s = theta_s+1 - theta_s becomes too negative, then the pivot would cross the line AB, which is not physically possible. Thus we must restrict the value of theta_s+1 - theta_s to a minimum.

Figure 4. The configuration of a pivot with the smallest value of phi_s = theta_s+1 - theta_s . The mass at the pivot lies on the line AB.

With reference to figure 4:

Also

The magnitude of the theta_s for the the values of A in 6.2, are very low. For the condition above to be maintained A less than 16 is needed for omega = kappa = pi and A less than 7 for omega = kappa = 2pi. Although, this is only true for the linear model, where the numerical method is stable under these forcings.

Relation to Real Life

The real lamprey, though, has around 100 segments, therefore to achieve the profiles given, consecutive angles do need to breech the condition we have imposed. Bowtell and Williams used much larger numbers of rods which meant they could get a much closer approximation to the real profile.

The Animated Profile

Conclusions

The linear model, however, is capable of producing stable and sensible results, that is, achieving the maximum allowable angles.

There is clearly a problem with our very simplified case of N=3, which proves that a larger value of $N$ is needed to be able to model the lamprey more satisfactorily.

References

Bowtell, G & Willimas, T.L 1991. Anguilliform body dynamics: modelling the interaction between muscle activation and body curvature. Phil. Trans. R. Soc. Lond B, 385-390.

Woodhouse, N.M.J. 1987. Introduction to Analytical Dynamics. Oxford University Press, 31-48.

Bibliography

Morgan, M. 1998. Using Java 1.2. Que.

Bell, D. & Parr, M. 1999. Java for Students. Prentice Hall.