This paper was awarded in the V International Competition (1996/97) "First Step to Nobel Prize in Physics" and published in the competition proceedings (Acta Phys. Pol. A 93 Supplement, S-55 (1998)). The paper is reproduced here due to kind agreement of the Editorial Board of "Acta Physica Polonica A".

Some Generalisations of Brachistochrone Problem

A.S. PARNOVSKY

Lyceum No. 145 for Natural Sciences, Sh. Rustaveli 46, Kiev-23, 252023, Ukraine
Abstract
300 years ago Johann Bernoulli solved the problem of brachistochrone (the problem of finding the fastest travel curve's form) using the optical Fermat concept. In the same way we solved some generalisations of this problem. We obtained the fastest travel curve's form in a gravitational field for a point-like mass. When moving in a uniform gravitational field we considered the influence of the dry friction or drag.


PACS numbers: 03.20.+i, 46.10.+z



1. Introduction


Last year there was the 300th anniversary of solving the problem of curve of the fastest travel. It was the first step in the development of variational calculus. In spring 1696 Johann Bernoulli proposed a brachistochrone problem (from the Greek $\beta\rho\alpha\chi\upsilon\sigma$ -- short, $\chi\rho o\nu
o\sigma$ -- time) in the famous scientific magazine Acta Eruditorum. A material particle moves without friction along a curve. This curve connects point A with point B (point A is placed above point B). No forces affect it, except the gravitational attraction. The time of travel from A to B must be the smallest. This brings up the question: what is the form of this curve?

Bernoulli announced that he had found an excellent solution of this problem, but he did not want to publish it in order to prompt prominent mathematicians to apply their skill to the mathematical problems of a new type. In particular, he challenged to a competition his elder brother Jacob. At that time they were in hostile relations, and Johann openly called Jacob an ignorant. Jacob Bernoulli published his solution of the problem in the same year. It is intriguing that the geometrical optics' laws were applied to a mathematical problem dealing with mechanical motion.

The main ideas of this solution are given in the book [1]. On this basis a differential equation of a brachistochrone is built and solved in the next section of this article. It will be shown that the fastest travel curve is an arc of a cycloid. Some generalisations of the problem are considered in the next sections. There is a friction affecting the particle in them. There can be the dry friction in accordance with the Coulomb-Amontons law, $F=\mu N$, as well as the fluid friction proportional to velocity. We also considered the fastest travel trajectories of a point-like or spherical mass moving in a gravitational field.


2. A solution of the classical brachistochrone problem


During the motion of the particle along brachistochrone its velocity v and angle $\alpha$ between the direction of this velocity and the vertical changes. In order to find the connection between these quantities, Jacob Bernoulli used the analogy with the discovered in the XVII century Fermat concept. It states that the light moves in the medium with varying refractivity in such a way that its travel time is minimal. In the case when the refraction index n depends only on y coordinate the Snellius (Snel van Royen) law states that the quantity $n(y)\sin\alpha$ is constant. Here $\alpha$ is the angle between the y axis and light path. The light velocity in medium is v=C/n and the relationship $v\propto\sin\alpha$ is fulfilled for ray's trajectory. Therefore, the equation

\begin{displaymath}v=C\sin\alpha,\qquad C\rm =const\eqno(1)\end{displaymath}

must be valid also for the brachistochrone. Let us prove that this condition suffices to find the form of the brachistochrone. We introduce the rectangular Cartesian coordinates with y axis directed vertically downwards as shown in Fig. 1. The coordinate origin is placed at point A. Point B has the

\includegraphics[width=8cm]{par1}

$\textstyle \parbox{0.9\textwidth}{{\small Fig.~1.}~The definition of the $(x,y)
$\ coordinate system and angles $\alpha$\ and $\varphi$.}$

coordinates (x0>0, y0>0). The curve inclination to horizon $\varphi$ is connected to the derivative y' (here and further we will denote the derivative with respect to x with a prime) with an equation $y'=\tan\varphi$. Moreover, $\varphi=\pi/2-\alpha $, whence $y'\rm =cotan\alpha$. Thus, we can write Eq. (1) in the following form:

\begin{displaymath}v={C\over\sqrt{1+(y')^2}}.\eqno(2)\end{displaymath}

From the condition of the conservation of energy

\begin{displaymath}v=\sqrt{2gy}.\eqno(3)\end{displaymath}

Substituting (3) in (2), we obtain the sought-for differential equation

\begin{displaymath}y={2C_1\over1+(y'
)^2},\qquad C_1={C^2\over4g}\rm =const,\eqno(4)\end{displaymath}

whence

\begin{displaymath}
y'=\sqrt{\frac{2C_1} y-1}\eqno(5)\end{displaymath}

and

\begin{displaymath}x=\int\frac{{\rm d}y}{\sqrt{2C_1/y-1}}.\eqno(6)\end{displaymath}

After substitution

\begin{displaymath}y=C_1(1-\cos\lambda)\eqno(7)\end{displaymath}

we obtain

\begin{displaymath}x=C_1\int\frac{\sin
\lambda\rm d\lambda}{\sqrt{\frac2{1-\cos...
...da)
{\rm d}\lambda=C_1(\lambda-\sin\lambda)\rm +const.\eqno(8)\end{displaymath}

From the initial conditions C2=0 and Eqs. (6), (8) define the brachistochrone equation in a parametrical representation. One can see that it is a cycloid with radius C1.

Value $\lambda=0$ corresponds to the point A and $\lambda=\lambda_0$ to the point B. We can obtain $\lambda_0$ from the conditions

\begin{displaymath}x_0=C_1(\lambda_0-
\sin\lambda_0),\qquad y_0=C_1(1-\cos\lambda_0),\eqno(9)\end{displaymath}

or, more precisely, from

\begin{displaymath}x_0/y_0=(\lambda_0-\sin\lambda_0)/(1-\cos\lambda_0).\eqno(10)\end{displaymath}

\includegraphics[width=8cm]{par2}

$\textstyle \parbox{0.9\textwidth}{{\small Fig.~2.}~The dependence of $\lambda_0...
...(\lambda_0,~x_0$, and $y_0$are the $\lambda,~x,~y$\ values for the point $B$).}$

Figure 2 shows the dependence of $\lambda_0$ on x0/y0 which was calculated using a computer. Parameter C1 can be obtained from (9) using the calculated value of $\lambda_0$. Using (3), (7) and (8), we obtain

\begin{displaymath}T=
\int_0^{y_0}\frac{{\rm d}y}{v_y}=\int_0^{y_0}\frac{{\rm d...
...phi}=
\int_0^{y_0}\frac{{\rm d}y\sqrt{1+(y')^2}}{\sqrt{2gy}y'}\end{displaymath}


\begin{displaymath}=\int_0^{\lambda_
0}\frac{C_1\sin\lambda{\rm d}\lambda\sqrt{...
...{\lambda_0}{\rm d}\lambda=\sqrt{\frac{C_1}g}\lambda_0.\eqno(11)\end{displaymath}


3. The generalisations of the brachistochrone problem


Below we will consider some generalisations of the brachistochrone problem. When generalising, we tried to adhere to the following rules:

As a result we choose the following problems:

The giant brachistochrone. If curve's size becomes comparable to Earth's one, one has to take into account the nonuniformity of the gravitational field. This problem will be considered in Sec. 4.

Brachistochrone with friction. Friction was added when moving. A dry friction with the friction index $\mu$ is considered in Sec. 5. The affection of medium's drag, proportional to the velocity (low Reynolds numbers) is considered in Sec. 6.


4. The giant brachistochrone


In this problem one has to find the form of the fastest travel curve of a point-like or spherical mass M in a gravitational field. It would be reasonable to use the polar coordinates with an origin placed at the centre of the attracting body. We denote the radial coordinate of the point A by r1; r2 is the same for the point B. The angle coordinate $\theta$ will be reckoned from the direction of point A. The condition of the conservation of energy will take the form

\begin{displaymath}v=\sqrt{2\gamma M(r^{-1}-r_1^{-1})
},\eqno(12)\end{displaymath}

instead of (3). In Eq. (12) we denoted the gravitational constant by $\gamma$. Let us use the Fermat concept and the analogy with the light propagation. Inasmuch the effective refraction index n depends on radius r, ray's trajectory can be given by the Bouguer law

\begin{displaymath}n(r)r\sin \alpha\rm =const,\eqno(13)\end{displaymath}

where $\alpha$ is the angle between the radius-vector and the direction of the ray's propagation and

\begin{displaymath}\sin\alpha=(1+ r,_\theta^2/r^2)^{-1/ 2},\eqno(14)\end{displaymath}

where with the $\theta$ index after the comma we denoted the derivative with respect to $\theta$.

Combining Eqs. (12)-(14) and setting n equal to C/v we obtain the sought-for differential equation

\begin{displaymath}r,_\theta=r\sqrt{\frac{C_2r^2}{r^{-1}-r_1^{-1}}-1},\qquad C_2\rm =const.\eqno(15)\end{displaymath}

Its solution with given initial conditions is

\begin{displaymath}\theta=\int_r^{r_1}\frac{{\rm d}r}r\sqrt{\frac{1-r/r_1}{C_2r
^3-1+r/r_1}}.\eqno(16)\end{displaymath}

To simplify this expression, we introduce the new variable z and the new constant z0, which is the value of the variable z when the denominator of the root expression vanishes. The previous constant C2 becomes equal to r1-3(z0-3-z0-2). Thus,

\begin{displaymath}\theta
=\int_z^1\frac{{\rm d}r}r\sqrt{\frac{1-z}{(z-z_0)[(1-z_0)z^2+z_0(1-z_0)z+z_0
^2]}}.\eqno(17)\end{displaymath}

We were able to solve the integral only numerically. (See the details of the numeric calculations in Sec. 7.)

\includegraphics[width=5cm]{par3}

$\textstyle \parbox{0.8\textwidth}{{\small Fig.~3.}~The set of particle's trajec...
...es starting from point $A$\ for the \lq\lq giant
brachistochrone'' generalisation.}$

\includegraphics[width=5cm]{par4}

$\textstyle \parbox{0.8\textwidth}{{\small Fig.~4.}~The body's movement trajectories for the \lq\lq brachistochrone with friction''
generalisation.}$

A set of trajectories which start from the same point, corresponding to the values of parameters equal to $0.1,0.2,\dots,0.9$ is shown in Fig. 3. One can see that the particle moves nearly radially at the start, deflecting increasingly during the motion. From Eq. (13) one can see that when $\alpha\neq0$ the particle cannot reach the centre. It must reach the minimal radius r1z0 and begin receding from the attracting body. If the minimal distance accords to the azimuthal angle $\theta_0$, the trajectory will be symmetrical about the $\theta=\theta_0$ ray. The trajectory ends on r1 distance from the centre of attraction in the point in which the azimuthal angle is equal to $2\theta_0$ and the particle's velocity vanishes. Let us consider the relation $\theta_0(z_0)$. Introducing new variable U=(z0/z)2/3, we obtain from (17)

\begin{displaymath}\theta_0=z
_0^{3/2}\int_{z_0}^1\frac{{\rm d}z}z\sqrt{\frac{1-z}{(z-z_0)[(1-z_0)z^2+z_0(
1-z_0)z+z_0^2]}}
\end{displaymath}


\begin{displaymath}\qquad=\frac23\int_{z_0^{2/3}}^1{\rm d}U\sqrt{\frac{1-z_
0U^{-2/3}}{1-z_0-U^2+z_0U^{4/3}}}.\eqno(18)\end{displaymath}

In Fig. 4 the plot of these relations is shown. It was obtained using the numerical integration. One can see that it differs from the straight line only slightly. But the most interesting property is the fact that when $z_0
\rightarrow0,~\theta_0\rightarrow\pi/3$. Actually, from (18) one can see that

\begin{displaymath}\theta_0\mathop{\longrightarrow}_{z_0\rightarrow0}\frac23\int...
...{\sqrt{1-U^2}}=\frac23\arcsin U\vert _0^1=\frac{\pi}3.\eqno(19)\end{displaymath}

Thus, the body moving along the trajectories depicted in Fig. 3, cannot change its azimuth more than by $120^\circ$. Therefore they are useless when point's B azimuth is from $120^\circ$ to $240^\circ$. In this case the fastest travel, apparently, is the falling to the centre of the attracting body and the consequent radial movement from the centre to point B. Let us find the time of falling to the centre. Using (12) and substituting U=[r/(r1-r)]1/2, we obtain

\begin{displaymath}T=\int_0^{r_1}\frac{{\rm d}r}v=(2\gamma M)^{-1/2}\int_0^{r_1}
\frac{{\rm d}r}{\sqrt{r^{-1}-r_1^{-1}}}\end{displaymath}


\begin{displaymath}\qquad=(\gamma M/2)^{-1/2}r_1^{3
/2}\int_0^\infty\frac{U^2{\...
...{(1+U^2)^2}=\frac{\pi r_1^{3/2}}{2\sqrt{2
\gamma M}}.\eqno(20)\end{displaymath}

However, in the centre of attraction the body's velocity with infinite norm undergoes an abrupt change of direction when $\theta_2\neq\pi$. It is physically prohibited. Apparently, among the physically permissible trajectories there is no curve with the minimal travel time. However, there are trajectories with travel time differing from the AOB travel time as small as possible.


5. Brachistochrone with friction


Let us return to the brachistochrone problem in an uniform gravitational field. We set the friction coefficient of the particle with the surface equal to $\mu$. Equations (1) and (2) will be still valid, but Eq. (3) will be violated because of the energy losses for the friction. We will find the equation, generalising it. The supporting force N consists of the weights component $mg\cos\varphi$ and the centrifugal force mv2/R=mv2 k, where R and k are the radius of curvature and the curvature of the trajectory

\begin{displaymath}
k=\frac{\vert y''\vert}{\sqrt{[1+(y')^2]^3}}.\eqno(21)
\end{displaymath}

Let us suppose that the trajectory is curved in the same way as the cycloid does (y''<0), so the centrifugal force is oriented in the same direction as the normal component of the weight is. Otherwise, the body can break away from the trajectory. For the infinitesimal shift dS along the trajectory the work of friction forces is equal to

\begin{displaymath}
{\rm d}A=-\mu{N\rm d}S=-\mu{m\rm d}S
(g\cos\varphi+kv^2)=-\mu{m\rm d}x(g+kv^2/\cos\varphi).\eqno(22)\end{displaymath}

The variation of kinetic energy is equal to

\begin{displaymath}{\rm d}(mv^2/2)=mg{\rm d}y{\rm +d}A=
m[g({\rm d}y-\mu{\rm d}x)-\mu{kv\rm ^2d}x/\cos\varphi].\eqno(23)\end{displaymath}

$\mu kv^2 /\cos\varphi>0$, therefore, the body cannot get to the area $y<\mu x$, bounded with the friction angle $\theta\rm =arctan\mu$. Substituting in (23) with Eqs. (2), (21) and $\cos^{-1}\varphi=[1+(y')^2]^{1/2}$, we obtain

\begin{displaymath}
\frac{C^2}2{\rm d}\{[1+(y')^2]^{-1}\}{\rm =d}x\left\{g(y'-\mu)+\frac{C^2\mu y
''}{[1+(y')^2]^2}\right\}.\eqno(24)\end{displaymath}

From this we obtain the differential equation of the sought-for trajectory

\begin{displaymath}C^2y''(\mu+y')=g(\mu-y')[1+(y')^2]^2. \eqno(25)\end{displaymath}

If we do the rearrangement of homothety with an origin at point A

\begin{displaymath}y^*=Ky,\qquad x^*=Kx,\eqno(26)\end{displaymath}

the first derivative y' will not change and the second derivative y'' will decrease K times. Setting K=C 2/g, we will rearrange (25) to the form

\begin{displaymath}
y''(\mu+y')=(\mu-y')[1+(y')^2]^2,
\eqno(27)\end{displaymath}

not containing parameter C. Thus, the sought-for trajectories can be mapped to each other with the scaling (26). Therefore it will be enough to find one of them, for example (27). Its first integral is equal to

\begin{displaymath}x=\int\frac{{\rm d}(y')(y'+\mu)}{(\mu-y')[1+(y')^2]^2}=\frac{y'(\mu^2-1)-2
\mu}{2[1+(y')^2](1+\mu^2)}
\end{displaymath}


\begin{displaymath}
\qquad+\frac{\mu^4+4\mu^2-1}{2(1+\mu^2)^2}\left
({\rm arcta...
...)-\frac\mu{(1+\mu^2)^2}\ln{(y'-\mu)^2\over1+(y
')^2}.\eqno(28)\end{displaymath}

The integration constant is found from the condition $y'(0)
=\infty$. The details of calculation are given in Sec. 7.

It is easy to test that for $\mu=0$ this expression reduces to the equation of the cycloid with radius 1/4. When $\mu\neq0$, the shape of the trajectories dramatically changes. From (25) one can see that at y''<0 expressions $y'+\mu$ and $y'-\mu$ must have the same signs, therefore $y'>\mu
$. Contrary to a cycloid, brachistochrone with friction cannot hook up and even becomes horizontal. In Eq. (28) x can tend to infinity only if one of terms in the right part diverges. It is possible only if $y'\rightarrow\mu$. Thus, if brachistochrone with friction can have infinite length, it tends asymptotically to the friction angle as

\begin{displaymath}y'-\mu\propto\exp[-x(1+\mu^2)^2/2
\mu].\eqno(29)\end{displaymath}

Expanding the right part of (28) into a series in reciprocal powers of y' at $y'\rightarrow\infty$, we obtain

\begin{displaymath}
x=(1/3)(y')^{-3}+\Lambda (y')^{-4}+
\dots,\end{displaymath}


\begin{displaymath}
\Lambda =\mu(1-6\mu^2+16\mu^3+9\mu^4-16\mu^5)/2(1+\mu^2)^2,
\eqno(30)\end{displaymath}

so for small x

\begin{displaymath}
y\approx\sqrt[3]9/2x^{2/3}+\Lambda x.\eqno(31)\end{displaymath}

The trajectories of the body at $\mu=0,~0.2,~0.5,~1$, obtained using the numeric integration of Eq. (28) are shown in Fig. 4. The details of calculations are given in Sec. 7.


6. Brachistochrone with drag


If the drag force is proportional to velocity

\begin{displaymath}
\mbox{\boldmath $F$}=-k\mbox{\boldmath $v$},\qquad k\rm =const,\eqno(32)\end{displaymath}

then we have

\begin{displaymath}
{\rm d}\left(\frac{mv^2}2\right)=mg{\rm d}y-kv{
\rm d}S=mg{\rm d}y-kv{\rm d}x/\cos\varphi\eqno(33)\end{displaymath}

instead of (23). Taking into account Eq. (2) we obtain

\begin{displaymath}C^2y'y''=g[1+(y')^2]^2(y'-\beta),\qquad\beta=kC/gm\rm =const.\eqno(34)\end{displaymath}

From this

\begin{displaymath}
-\frac{gx}{C^2}=\int\frac{y'{\rm d}(
y')}{[1+(y')^2](y'-\beta)}=\frac{y'+\beta}{2[1+(y')^2](1+\beta^2)}
\end{displaymath}


\begin{displaymath}\qquad+
\frac{(1-\beta^2)}{2(1+\beta^2)^2}\left[{\rm arctan}...
...^2)^2}\ln\left[\frac{(y'-\beta)^2}{1+(y')^2}\right].\eqno
(35)\end{displaymath}

The integration constant is found from the condition $y'(0)=\infty$. The details of calculation are given in Sec. 7.

One can see that the $y'>\beta$ condition must be valid. For $x\rightarrow
\infty,~y'\rightarrow\beta$ as

\begin{displaymath}
y'-\beta\propto\exp\left[-\frac{g(1+\beta^2)^2}{\beta C^2}x\right],\eqno(36)\end{displaymath}

as is the case (29). The role of the $\mu$ coefficient is played by the variable $\beta$. Let us consider its physical sense. The velocity v0 of the steady motion along the inclined plane with the inclination angle $\varphi_0$ is equal to

\begin{displaymath}
kv_0=mg\sin\varphi_0.\eqno (37)\end{displaymath}

From (1), (2) we have

\begin{displaymath}v_0=C\cos\varphi_0.\eqno(38)\end{displaymath}

Eliminating v0, we obtain

\begin{displaymath}\tan\varphi_0=kC/mg=\beta.\eqno(39)\end{displaymath}

Expanding the right part of (35) into a series in reciprocal powers of y', we obtain

\begin{displaymath}gx/C=(1/3)(y')^{-3}+\Lambda (y')^{-4}+\dots\eqno(40)\end{displaymath}

This equation differs from (30) with the expression for $\Lambda $ and the change of scale of x. From this one can see that the brachistochrone with drag differs from a cycloid very slightly in the start of falling, so we do not give its plot.

Let us demonstrate that the qualitative behaviour is valid in the case when the forces of friction and of drag F(v), which increases when the velocity grows, affect the particle. The acceleration of particle is equal to

\begin{displaymath}a=g\sin \varphi-[\mu N+F(v)]/m.\eqno(41)\end{displaymath}

At the start of falling the brachistochrone inclination is close to vertical and N is small. The velocity and the drag force are also small. The falling will differ from the movement without drag. Let us suppose that the trajectory of the fastest travel tends asymptotically to the straight line. Then the contribution of the centrifugal forces to N decreases and $N\rightarrow\mu mg\cos\varphi$. Taking into account (1) we obtain for the velocity of the varying of the inclination angle the asymptotic equation

\begin{displaymath}-C\sin\varphi{\rm d\varphi/d}t
\rightarrow g\sin\varphi-\mu g\cos\varphi-F(C\cos\varphi)/m=Q(\varphi).\eqno
(42)\end{displaymath}

The $Q(\varphi)$ increases when $\varphi$ grows and the d$\varphi/$dt derivative decreases. It is obvious that there exists such $\varphi_0$ angle that $Q(\varphi_0)=0$. The brachistochrone inclination will tend to it. It proves our supposition. Near $\varphi_0$ we have $Q(\varphi)\approx\eta(
\varphi-\varphi_0),\linebreak\eta\rm =const$ function and the $\varphi-\varphi
_0$ value will decrease exponentially.


7. Calculations


Let us give the details of the calculation of integral (28) which has the form

\begin{displaymath}
I=\int F(z)dz,\qquad F(z)=\frac{z+\mu}{(\mu-z)(1+z^2)^2}.\eqno(43)\end{displaymath}

We expand F(z) into a sum of elementary fractions

\begin{displaymath}F(z)=\frac{Az+B}{(1+z^2)^2} +\frac C{z-\mu}+\frac{Dz+E}{1+z^2}.\eqno(44)\end{displaymath}

From this

\begin{displaymath}
-(z+\mu)=(C+D)z^4+(E-D\mu)z^3+(A+2C+D-E\mu)z^2
\end{displaymath}


\begin{displaymath}\qquad+(B+E-A\mu-D\mu)z+C-B\mu-E\mu.\eqno(45)\end{displaymath}

Setting coefficients of corresponding powers of z equal, and solving the obtained system of equations, we obtain

\begin{displaymath}A=2\mu/(1+\mu^2),\qquad B=(\mu^
2-1)/(\mu^2+1),\qquad C=-2\mu/(1+\mu^2)^2,
\end{displaymath}


\begin{displaymath}
D=2\mu/(1+\mu^2)^2, \qquad E=2\mu^2/(1+\mu^2)^2.\eqno(46)\end{displaymath}

After expanding (44), integral (43) reduces to an elementary one

\begin{displaymath}
I=(Bz-A)/2(1+z^2)+(E+B/2{\rm )arctan}z\end{displaymath}


\begin{displaymath}
\qquad+C\ln(z-\mu)+(D/2)\ln(1+z^2)\rm +const.\eqno(47)
\end{displaymath}

The integral (35) can be calculated in the same way

\begin{displaymath}
\int{z{\rm d}z
\over(1+z^2)^2(z-\beta)}={1\over1+\beta^2}\i...
...1+z^2)^2}-{\beta
\over1+\beta^2}\int{z{\rm d}z\over(1+z^2)^2}
\end{displaymath}


\begin{displaymath}
\qquad+{\beta\over(1+\beta^2)
^2}\int{{\rm d}z\over z-\beta...
...^2}-{\beta\over(1+\beta^2)^2}\int{{\rm d}z\over1+z^2}.\eqno(48)\end{displaymath}

The numeric calculations were done on the computer using the programming language FORTRAN. When studying the form of the trajectories at the start of movement we used the approximate formula (31) or similar to it. We also used expansion of the integrand near its singularity point z=z0 when calculating the integrals (17) and (18).


References


1. R. Courant, H. Robbins, What is Mathematics?, Oxford University Press, London 1948

2. G. Korn, T. Korn, Mathematical Handbook, McGraw-Hill, New York 1968



This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.47)
Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.