Ten programs for
    exploring dynamics

    by

    David Acheson





    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

      dx/dt=f(x,t)

    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

      dx/dt=f(x,y)
      dy/dt=g(x,y)

    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

      dx/dt=f(x,y,t)
      dy/dt=g(x,y,t)

    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.


  • PENDANIM.BAS

    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:

      .2,0
      3,0
      3.1415,0
      0,2
      0,4

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


  • PENDOUBL.BAS

    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:

      m=m2/(m1+m2)

    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).

    The 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".


    Click here to download all the programs at once , in compressed (zip) form.



    Back to Home Page