Is continuity inherently more powerful than discreteness?
(Or how I built my first (proper) computer...)
Tuesday, June 9, 2015
Sunday, June 7, 2015
Coupled FitzHugh - Nagumo system and the Fuss - Catalan numbers
After some mathematical noodling (whilst awaiting a replacement oscilloscope), I've come up with two bizarre expressions for the exact equilibrium position of the coupled FitzHugh-Nagumo attractor (as per previous posts), for the cases in which coupling parameter k = -1 and 0 respectively:
C3(n) is the nth Fuss-Catalan number of order 3, or ternary number [1], the first few terms being 1, 1, 3, 12, 55, 273, 1428, 7752...
Given x1, we can readily obtain y1 from (as per previous post):
Example
Using the above equations for x1, with a = 0.4, b= 0.5, c = 0.01, as per Sprott's figure 6.7, we calculate
x1 = -0.7901341941135... (k = -1 case), and
x1 = -0.8056544203709... (k = 0 case).
These agree with the values for x1 given by equations (7) and (8) in the previous post.
Reference
[1] Sequence A001764, N. J. A. Sloane, A Handbook of Integer Sequences, New York and London: Academic Press (1973). Online at https://oeis.org
C3(n) is the nth Fuss-Catalan number of order 3, or ternary number [1], the first few terms being 1, 1, 3, 12, 55, 273, 1428, 7752...
Given x1, we can readily obtain y1 from (as per previous post):
Example
Using the above equations for x1, with a = 0.4, b= 0.5, c = 0.01, as per Sprott's figure 6.7, we calculate
x1 = -0.7901341941135... (k = -1 case), and
x1 = -0.8056544203709... (k = 0 case).
These agree with the values for x1 given by equations (7) and (8) in the previous post.
Reference
[1] Sequence A001764, N. J. A. Sloane, A Handbook of Integer Sequences, New York and London: Academic Press (1973). Online at https://oeis.org
Coupled FitzHugh - Nagumo system - equilibrium position
The coupled FitzHugh-Nagumo oscillators produce an attractor (previous post) which has an obvious resting / equilibrium point. The trajectory spends a long time spiralling about this point before making one or more fast excursions, before returning to the equilibrium point, and so on.
The following shows the location of the centre of the equilibrium region for Sprott's figure 6.7 [1]. I've added my estimates as to the coordinates of the equilibrium centre:
Given that the trajectory moves slowly around the equilibrium point, it seems reasonable that one could work out an analytic solution as to the whereabouts of this point. In pursuit of this, I took Sprott's coupled FitzHugh-Nagumo system (equation (6.7)), and set the right hand sides to zero (i.e. the game being to find where the speed of the trajectory is zero):

Note
The ninth order (degree) polynomial appears to have just one real root. This is equation (7) above.
Reference
[1] Sprott J.C. Elegant Chaos: Algebraically Simple Chaotic Flows, pp. 123 - 125.
The following shows the location of the centre of the equilibrium region for Sprott's figure 6.7 [1]. I've added my estimates as to the coordinates of the equilibrium centre:
Given that the trajectory moves slowly around the equilibrium point, it seems reasonable that one could work out an analytic solution as to the whereabouts of this point. In pursuit of this, I took Sprott's coupled FitzHugh-Nagumo system (equation (6.7)), and set the right hand sides to zero (i.e. the game being to find where the speed of the trajectory is zero):
After a fair bit of algebra, I worked out the following rather compact equations for the coordinates of the equilibrium position (see box below). These equations are valid for small values of c - consistent with Sprott's comment that 'chaotic solutions require that c be rather small' [1]. These equations were derived from a power series development of an exact solution for the coordinates, which was worked out, the derivation of which is given also, below.
It is evident that the y-coordinate does not depend on parameter c, when c is small.
The following tables gives the calculated coordinates of the equilibrium point, for k = -1 and k = 0 (no coupling). All other values are as per Sprott's figure 6.7 (i.e. a = 0.4, b = 0.5, c = 0.01). The first table gives the results obtained using the exact solution, the second gives results obtained using the above, boxed, equations.
| Coordinates of equilibrium point calculated using exact derivation. |
![]() |
| Coordinates of equilibrium point calculated using approximate derivation. |
These agree rather nicely with the values I measured directly off Sprott's figure 6.7 - namely -0.78 for x1 and +0.51 for y1.
What's interesting (to me) is that the equilibrium point actually depends on the value of the coupling parameter k. I (perhaps naively) assumed that the equilibrium point for a single decoupled oscillator would be the same as that for the coupled case - but experiments with the analog computer indicated otherwise:
![]() |
| Analog computer results for coupled (k = -1) and uncoupled (k = 0) cases. It is evident that the equilibrium points have very different y-values. |
In particular the y-coordinate of the equilibrium point varies strongly with k.
In contrast, the x-coordinate of the equilibrium position does not depend strongly on the coupling parameter k. To first order, the x-coordinate is dominated by the -a/b term which, in this case is equal to -0.8.
For reasons which remain unexplored (offset voltages?), the analog computer's results are slightly different to Sprott's. My equilibrium positions are shifted with respect to the origin of the plot.
![]() |
Coordinates of equilibrium point measured from
analog computer output plot.
|
Derivation of Exact Coordinates of Equilibrium Point

The ninth order (degree) polynomial appears to have just one real root. This is equation (7) above.
Reference
[1] Sprott J.C. Elegant Chaos: Algebraically Simple Chaotic Flows, pp. 123 - 125.
Sunday, May 31, 2015
Coupled FitzHigh - Nagumo System
During the 1950s Hodgkin and Huxley came up with a set of differential equations which mimic the physics of neural activity - based on experiments with the giant squid axon.
In essence, a neuron doesn’t fire until a threshold input stimulus is reached and, after firing, neurons exhibit a relaxation time, which prevents them from firing again for a short time.
A simplified version of the (four dimensional) Hodgkin-Huxley model is the (two dimensional) Fitzhugh -Nagumo [1] oscillator:
In the above, x is the neuron membrane voltage, a, b, c and k are constants and y is the 'recovery variable' for the neuron membrane potential. The membrane voltage (x) has a cubic non-linearity.
Fitzhugh-Nagumo oscillators can be coupled together - Sprott considers the following [2] coupling of two Fitzhugh-Nagumo oscillators - the stimulus for the first is some constant (k) times the membrane voltage of the second (x2); the stimulus for the second is some constant (k) times the membrane voltage of the first (x1):
Note: there appears to be a typographical error in Sprott's equation (6.10) - the third equation should (I think) be as I've written above - i.e. y2 rather than y1 appearing on its right hand side.
Here's my patching for the coupled equations:
The following plots were obtained using the parameter values / initial conditions given by Sprott [2] (i.e. a = 0.4, b = 0.5, c = 0.01 and k = -1). The plots show the first ten minutes of the trajectories - using the 0.2 µF integration capacitors throughout. The red plot shows y2 vs. x2; the black plot shows y1 vs. x1.
With the computer in FAST mode I captured the following (y1 vs. x1)...just before my Tektronix 465 blew up(!).
It's interesting to watch the trajectory unfold with the computer running in slow time (0.2 µF integration capacitors = 0.2 second time constant), on the plotter. The pen whizzes around the perimeter from the initial condition position, settling down at (almost) a point (towards bottom left - resting state - bright splodge on oscilloscope display); the trajectory then slowly spirals outwards (initially very slowly) from this resting point, eventually gaining speed (and eventually reaching the firing threshold / stimulus) and going around the perimeter again (quickly), and so on.
Sometimes the fast part of the trajectory follows the outer path / orbit, followed immediately by a second slightly inner fast path / orbit, before returning to the 'resting spiral' (equilibrium resting potential). The fast moving outer parts of the trajectory corresponds to the neuron firing; the spiral part corresponds to the equilibrium / resting part. I guess the second fast orbit which is, on occasion, observed corresponds to a relative refractory (see below) part of the trajectory.
Reference [3] gives a good depiction of the trajectory and its elements, reproduced below. The neuron fires (active part of trajectory); immediately followed by a part in which it cannot be stimulated, no matter how great a stimulus is applied (absolute refractory), followed by a 'relative refractory' part (in which a second firing is possible - but only given enough stimulus) - followed by the resting phase.
References
[1] FitzHugh R. (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophysical J. 1:445-466.
[2] Sprott J.C. Elegant Chaos: Algebraically Simple Chaotic Flows, pp. 123 - 125.
[3] http://www.scholarpedia.org/article/FitzHugh-Nagumo_model
In essence, a neuron doesn’t fire until a threshold input stimulus is reached and, after firing, neurons exhibit a relaxation time, which prevents them from firing again for a short time.
A simplified version of the (four dimensional) Hodgkin-Huxley model is the (two dimensional) Fitzhugh -Nagumo [1] oscillator:
In the above, x is the neuron membrane voltage, a, b, c and k are constants and y is the 'recovery variable' for the neuron membrane potential. The membrane voltage (x) has a cubic non-linearity.
Fitzhugh-Nagumo oscillators can be coupled together - Sprott considers the following [2] coupling of two Fitzhugh-Nagumo oscillators - the stimulus for the first is some constant (k) times the membrane voltage of the second (x2); the stimulus for the second is some constant (k) times the membrane voltage of the first (x1):
Note: there appears to be a typographical error in Sprott's equation (6.10) - the third equation should (I think) be as I've written above - i.e. y2 rather than y1 appearing on its right hand side.
Here's my patching for the coupled equations:
![]() |
| Richard FitzHugh with analog computer, ca. 1960, from [3]. |
| Analog computer, ca. 2015. Patched for coupled FitzHugh - Nagumo system. |
The following plots were obtained using the parameter values / initial conditions given by Sprott [2] (i.e. a = 0.4, b = 0.5, c = 0.01 and k = -1). The plots show the first ten minutes of the trajectories - using the 0.2 µF integration capacitors throughout. The red plot shows y2 vs. x2; the black plot shows y1 vs. x1.
![]() |
| 10 minutes' worth of trajectory. Green dots mark start / initial condition point in each case (red pen took time to get going!). |
With the computer in FAST mode I captured the following (y1 vs. x1)...just before my Tektronix 465 blew up(!).
| Attractor for the coupled FitzHugh - Nagumo system, as per [2], figure 6.7. (1 V/div both axes.) |
It's interesting to watch the trajectory unfold with the computer running in slow time (0.2 µF integration capacitors = 0.2 second time constant), on the plotter. The pen whizzes around the perimeter from the initial condition position, settling down at (almost) a point (towards bottom left - resting state - bright splodge on oscilloscope display); the trajectory then slowly spirals outwards (initially very slowly) from this resting point, eventually gaining speed (and eventually reaching the firing threshold / stimulus) and going around the perimeter again (quickly), and so on.
Sometimes the fast part of the trajectory follows the outer path / orbit, followed immediately by a second slightly inner fast path / orbit, before returning to the 'resting spiral' (equilibrium resting potential). The fast moving outer parts of the trajectory corresponds to the neuron firing; the spiral part corresponds to the equilibrium / resting part. I guess the second fast orbit which is, on occasion, observed corresponds to a relative refractory (see below) part of the trajectory.
Reference [3] gives a good depiction of the trajectory and its elements, reproduced below. The neuron fires (active part of trajectory); immediately followed by a part in which it cannot be stimulated, no matter how great a stimulus is applied (absolute refractory), followed by a 'relative refractory' part (in which a second firing is possible - but only given enough stimulus) - followed by the resting phase.
![]() |
| Phase portrait and physiological state diagram of FitzHugh-Nagumo model, taken from [3]. |
References
[1] FitzHugh R. (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophysical J. 1:445-466.
[2] Sprott J.C. Elegant Chaos: Algebraically Simple Chaotic Flows, pp. 123 - 125.
[3] http://www.scholarpedia.org/article/FitzHugh-Nagumo_model
Wednesday, May 27, 2015
100 weeks
Almost exactly 100 weeks since I started, I've finished. Never again!
| Baltic birch plywood side panels. Recessed handles by Mentor. |
| Rear of the computer. Two 120 mm fans at top. |
| Funny - doesn't look like a USB connector - because it isn't - this allows a second computer (or, more realistically, part thereof) to be connected... |
| Hour counter, IEC mains connector and fuse. |
Sunday, May 17, 2015
Conjecture: A tilted circular billiard table is chaotic
Here we have a point-like ball moving inside a (soft, elastic) circular boundary i.e. a circular billiard table. Interestingly, Charles L. Dodgson, (Lewis Carroll of Alice's Adventures in Wonderland fame), once published a set of rules for a two-player game of circular billiards [1]. However, this is predated somewhat by Alhazen / Ptolemy [2].
This post follows the same approach as the previous post, but with the V-wedge replaced with a full circle of radius R:
q is the distance from a point (x,y) to the circular surface, measured along the normal to the surface.
Hence if a point-like ball is at a point (x,y) outside the circular perimeter, there is an elastic restoring force (acting towards the centre of the circle) given by kq, where k is the spring constant of the wall of the circle. It is easy enough to work out q, and thus the accelerations along x and y:
where m is the mass of the ball.
Acceleration due to gravity is directed purely along the negative y-axis. We assume that the radius of the circle is small compared to the radius of the earth (i.e. g is constant).
The computer's squarers yield the square of the input voltage divided by ten - both quantities Q and q end up divided by the square root of ten; this comes out in the wash of the divider, however. The only point to remember is that the voltage corresponding the the radius R must also be divided by √10 also (before it gets subtracted from Q/√10 to yield q/√10). Here's the patching:
No gravity case
We simply set g = 0.
Now the ball doesn't move (there is no initial force). However if we set the initial velocity to something non-zero the point-like ball moves across to the circular wall, and rebounds. This is most interesting perhaps if we start the ball from someplace other than the origin - otherwise the ball simply oscillates back and forth through the origin. The following plots show the trajectory for the ball starting at (1,1) with an initial velocity of 3 m/s along the positive x-axis. -R√10 is set to -0.632 V, corresponding to a billard table radius R of 2 (metres).
In all cases, the trajectory was plotted for eight minutes (integrators have 1 second time constants). Different coloured plots are for different wall spring constants: (a) k = 0.2 N/m (black line), (b) k = 1 N/m (red line), (c) k = 10 N/m (these spring constants assume a 1 kg mass).
The results indicate that is indeed not possible for the trajectory of a moving ball within a circular billiard table to cover the whole table completely. Unless of course...
Gravity Case
...we tilt the table. In other words we have a non-zero gravitational force (here along the negative y-axis).
We start with the ball at (1,1) as before, but with zero initial velocity. Acceleration due to gravity (g) is 10 m/s2. The following result is for k = 10 N/m (eight minutes of trajectory).
A lesser tilt also gives extensive coverage of the (lower half) of the table...here g = 2 m/s2 (about an 11° slope). Everything else as before:
An even lesser tilt (g = 0.5 m/s2 - about a 2.8° slope on planet Earth) still yields good coverage - albeit at a slower pace:
Finally, g = 0.1 m/s2 (about a 0.57° slope on planet Earth - and surely hardly noticeable in a darkened pool hall) yields the result below. This result is suspect to some degree however, since the ball passes higher than its starting height - suggests some issue with the computation - e.g. an offset null somewhere (e.g. squarer outputs), in effect inputting energy into the system:
References
[1] Dodgson, Charles Lutwidge ('Lewis Carroll'). Circular Billiards for two players. 1890.
[2] Alhazen's problem corresponds to finding the point on the edge of a circular billiard table at which a cue ball at a given point must be aimed in order to bounce off the edge of the table and strike another ball at a second given point, named after the Arab scholar Abū ʿAlī al-Ḥasan ibn al-Ḥasan ibn al-Haytham (c. 965 AD – c. 1040 AD), although it was first formulated by Ptolemy even earlier, in 150 AD.
This post follows the same approach as the previous post, but with the V-wedge replaced with a full circle of radius R:
q is the distance from a point (x,y) to the circular surface, measured along the normal to the surface.
Hence if a point-like ball is at a point (x,y) outside the circular perimeter, there is an elastic restoring force (acting towards the centre of the circle) given by kq, where k is the spring constant of the wall of the circle. It is easy enough to work out q, and thus the accelerations along x and y:
where m is the mass of the ball.
Acceleration due to gravity is directed purely along the negative y-axis. We assume that the radius of the circle is small compared to the radius of the earth (i.e. g is constant).
The computer's squarers yield the square of the input voltage divided by ten - both quantities Q and q end up divided by the square root of ten; this comes out in the wash of the divider, however. The only point to remember is that the voltage corresponding the the radius R must also be divided by √10 also (before it gets subtracted from Q/√10 to yield q/√10). Here's the patching:
No gravity case
We simply set g = 0.
Now the ball doesn't move (there is no initial force). However if we set the initial velocity to something non-zero the point-like ball moves across to the circular wall, and rebounds. This is most interesting perhaps if we start the ball from someplace other than the origin - otherwise the ball simply oscillates back and forth through the origin. The following plots show the trajectory for the ball starting at (1,1) with an initial velocity of 3 m/s along the positive x-axis. -R√10 is set to -0.632 V, corresponding to a billard table radius R of 2 (metres).
In all cases, the trajectory was plotted for eight minutes (integrators have 1 second time constants). Different coloured plots are for different wall spring constants: (a) k = 0.2 N/m (black line), (b) k = 1 N/m (red line), (c) k = 10 N/m (these spring constants assume a 1 kg mass).
The results indicate that is indeed not possible for the trajectory of a moving ball within a circular billiard table to cover the whole table completely. Unless of course...
Gravity Case
...we tilt the table. In other words we have a non-zero gravitational force (here along the negative y-axis).
We start with the ball at (1,1) as before, but with zero initial velocity. Acceleration due to gravity (g) is 10 m/s2. The following result is for k = 10 N/m (eight minutes of trajectory).
A lesser tilt also gives extensive coverage of the (lower half) of the table...here g = 2 m/s2 (about an 11° slope). Everything else as before:
An even lesser tilt (g = 0.5 m/s2 - about a 2.8° slope on planet Earth) still yields good coverage - albeit at a slower pace:
Finally, g = 0.1 m/s2 (about a 0.57° slope on planet Earth - and surely hardly noticeable in a darkened pool hall) yields the result below. This result is suspect to some degree however, since the ball passes higher than its starting height - suggests some issue with the computation - e.g. an offset null somewhere (e.g. squarer outputs), in effect inputting energy into the system:
References
[1] Dodgson, Charles Lutwidge ('Lewis Carroll'). Circular Billiards for two players. 1890.
[2] Alhazen's problem corresponds to finding the point on the edge of a circular billiard table at which a cue ball at a given point must be aimed in order to bounce off the edge of the table and strike another ball at a second given point, named after the Arab scholar Abū ʿAlī al-Ḥasan ibn al-Ḥasan ibn al-Haytham (c. 965 AD – c. 1040 AD), although it was first formulated by Ptolemy even earlier, in 150 AD.
Thursday, May 14, 2015
Gravitational Wedge Billiard
Here we have a V-shaped barrier in which a point mass (ball) bounces about, under the influence of a constant vertically downward force (e.g. gravity). There seemingly being nothing new in chaos, this is (perhaps well) known as the gravitational wedge billiard [1].
We consider a right angle wedge as shown below:
As in the previous post, I've assumed that if the ball is below slope there is an outward (towards positive y-axis) restoring force along the normal to the slope. As before, this force is given by Hooke's law. In the previous post a value for q (the normal distance from ball to slope) was worked out, applicable to the right hand side (i.e. x > 0). Extending the thing to both sides of the wedge gives the result:
and, as before m is the mass of the ball and k the spring constant. Hence the above gives the motion of the ball as it hits, descends into, and eventually bounces off the slope. The bounce is elastic. The bigger the value of k, the less spongy the slope is. Here's the patching:
Each of the following results shows the first four minutes' worth of the trajectory (1 second time constant on the integrators), with the ball starting at (0,6). Acceleration due to gravity (acting vertically downwards) is 10 m/s2. The quoted spring constants correspond to a 1 kg mass. Note that, as spring constant k is reduced, the ball sinks further below the x axis, at least for the first few bounces. Although, in theory, one would expect the ball to simply bounce up and down along the y-axis forever, in practice it deviates to one side or the other after two or three bounces, enough to flick it on to one of the sides of the wedge, high up enough to lead to further deviations and ensuing chaos.
An interesting observation: the 'effective' wedge angle reduces as the spring constant reduces: the more spongy the slope, the smaller the wedge angle. What's also interesting / unexpected (perhaps) is that the sides of the triangular wedge trajectory space remains sharply defined.
Finally, if, rather than dropping the ball straight down the wedge's vertical axis, we launch it from part way up the left hand side, towards the right hand side, we achieve a more stable trajectory / orbit. For example, with the initial velocities along x and y as +8 m/s and +6 m/s respectively, and starting at (-4,4) we get:
This looks something like a precessing(?) period two orbit.
References
[1] Classical and quantum chaos of the wedge billiard. I. Classical mechanics, T. Szeredi and D. A. Goodings, Physical review. E, November 1993; 48(5).
We consider a right angle wedge as shown below:
As in the previous post, I've assumed that if the ball is below slope there is an outward (towards positive y-axis) restoring force along the normal to the slope. As before, this force is given by Hooke's law. In the previous post a value for q (the normal distance from ball to slope) was worked out, applicable to the right hand side (i.e. x > 0). Extending the thing to both sides of the wedge gives the result:
and, as before m is the mass of the ball and k the spring constant. Hence the above gives the motion of the ball as it hits, descends into, and eventually bounces off the slope. The bounce is elastic. The bigger the value of k, the less spongy the slope is. Here's the patching:
Each of the following results shows the first four minutes' worth of the trajectory (1 second time constant on the integrators), with the ball starting at (0,6). Acceleration due to gravity (acting vertically downwards) is 10 m/s2. The quoted spring constants correspond to a 1 kg mass. Note that, as spring constant k is reduced, the ball sinks further below the x axis, at least for the first few bounces. Although, in theory, one would expect the ball to simply bounce up and down along the y-axis forever, in practice it deviates to one side or the other after two or three bounces, enough to flick it on to one of the sides of the wedge, high up enough to lead to further deviations and ensuing chaos.
![]() |
| k = 200 N/m |
![]() |
| k = 20 N/m |
![]() |
| k = 10 N/m |
Finally, if, rather than dropping the ball straight down the wedge's vertical axis, we launch it from part way up the left hand side, towards the right hand side, we achieve a more stable trajectory / orbit. For example, with the initial velocities along x and y as +8 m/s and +6 m/s respectively, and starting at (-4,4) we get:
References
[1] Classical and quantum chaos of the wedge billiard. I. Classical mechanics, T. Szeredi and D. A. Goodings, Physical review. E, November 1993; 48(5).
Subscribe to:
Posts (Atom)


























