ALINK="#FF0000">
LINUX GAZETTE
[ Prev ][ Table of Contents ][ Front Page ][ Talkback ][ FAQ ][ Next ]

"Linux Gazette...making Linux just a little more fun!"


Numerical Workbenches

By Christoph Spiel


Some people think GNU/Linux is a good operating system, but has not enough applications to make it succeed in the market. Although this might be true for the desktop area, it is certainly wrong for numerical workbenches. In this field GNU/Linux users have many different (and excellent) choices -- in fact too many to introduce them all. Therefore, this series of articles introduces three outstanding applications:

GNU/Octave 2.1.34
http://www.che.wisc.edu/octave/
Scilab 2.6
http://www-rocq.inria.fr/scilab/
Tela 1.32
http://www.geo.fmi.fi/prog/tela.html

To find out about more numerical workbenches, check out http://sal.kachinatech.com/A/2/

Introduction

What can these programs do? Isn't paper and pencil -- er -- a spreadsheet program enough?

The main application areas of numerical workbenches are:

However, because all of them provide complete programming languages to the user and, moreover, are designed to be extended, the number of numerical problems they can solve is almost limitless.

Numerical Mathematics

Now, what the heck is numerical math anyhow? Numerical Mathematics is the branch of math that develops, analyzes, and applies methods to compute with finite precision numbers. Computer hardware, for example, uses numerical math.

Why do computers work with finite precision numbers? Why has nobody developed a scheme that allows for the storage of exact numbers?

  1. Solving problems with finite precision numbers is faster -- much, much faster. Take for example the sum of all square roots from one to one million. On my computer, doing the exact computation with MuPAD-1.4.2 available at http://math-www.uni-paderborn.de/MuPAD/index.html
        time(sum(sqrt(i), i = 1..10^6));
    

    takes about 40 seconds, whereas getting the approximate result with Tela-1.32

        tic(); sum(sqrt(1:10^6)); toc();
    

    takes 0.31 seconds, that is, the answer in finite precision is returned over 100 times faster! Put another way, we can crunch hundred times more data with finite precision numbers in the same time slice.

  2. When using good algorithms -- the ones suggested by numerical mathematicians -- and being careful one can get surprisingly precise answers even with finite precision numbers.

  3. Admit it, most of the time users do not need exact results! A good approximation -- with ``sufficiently many correct digits'' -- will do.

Article Organization

In this article series, we point out the similarities among the three applications that we are going to discuss. We will use GNU/Octave in most of the examples. Where there are important differences you should be aware of, we have put a Differences paragraph at the end of the section.

Technical details for the terminally curious have been put in Details sections.

Getting In and Out

To give you a hands-on experience, let us start each of the applications, request help on a function, and then quit.

GNU/Octave
    cspiel@hydra:~/articles/numerical-workbenches $ octave
    GNU Octave, version 2.1.34 (i686-pc-linux-gnu).
    Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 John W. Eaton.
    This is free software with ABSOLUTELY NO WARRANTY.
    For details, type `warranty'.
    *** This is a development version of Octave.  Development releases
    *** are provided for people who want to help test, debug, and improve
    *** Octave.
    ***
    *** If you want a stable, well-tested version of Octave, you should be
    *** using one of the stable releases (when this development release
    *** was made, the latest stable version was 2.0.16).
    octave:1> help diag
    diag is a built-in function
     - Built-in Function:  diag (V, K)
         Return a diagonal matrix with vector V on diagonal K.  The second
         argument is optional.  If it is positive, the vector is placed on
         the K-th super-diagonal.  If it is negative, it is placed on the
         -K-th sub-diagonal.  The default value of K is 0, and the vector
         is placed on the main diagonal.  For example,
              diag ([1, 2, 3], 1)
              =>  0  1  0  0
                       0  0  2  0
                       0  0  0  3
                       0  0  0  0
    octave:2> quit
    cspiel@hydra:~/articles/numerical-workbenches $

Alternatively use exit or press C-d to quit GNU/Octave.

GNU/Octave offers the user function-name completion, this is, when only part of a function's name is entered and the user hits Tab, the partial name is completed as much as possible. A second Tab displays the list of remaining choices.

Scilab
After starting Scilab thus:
    cspiel@hydra:~/articles/numerical-workbenches $ scilab

we get a new X-window in which the Scilab interpreter runs. Asking for help opens an xless(1x) window. (Both these links are screenshots, so click on them.)

To exit Scilab, enter quit or exit.

Scilab can also be launched in non-window mode by passing the -nw option to it:

    cspiel@hydra:~/articles/numerical-workbenches $ scilab -nw
                               ===========
                               S c i l a b
                               ===========
                              scilab-2.6
                      Copyright (C) 1989-2001 INRIA
    Startup execution:
      loading initial environment
    -->help diag

The help system then uses the text output, too.

Tela
Tela's banner is quite terse, nonetheless, the help system is as comprehensive as necessary. Note that Tela offers function name completion as GNU/Octave does.
    cspiel@hydra:~/articles/numerical-workbenches $ tela
    This tela is a tensor language, Version 1.32.
    Type  ?help  for help.
    ->TAB completion works; try docview() and source("demo")
    >help diag
    diag(V, K) (V is a vector) returns a square diagonal matrix, with
       vector V on the main diagonal (K == 0, default), the K-th super
       diagonal (K > 0) or the K-th sub-diagonal (K < 0).
       diag(M, K) (M is a matrix) returns the main diagonal (K == 0,
       default), the K-th super diagonal (K > 0), or the K-th sub-diagonal
       (K < 0) of M as a vector.  M need not necessarily be square.
    >quit()
    63 instructions, 0 arithmetic ops.
    0.0063 MIPS, 0 MFLOPS.
    cspiel@hydra:~/articles/numerical-workbenches $

Tela can also be exited by pressing C-d.

Better Than a Pocket Calculator!

Now that we know how to start and exit the programs, let us look at them in action.

Simple Expressions

We want to see:

  1. Whether we can write mathematical expressions the way we are used to from school. Ugh!
    1 + 2 * 3 ^ 4 should be treated as 1 + (2 * (3 ^ 4)), yielding 163. It should not be treated as ((1 + 2) * 3) ^ 4, which equals 6561,
  2. How many bits are necessary to store 10^6, and
  3. How steep is our driveway? (measured in degrees) Our garage is 7 meters away from the street and half a meter above it.

Here we go.

    cspiel@orion:~/articles/numerics $ octave

All three programs are console-based. That is, the user gets a prompt whenever the application is ready to accept input. We enter our first question as we write it on paper. Hitting return terminates the line, the program evaluates it, and returns the result in variable ans (more on variables later).

    octave:1> 1 + 2 * 3 ^ 4
    ans = 163

Aha, obviously GNU/Octave knows elementary-school math!

Our second question requires the logarithm function log, which returns the natural logarithm of its argument; this is, the logarithm to base e.

    octave:2> log(10^6) / log(2)
    ans = 19.932

We conclude that 1,000,000 needs 20 bits to be stored.

Finally, how steep is our driveway? What we need here is an angular function, namely the arctangent, written as atan(argument).

    octave:3> atan(0.50 / 7.0)
    ans = 0.071307

Hmm, ain't that a bit too flat? Digging out the wisdom of long forgotten math classes, we remember that the arctangent of 1 is 45 degrees. Let us check this!

    octave:4> atan(1)
    ans = 0.78540

Ouch, off by a factor of 57! Do we have to throw the program away? Wait -- 57 equals almost 180 over pi. This means GNU/Octave has returned the result in radians, not in degrees. All angular functions work in units of radians, this is, an angle of 360 degrees is equivalent 2 pi radians.

We try again, supplying the correct conversion factor:

    octave:5> atan(0.50 / 7.0) * 360/(2 * 3.14)
    ans = 4.0856

Approximately 4 degrees, that looks good. Our garage certainly won't get flooded in the next deluge.

Details

Differences

Variables

In the last section we have not gained much in comparison with a pocket calculator, have we? The first feature where our programs beat pocket calculators and spread-sheets are names that we can give parameters or results; these are called variables.

Assume our better half wants us to build a garden in the yard, but we want to watch basketball. Therefore we quickly need a hard figure that proves we don't have enough compost for the desired size. Ha -- brilliant idea!

[image: plan for a flower bed]

From our little plan we take the following lengths in feet:

    houseside_length = 10
    creekside_length = 6
    width = 2

Our better half also said the layer of new soil ought to be at least five inches, so

    height = 5 / 12

GNU/Octave to the rescue!

    octave:1> houseside_length = 10
    houseside_length = 10
    octave:2> creekside_length = 6
    creekside_length = 6
    octave:3> width = 2
    width = 2
    octave:4> height = 5 / 12
    height = 0.41667
    octave:5> volume = (houseside_length + creekside_length) * width * height
    volume = 13.333

The compost box is 6' x 4' and currently holds eight inches of usable compost.

    octave:6> box_width = 6
    box_wight = 6
    octave:7> box_depth = 4
    box_depth = 4
    octave:8> compost_height = 8/12
    compost_height = 0.66667
    octave:9> compost_volume = box_width * box_depth * compost_height
    compost_volume = 16

Oh no, we have just dug our own grave. We have got enough compost! What about taping the match on the VCR?

Details

Structured Data

Until now we have not exploited where computers are really good at: repetitive work.

Vectors

Say we got a long receipt from the grocery store. [Your ad here!] How can we get the VAT in Dollars on each item given the gross amount and the VAT rate in percent? The formula

            vat_percent / 100
    vat = --------------------- * gross_amount
          1 + vat_percent / 100

is trivial, but we want to save us repeated typing.

The list of all gross amounts in the receipt forms what numerical programs call a vector. Vectors are built from values by enclosing these values in square brackets and separating them with commas like this:

    octave:1> gross = [1.49, 4.98, 0.79, 5.49, 0.96, 0.96, 0.96, 0.96]
    gross =
      1.49000  4.98000  0.79000  5.49000  0.96000  0.96000  0.96000  0.96000

The vector is built from left to right using our supplied numbers in the same order that we enter them.

Wouldn't it be wonderful if we simply wrote: gross * (vat_percent/100) / (1 + vat_percent/100) and get the VAT of each item? It really is that simple.

    octave:2> vat_percent = 7
    vat_percent = 7
    octave:3> a = (vat_percent/100) / (1 + vat_percent/100)
    a = 0.065421
    octave:4> vat = a * gross
    vat =
      0.097477  0.325794  0.051682  0.359159  0.062804  0.062804  0.062804  0.062804

Wow -- it works! For the first time we have really gained convenience and expressiveness: a single multiplication sign performs eight multiplications in a row.

What has happened? vat_percent is a single value, which is called scalar in numerics to distinguish it from vectors. Well, if vat_percent is a scalar, then vat_percent/100, 1 + vat_percent/100, and a are scalars, too. Finally, scalar a must be multiplied with vector  gross. What we wanted and what happened was that a was multiplied in turn with every element of gross. This holds for every operator, not only multiplication! In general

     vector  op  scalar

and

     scalar  op  vector

apply scalar to every element in vector according to operation op. In our example, this is as if we had written the following

    vat(1) = a * gross(1)
    vat(2) = a * gross(2)
    ...
    vat(8) = a * gross(8)

where we have introduced a new piece of syntax: vector indexing. Each element (a scalar) of a vector can be accessed by its index, which is the number of it's place in the vector. The index is written in parenthesis after the vector. For example, to get the second element in gross, we write

    octave:5> gross(2)
    ans = 4.9800

Elements in vectors can be assigned to with the same syntax. Just place the indexed vector to the left of the assignment sign, for example, gross(2) = 5.12.

What else can be thought of as a vector of numbers besides our receipt? Any series of values! Most of the time the values will be related, like the temperature measured each day at 8am, the diameters of a batch of metal rods, the velocities of all westbound traffic across Buffalo Street at the corner of West Clinton Street on Wednesday April 18, 2001. As we are living in the Digital Age, many more series of data fit the notion of a vector: every piece of music on a CD is a vector of sound amplitudes and the indices mark discrete moments in time.

Details

Differences

Next Month

Christoph Spiel

Chris runs an Open Source Software consulting company in Upper Bavaria/Germany. Despite being trained as a physicist -- he holds a PhD in physics from Munich University of Technology -- his main interests revolve around numerics, heterogenous programming environments, and software engineering. He can be reached at cspiel@hammersmith-consulting.com.


Copyright © 2001, Christoph Spiel.
Copying license http://www.linuxgazette.net/copying.html
Published in Issue 69 of Linux Gazette, August 2001

[ Prev ][ Table of Contents ][ Front Page ][ Talkback ][ FAQ ][ Next ]