programs/!PROG1.JPETen programs for
exploring dynamics


David Acheson

UPDATE (2012)

When 'From Calculus to Chaos' was published, in 1997, I intended the Ten QBasic programs to run on virtually any PC, without the need for additional software.

Needless to say, a lot has happened since then, and this is no longer the case.

The simplest way I know of running them, now, is to use qb64.

To this end I have taken 6 of the programs, namely
tweaked them slightly for qb64, and put them in a zip file calchaos.

To proceed, download qb64, unzip, extract all the files, and put everything in a folder on your computer called qb64, say. Then download calchaos, unzip and place all six .BAS programs in the same folder as qb64.

To start a program, double click on the main qb64 file (qb64.exe), and open the program you want. By clicking on 'Run', the program should then run in a window. Consult the notes below (and the book) while running the program.

To run the program full screen instead, modify the program by adding _FULLSCREEN near the beginning.
To STOP a program running at any time, if all else fails, press Ctrl-Break.

I include below the original notes about all ten programs, because they are still relevant.
(And some users might conceivably wish to try a different approach, such as DosBox.)


The ten QBasic programs which follow are the core programs in my book

From Calculus to Chaos.

Further information about them, including full listings, can be found in the book, together with an introduction to elementary programming in QBasic.

QBasic is, so far as I know, distributed with every version of Microsoft system software, i.e. Windows 95, Windows NT and the later versions of DOS. If it is not on your hard disk try looking at the installation disks or CD-ROM.

I have tried to make the programs powerful and effective but simple to understand. They make no attempt, then, to compete with sophisticated professional software, and have a really rather different purpose.

To download any of the programs individually, click on a file name, and then save the file (using the File menu in Netscape for example). Alternatively, download all the files at once, in compressed form, by clicking at the bottom of this page.

To stop program execution at any time, press Ctrl-Break. In some programs, pressing q for "quit" will cause a request for a new set of initial conditions to appear on screen.

The first six are general-purpose programs for the numerical integration of systems of ordinary differential equations.


·        1XT.BAS

This program numerically integrates a single differential equation of the form


after plotting the direction field.

It is currently set up for equation (4.1), and can be used to confirm Fig. 4.1.


·        2PHASE.BAS

This program is for an autonomous 2nd-order system


It is currently set up for the equations (5.46a,b) describing simple pendulum motion, and can be used to confirm Fig. 5.17(a).


·        2XTPHASE.BAS This program is for


so the system need not be autonomous. In the autonomous case, however, one merit of this program is that the solution can be seen developing in both the x-t plane and the phase plane simultaneously.

The program is currently set up for the damped linear oscillator (Fig. 5.18).


·        NPHASE.BAS

This program is the backbone for many others in the book, and is for a general nth-order autonomous system.

As it stands, it is set up for the Lorenz equations (11.7), and on inputting r=28, say, we obtain the distinctive `butterfly' of Fig. 11.10.


·        NXT.BAS

This program is similar to NPHASE, but plots one of the variables against time t.

It can be used to investigate

(a) sensitivity to initial conditions and

(b) whether the numerical integration can be `trusted' at all, i.e. whether it produces the same results if the time-step is either halved or doubled. This is because the user is asked, repeatedly, to input (two) initial conditions, a time-step h, and a colour (integer from 1 to 15) in which to present the resulting x-t graph. (To get out of this loop, press Ctrl-Break.)

As it stands, the program is set up for the forced cubic oscillator equation (11.1), converted into an autonomous 3rd-order system. It can be used to produce the chaotic outcome in Fig. 11.1.


·        NVARY.BAS

This program is, again, similar to NPHASE, but it allows us to investigate the effects of gradually changing a parameter which occurs in the equations.

The program is currently set up for the forced Duffing equation (11.14), converted into an autonomous 3rd-order system. The parameters are as in Fig. 11.16.

Try inputting w=1, and when the system has settled into a closed curve, corresponding to a periodic oscillation, press q for "quit". You can then change the parameter w to a new value, and on pressing RETURN the numerical integration will continue from where the system had just got to.

Gradually increasing w from say 1.0 to 1.6 results eventually in a sudden change in the oscillation, and on then gradually decreasing w again the oscillation again jumps back to its original form, but at a lower critical value of w. This is hysteresis, which is quite a common feature of nonlinear systems.

Each of the next three programs is intended for a specific mechanical system.



This program provides an animation of simple pendulum motion, by numerically integrating the relevant equations (see Ex. 5.5). There is frictional damping in the model, and as it stands the damping coefficient in the program is k=0.1 . Try the following input values for the initial angle and the initial angular velocity:


A number of simple variations on PENDANIM are used in the book.



This program provides an animation of a double pendulum which has its pivot vibrating up and down (Sections 12.4, 12.5).

The amplitude and frequency of this up-and-down pivot motion are denoted by a and w respectively. As in PENDANIM, the frictional damping coefficient is set at k=0.1 in the program. The parameter m, currently set as 0.1 in the program, is defined by:


where m1 is the mass of the `first' pendulum bob, attached directly to the pivot, and m2 is the mass of the second.

For chaotic motion, try inputting

a=0.35, w=2, ang1=1, ang2=1

For instability of the unvibrated `inverted' state and consequent `collapse' try

a=0, w=2, ang1=3.1, ang2=3.1

For stabilization of the inverted state by pivot vibration, as in my upside-down pendulum theorem), try

a=0.1, w=20, ang1=2.8, ang2=2.8


·        THREEBP.BAS

This program corresponds to the famous `three-body problem', in which 3 point masses attract one another according to the inverse square law of gravitation (Section 6.8).

programs/!PROGTHR.JPEThe program is essentially NPHASE, with a fixed time step h. The step-by-step method loses accuracy during `close encounters' between masses. One indication of computation accuracy is provided by the total energy, which is continually updated and displayed on screen, and should, in principle, be constant.

As the program stands, the three masses are equal, and all the initial positions and velocities are fixed, except for the initial coordinates of the `third' mass m3.

Try, for example, the following values for x3, y3:

0, 1000
-0.1, .75
0.1, 0.15
0.1, 0.6


The final program is rather different from the rest and is for a partial differential equation.


·        HEAT.BAS

This program uses a step-by-step method to numerically integrate the heat equation (7.16), with a given initial temperature distribution T=f(x), specified in the program.

If we input, say, m=50 (number of grid points in x-direction) and k=.0001 (time step), we see the initial hot spot spread out, quickly at first and more slowly later on.

Warning: inputting a larger value of k than is advised in HEAT will lead to a dramatic numerical instability, which can, in my experience, `crash' the computer system, so be ready to press q for "quit".

Back to Home Page