This document describes several programs which illustrate the kinds of calculations the Analytical Engine might have been called upon to perform. These programs are more complicated than those which appear in the Programming Cards and Mathematical Function Library documents, although some are based upon programs which appear there. All of the following programs are intended to be run in the Web Emulator; they may use features which are not supported in the Command-Line Emulator and require modification in order to work with it. You can run any of these sample programs by clicking the “gears” icon to the right of the text describing the program.

The Fibonacci sequence is an integer sequence defined by the rule that each term in the sequence is the sum of the preceding two terms. The first two terms can be defined as either 1 and 1 or 0 and 1. Using the first definition, the sequence is:

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …

where each term starting with the third is given by:
*F*_{n} =
*F*_{n−1} +
*F*_{n−2}

The Fibonacci series appears in many contexts in nature, such as the arrangements of seeds in a sunflower and leaves on the stems of plants. The ratio of terms in the sequence divided by their immediate predecessor converges upon the golden ratio:

which is found in many aspects of nature and, pleasing to the human eye, commonplace in architecture and design.

The program evaluates the number of terms specified in column 6 of the
store (change the value on the `N006` number card at the
top to compute fewer or more terms), and for the second and successive
terms, the ratio of the term to its predecessor to ten decimal places.
You will see the ratio value quickly converge toward φ, the golden
ratio.

Newton's method (or the Newton-Raphson method) is a numerical technique for finding roots to real-valued functions by making an initial guess which is then refined based upon the value of the function at the guess and its discrepancy with the value sought. Newton's method is one of the most fundamental algorithms in numerical computation.

Here we demonstrate it by evaluating the square root of a
number *N*. We make an initial guess of one half of the number
(we could make a closer guess with a more complicated
expression, but Newton's method converges so quickly there
isn't any reason to go to the trouble), then refine our guess
in subsequent iterations through the process:

and continue to iterate until the estimation of the square root,
*x*, does not change from one iteration to the next to
the precision set for the calculation.

The program takes the square root of the number initially placed in
column `000` of the Store, leaving the result in the same column.
Computation is performed in fixed-point arithmetic, with the number
of decimal places specified by the “A set decimal places”
card at the top of the program. By default, the square root of 2 is
evaluated, but you can compute the square root of any non-negative
number by changing the value on the `N000` number card.
The argument, successive approximations for each iteration, and
result are displayed on the Printer. The result computed for the
square root of two to twenty-five decimal places will be
1.4142135623730950488016887. This value would normally be
rounded to one fewer decimal place, yielding the correct
value to twenty-four decimal places of
1.414213562373095048801689. At the end of the computation a
summary of Engine operations and estimate of the time to complete
them is displayed on the Attendant's Log.

The limaçon, or limaçon of Pascal, is a quartic curve defined in Cartesian coordinates as the locus of solutions to the equation:

(*x*^{2} + *y*^{2} − *ax*)^{2} = *b*^{2}(*x*^{2} + *y*^{2})

where *a* and *b* are parameters which control the size
and shape of the curve. Transforming to polar coordinates, the curve
may be expressed in parametric form as:

*x* = (*b* + *a* cos *t*) cos *t*

*y* = (*b* + *a* cos *t*) sin *t*

It is this form which is used in the program, with the parameter varying from 0 to 2π
to plot the curve. Calculations are performed with 25 decimal places of precision,
using the
Mathematical Function Library to evaluate the sine
and cosine functions at each iteration. The computation is entirely
straightforward and takes no advantage of the curve's symmetry with respect
to the *x* axis. The curve is plotted on the Curve Drawing Apparatus,
after first plotting the *x* and *y* axes in light blue. At
the end of the plot, the Attendant's Log displays the (breathtaking) number
of operations performed and how long the Engine might have taken to
execute them. Parameters *a* and *b* are specified at
the top of the program; try varying them to see how they affect the
plotted curve.

One of the first problems to which mechanical, electromechanical, and electronic computers were applied was calculating tables for naval gunnery and army artillery. Indeed, ENIAC, one of the first electronic computers, was funded by the U.S. Army with that application in mind. The flight of an artillery shell, and hence how a gun must be aimed for it to hit its target, depends in a complex manner upon a multitude of factors: the muzzle velocity with which it leaves the gun, its mass, its size and shape (which determines how much air resistance it will encounter in flight and the degree wind may cause it to diverge from a pure ballistic trajectory), the barometric pressure and humidity of the air through which it is passing, the elevation of the gun's barrel, and more. Naval gunnery further complicates the problem, since both the ship from which the gun is fired and the target may be moving with respect to one another, and pitching and rolling due to rough seas.

So formidable is this problem that in the 20th century naval fire control systems were developed: mechanical or electromechanical analogue computers which could solve the problem in close to real time. Before the advent of such machines, gunners relied upon range tables, compiled by a combination of experiment and manual calculation. Given the importance of naval gunnery to the British Empire, it is almost certain that one of the first tasks for which The Analytical Engine would have been employed was the computation of artillery tables for warships and shore batteries.

Problems such as this are best solved through the process of numerical integration. The flight of the artillery shell is traced along its trajectory from muzzle to target in discrete time steps, with all of the factors affecting its path (momentum, gravitational force, air resistance, etc.) computed at each step. This process is lengthy and tedious to do by hand, but it's ready-made for a computing machine, which can work through the successive steps without human intervention.

This sample program illustrates numerical integration by modeling the trajectory of a projectile fired from a British BL 15 inch Mark I naval gun, which served the Royal Navy from 1915 through 1959, including both World Wars. As a sample program, the computation is much simpler than would have been used to create real artillery tables, but the basic principle of numerical integration is the same.

We start by defining the parameters of the problem, physical constants, and options for the computation as number cards at the top of the program. The following table lists these parameters, giving the symbol by which they are referred to below and the column in the Store into which each is placed. Consistent with the Victorian vintage of The Analytical Engine, these definitions and all calculations in the program will be done in quaint imperial units.

From the calibre (inside barrel diameter) *c* of the gun, we
can calculate the frontal area *a* of the projectile, which will figure
in the calculation of atmospheric drag, from the formula for the area
of a circle.

At each step *n* in the calculation, we start with the current velocity of
the projectile *v _{n}*, its distance downrange

We start by calculating the effect of air resistance, *d*. This
is proportional to the product of the air density ρ, the
coefficient
of drag *C _{d}* (a measure of how the projectile's
shape affects the drag it will encounter),
the frontal area of the projectile

The horizontal and vertical velocity must be treated separately
since the acceleration of gravity only acts upon the
vertical component of velocity, decreasing it over time
*s* by *g s*. Taking the square root of
the sum of the squares of the horizontal and gravity-affected
vertical velocities gives the new velocity
*v*_{n+1} at the end of the time step,
and the arc tangent of the vertical velocity divided by the
horizontal velocity gives *θ*_{n+1}, its
new angle of motion. Multiplying the horizontal and vertical
velocities by the time step and adding them respectively to
the range and altitude updates these quantities:
*x*_{n+1} and *y*_{n+1}.
If the new altitude is negative, the projectile has impacted and we're
done. Otherwise, we continue for another step, starting with the
updated position, velocity, and angle just computed.

At the end of the computation, a summary appears on the Printer, both
with and without the effects of atmospheric drag. The latter is not
calculated because the Royal Navy anticipated naval battles on the
Moon, but to illustrate the degree air resistance affects the flight
of an artillery shell. Try increasing the coefficient of drag on the
`N104` card to see how it reduces the range and altitude reached
by the projectile. A massive projectile
moving faster than the speed of sound delivers a large amount of
kinetic energy
to its target, quite apart from any explosives with which it
may have been filled. This kinetic energy is calculated in
“foot pounds” (another example of sloppy
terminology in imperial units, users of which often fail to
make the distinction between weight and mass—the correct
term is “foot lbf” where “lbf” stands
for “pounds of force”). The energy is then compared
to the quantity in pounds of the
high explosive TNT
which would release the same energy.

Computation with atmospheric drag Time of flight: 64.0 seconds Range: 32030.1 yards Max altitude: 5439.1 yards Velocity at impact: 1468.6 ft/sec Impact kinetic energy: 64957619.1 foot pounds (47.3 pounds of TNT) Computation without atmospheric drag Time of flight: 76.0 seconds Range: 53926.8 yards Max altitude: 7416.6 yards Velocity at impact: 2482.6 ft/sec Impact kinetic energy: 185629564.4 foot pounds (135.2 pounds of TNT)

The Curve Drawing Apparatus will display the trajectory of the
projectile, both taking into account air resistance (the black
curve) and without air resistance (blue curve).
Note how increasing the drag
coefficient (the `N104` number card at the top of the
program) affects the trajectory, including
deviation from the parabola expected absent the effects of
atmospheric drag. With very high drag, the projectile almost comes
to a stop and falls from the sky

Once we have a deterministic and high-precision computing device such as The Analytical Engine, it may seem distinctly odd, if not somewhat scandalous, to use it to generate random numbers, but there are a number of applications in mathematics, statistics, and the sciences where we wish to randomly sample, for example, in computations using the Monte Carlo method.

A pseudorandom sequence generator is an algorithm which, although
deterministic, produces a sequence of outputs which behave like
genuinely random sequences produced through natural processes such as
radioactive decay.
This program implements and demonstrates a pseudorandom generator
which uses two multiplicative congruential generators whose values
are combined to increase the period of the results. The generators
are “seeded” by a value initially supplied in column
0 of the store when the generator is initialised by including
the cards for `prand_init` from the library; this may be
any number between 1 and 2147483647. Subsequently, each time the
`prand_next` cards from the library are executed, a fixed point
pseudorandom number between 0 and 1, with the number of decimal places
specified by the “`A set decimal places to `*n*”
card, will be returned in column 0 of the store.

The sample program uses the library cards to generate and print
50 (change the `N011` card at the top of the program if you prefer a
different number) pseudorandom values and performs one of the simplest
randomness tests on them: computing their mean value. For a uniformly
distributed sequence of random values between 0 and 1, the mean value
should approach 0.5 as the length of the series increases. The seed
value for the generator is specified by the `N000` card
at the top; changing it will produce entirely different values in
the sequence, but have little effect on the mean value. The technique
used to generate the pseudorandom numbers has around 10 decimal places
of precision: nothing will be gained by setting the fixed decimal point
to a larger value.

Now that we have a pseudorandom number generator, let's take it to the casino and have some fun. Using the generator introduced in the previous example, it is possible to (very, very slowly) approximate the value of π.

Imagine a square with edges of length *l* = 2 (which happens
to be the shape of the drawing surface of the Engine's
Curve Drawing Apparatus). Now inscribe a circle of radius *r*
within this square; we'll draw it with the centre of the circle at
the square's centre to simplify the calculations, but it doesn't
matter where it's drawn. Now plot a number *n* of points at random
locations within the square and count how many points *i* fall inside
the circle. The area of the circle is
*a* = π*r*², so the probability
any given random point will fall within the circle is
*p* = *a*/*l*² = π*r*²/4.
A little algebra lets us solve for π = 4*p*/*r*².
From a number of random points, we can estimate the probability of points
falling within the circle as
*p* ≅ *i*/*n*. Finally, we can
estimate π ≅ (4 *i*/*n*)/*r*².

The sample program uses the pseudorandom generator to produce random
coordinate pairs in the range ±1, determines whether each point
falls inside or outside the circle (its radius, 0.5 by default, may
be changed by altering the value on the `N106` number card
at the top), and counts points inside and total points. Periodically
(the frequency it set by `N120`), an estimate for π is
calculated as described above and printed along with the number of
points generated so far. The circle and points are plotted on the
Curve Drawing Apparatus, with points inside the circle in green and
those outside blue. The plot is cleared and restarted every 2500
(set by `N122`) points due to limitations in the
implementation of the Curve Drawing Apparatus. Set `N150` to
0 to suppress the plotting of points entirely.

Here are the first ten lines output from the program, showing the estimates for the first thousand random points. The first column is the number of random points plotted, the second the estimate of π, and the third the error of the estimate with regard to the correct value.

100: 3.6800000000 17.13% 200: 3.5200000000 12.04% 300: 3.6800000000 17.13% 400: 3.6800000000 17.13% 500: 3.5200000000 12.04% 600: 3.6266666664 15.44% 700: 3.5428571428 12.77% 800: 3.4000000000 8.22% 900: 3.2711111108 4.12% 1000: 3.2160000000 2.36%

As the number of points crosses 30,000, the error in the estimated value has declined substantially. Because this is a random, statistical process and the resolution of the generator and plot are finite, the estimate will continue to oscillate around the correct value as more and more points are plotted.

29700: 3.1380471380 0.11% 29800: 3.1366442952 0.15% 29900: 3.1363210700 0.16% 30000: 3.1365333332 0.16% 30100: 3.1394019932 0.06% 30200: 3.1401324500 0.04% 30300: 3.1419141912 0.01% 30400: 3.1415789472 0.00% 30500: 3.1412459016 0.01% 30600: 3.1393464052 0.07%

Several other documents related to The Analytical Engine contain sample programs. Here are links to them.

- Mathematical Function Table in The Mathematical Function Library
- Inverses and Identities in The Mathematical Function Library
- Plot Cubic Curve in Programming Cards
- Commercial Computing in Sterling in Programming Cards

by John Walker