The following are links to the two solution graphs.
This project aims to explore a model for the body dynamics of a swimming lamprey. It has been based on Anguilliform body dynamics: modelling the interaction between muscle activation and body curvature by Graham Bowtell & Thelma L. Williams (1991).
The model is a simple two-dimensional rod and pivot system. The model is linearised and both the linear and non-linear models are solved numerically. These are programmed in an interactive manner, and presented on a web page.
Various parameters for the models are explored. The stability of the schemes and the validity of the results are discussed.
The numerical method for the linear model is shown to be stable, and capable of yielding reasonable results. The numerical method for the non-linear model proves not to be able to produce the same quality of results, and it's inaccuracies make it inconclusive.
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.
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 (Ms .) The fixed ends of each Ms are a distance w from the (s)th and (s+1)th rod. The (s)th rod makes an angle thetas with the positive x-axis.The muscle unit is modeled by a linear spring providing the stiffness element, a linear velocity-dependent dashpot providing the damping and a power unit capable of producing time dependent tension and thrust, in line with the understanding of the real animal (figure 2).
Figure 2. Each muscle Ms consists of a spring (stiffness mu), a dashpot (damping coefficient gamma) and a power source producing contractile tension Fs (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.With reference to figure 3, then AB can be calculated thus:
The total thrust tending to increase the length AB is given by
where xs / l is the extension per unit length which, from (1) is given by
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 theta1 , theta2 , theta3 , representing the angles between the 3 rods and the positive x-axis, and the position of the head defined as (x0 ,y0 ), 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:
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-thetas is the generalized force corresponding to the generalized coordinate thetas .
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:
Substituting in (5), we have:
Again using the normal identities, this reduces to:
Similarly, we have:
From the paper we have:
And we choose the forcing term:
And introduce a coefficient to the forcing A, to control the amplitude, giving
Back to Contents Applet 1 Applet 2
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
Back to Contents Applet 1 Applet 2
Consider the scheme
by forward difference approximations, where h is the step length.
Re-arranging and putting into matrix form, we have:
This can be solved by Gaussian Elimination, with from (11):
from equation (4).
From equations (15,16,17):
So thetai can be found by the backward difference
Back to Contents Applet 1 Applet 2
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 theta1 ,theta2 and theta3 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 (x0 ,y0 ) at the origin. With these initial conditions, the evolution is given by (as in the paper with N=3):
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.
Firstly considering omega = kappa = pi. For a value of the forcing coefficient, A=3, the agreement is fair between the two models, with theta1,3 only differing by approximately 0.008. There was excellent agreement for theta2. 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, F1 = F2 , equations (13) and (14), as the period of sin is 2pi, and it follows that Gs = Gs-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 theta2 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 theta2 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.
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 phis = thetas+1 - thetas becomes too negative, then the pivot would cross the line AB, which is not physically possible. Thus we must restrict the value of thetas+1 - thetas to a minimum.
Figure 4. The configuration of a pivot with the smallest value of phis = thetas+1 - thetas . The mass at the pivot lies on the line AB.
With reference to figure 4:
The magnitude of the thetas 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.
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 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.
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.
Morgan, M. 1998. Using Java 1.2. Que.
Bell, D. & Parr, M. 1999. Java for Students. Prentice Hall.