SCILAB Introduction for APAC-ANU Teaching Module
Steve Roberts & Graeme Chandler SCILAB Numerical Methods Tutorials A pdf version of this tutorial is available

1  Introduction

These exercises are intended to introduce you to the computer environment SCILAB. By the end of this tutorial you should have checked that your computer account is working and you will have practiced using an editor and SCILAB. One way or another you should try to get though all of this tutorial before the next tutorial starts, as it will require some programming work, so do not hesitate to ask for help. In this session you will mostly gain some experience with the SCILAB package. We will consider some simple SCILAB expressions and use them to consider some problems related to the use of floating point arithmetic. Along the way, questions will be asked. You should provide answers to these questions in your lab book. Marks for the tutorial will be awarded based on the work presented in your lab book.

2  Login

The startup procedure for SCILAB is a little different depending on what type of machine you are on and how it has been setup. But essentially you need to startup a SCILAB session, as well as an editor. For small calculations you can just use the SCILAB interactive command line, i.e. just type commands straight into SCILAB. For anything more complicated it is much more efficient to type your commands into a file (using your favorite text editor), save the file and then execute the commands from the editor (or using the exec filename command from the SCILAB command line). I will go through the procedure in a little more detail for Windows and Mac machines. So choose your system:

3  Windows Login

The integration of SCILAB into Windows is quite good. For this operating system SCILAB has an integrated editor. SCILAB is simply run from the application menu off the start menu. At the ANU, SCILAB can be found under Mathematics Software menu item. Once SCILAB is running you are met with the command window. Startup the editor (from the menu bar) and you are ready to start playing with SCILAB, either from the command line or from the editor.

4  Mac Login

Unfortunately the integration of SCILAB into the Mac world is still quite primitive. An editor is not integrated into SCILAB. You must choose your own. And SCILAB itself needs the Xserver running. At ANU we have setup a text editor, "text commander" in the application folder to help with the integration. If you use this application, you can edit files, and also startup a SCILAB session. This is not as good as the windows environment, but it is a start. If you don't have access to the ANU system, the procedure for starting SCILAB is a little more complicated. You will first need to ensure that an Xserver is running. This will depend on how you have setup your system. Once you have your Xserver running you will be able to start up an xterm. From the xterm command line, type scilab&. This will startup SCILAB. Startup your favorite editor and you are ready to go.

5  Talking between SCILAB and the Editor

From the SCILAB command window you can run SCILAB commands. In particular the command pwd will tell you the directory you are present working in. This is were SCILAB will look for files to execute. You can change this using CHANGE DIRECTORY item under the file menu or from the command line using the command chdir directory. From your text editor create some SCILAB commands such as

A = [ 1 2 3 ; 4 5 6 ; 7 8 9]

and then save the file (say into file test-01.sce) into the same directory in which SCILAB is looking (as shown by the pwd command). Now we can run the commands in the file via the command at the SCILAB prompt

exec test-01.sce

As you make changes to your file, you will need to save those changes and then run the exec command again (you can save typing by using the up arrow to recall previous commands). If you are using the editor integrated into SCILAB on the windows system, then the editor automatically asks to save the file you are working on, and automatically executes the file.

6  Basic Commands

Now let us play with SCILAB. Move the pointer to the SCILAB window and select the window (press the left mouse button). Now you should be able to type commands into the SCILAB window. SCILAB is made to easily work with vectors and matrices. First input a vector

 x = [ 0 ; 2 ; 5]

and then a matrix

A = ⎡ A2; A2; A2; ⎣ 1 3 4 −1 2 5 4 −3 5 ⎤ ⎥ ⎥ ⎥ ⎥

via the command

 A = [ 1  3  4 ; -1  2  5 ; 4  -3  5 ]

The semi-colons may be replaced by returns. SCILAB is case sensitive so that for example a and A will represent different variables. All the built-in functions in SCILAB like the rank or the eigenvalues for a matrix, are denoted by lower-case names, as in rank(A) or spec(A). If you cannot remember how to use a particular SCILAB command try the help command in SCILAB. Try

 help spec

SCILAB provides many vector and matrix operations. For instance we can multiply matrices and vectors. Try the simple commands

 A*x
 A*A

Lab Book: Record the results of these calculations.
SCILAB does all its calculation in double precision but formats the output using a short format. The format can be changed to show full precision with the SCILAB command

 format(20,'e')

Try some calculations to see the different format. The format can be changed back using

 format(10,'v')

A semi-colon at the end of a line will stop SCILAB output for that line. This is useful if you are working with large vectors or matrices. Try

 y = A*x ;

Note that the answer is not printed out. To see y just type

 y

Let's make some larger vectors. It is easy to produce vectors that have components which increment nicely. Try

 w = 0:0.1:7

This will produce a vector starting at 0 and increasing by 0.1 until 7 is reached. Note that a semi-colon at the end of the line would be helpful here. We can take functions of vectors (or matrices). For instance try

 z = sin(w)

Now we have two vectors w and z of the same length. We can plot the values of w versus the values of z. Try

 plot2d(w,z)

A plot of w versus z should be produced. You can save the graph to a file by using the file - export menu item on the graphics window.
Lab Book: Print this plot and paste into your lab book

7  Matrix Vector problems

We can also solve matrix problems. Try

 y = A \ x

The vector y should solve the linear equation x = A*y (check this). The inverse of a matrix can also be calculated using the inv command. Use the inv command to solve the matrix equation A y = x.
Lab Book: Compare the results of A\x and inv(A)*x. Are they equal? Why?
The SCILAB command rand(n,m) produces a random n ×m matrix. Redefine A to be a random 5 ×5 matrix and x a random 5 ×1 vector. The \ command is more general and can be used to solve over-determined systems (systems with more equations than unknowns) by finding a "least squares solution". In fact under-determined systems can also be accommodated.
Lab Book: Try solving an under-determined system (Matrix has less rows than columns). Try an over-determined system (more rows than columns). Document your experiments.

8  Loops and Conditionals

SCILAB provides for, while and if statements to control flow within a program. Consider the example program for calculating e^x.

ans = 0; n = 1; term = 1;
while( ans + term ~= ans )
        ans = ans + term;
        term = term*x/n;
        n = n + 1;
end
ans

Here ~= means "not equal to". We can avoid retyping by saving the program in a file. Open up a text editor. Type in the program. Save to a file using the file menu. You will be prompted for a file name. Use the name ex.sci, the .sci is a standard extension for SCILAB function files. Go back to the SCILAB window and test your program ex. Type in a value for x and then issue the command exec('ex.sci'). For instance try

 x = 1.0
 exec('ex.sci')

to let SCILAB run your program with x = 1.0.
Lab Book: Using your program, what is the approximation to e¹, e¹⁰ and e⁻¹⁰.
We can make our program into a function. Return to the editing window and change the program as follows:

function y = ex(x)
// EX A simple function to calculate exp(x)
y = 0; n = 1; term = 1;
while( y + term ~= y )
        y = y + term;
        term = term*x/n;
        n = n + 1;
end
endfunction

The lines starting with // are comment lines. Save your changes. Now you can use the command ex(x). But you must load in your changes. Try

 exec('ex.sci')
 ex(1.0)

Alternatively if you are using the editor supplied with SCILAB you can just use the menu item "execute" to load in your file. This algorithm is useless for large negative values of x. On the other hand, as in the lectures we can use this algorithm as a basis for a more robust algorithm. First though we need to become familiar with the if statement. The SCILAB if statement has the simple form

if expression then
   Scilab statements
else
   Scilab statements
end

By the way, SCILAB has a for statement of the form

for i=v
  Scilab statements
end

Here i is a dummy variable to be used in the loop, and v is a vector, usually a range of numbers defined 1:20, that the dummy variable will iterate through. So

for j=-4:2:6
  disp(j**2)
end

will print out the squares of the even numbers from -4 to 6. Use your function ex and the if statement to produce a new function newex which produces reasonable approximations of e^x for x positive or negative (as demonstrated in the lectures). You can use the same file to define both the ex and the newex functions, but you will have to use exec again to load in the new function newex.
Lab Book: Produce a table of the results of the ex, and newex functions compared to the inbuilt exp SCILAB function. Do the table for a range of large negative through to large positive numbers.

9  The Help Command

The help command is the easiest way to find out more about specific Scilab commands. If you have forgotten small details, for example. The command help name gives information about the Scilab command name.

help sin     // Information about sin.
help +       // Gives links to help on Scilab operator names
help log     // This is enough information about log
             //     to show  log  means log to the base  e .

Note that often the help information provides an example of how the command is used. This can be particularly useful when experimenting with a new command. You can use cutting and pasting to run the example commands. In Scilab the apropos command can also be used to search for relevant information.

apropos logarithm  // Provides a list of
                   // functions related to logarithms

10  The Help Window

It is also possible to get help by clicking on the HELP menu above the command window. From the HELP menu, select the HELP DIALOG item and the list of help topics is displayed. The advantage of this method is that it is possible to navigate around different topics and zero in on useful commands. For example, in the help dialog window click on the chapter 'Linear Algebra'. Then find the item 'linsolve'. Once found you click on 'show' to see the relevant help page. Note that double clicking doesn't work in windows Scilab version 2.6 and there is a slight bug when you change chapters and scroll to new items. Choosing the item a second time seems to work. Generally the help window is a good way to explore Scilab's commands.

Exercises

   

1. Is the inverse sine function one of Scilab's elementary functions? (The inverse sine is also known as sin⁻¹, arcsin, or asin.)
   1. How do you find sin⁻¹(.5) in Scilab? Use help to find out.
   2. If x=.5, is sin(sin⁻¹(x))−x exactly zero in Scilab?
   3. If x = π/3, is sin⁻¹(sin(x ))−x exactly zeros in Scilab? What about x = 5π/11?
2. Does SCILAB have a function to convert numbers to base 16, i.e. to hexadecimal form? (Hint: Use apropos to find a way to convert a decimal number to hexadecimal.) What is 61453 in base 16? Computers almost always represent numbers internally as hexadecimals.
3. Look for information about logarithms. Note that searching for logarithms fails where as logarithm succeeds. There are 8 entries, including the command, logm, for calculating the log of a matrix!
4. To wind down, use the SCILAB menu command DEMOS. This brings a menu of demonstrations and examples to explore.
   1. Visit the graphics to see a number of attractive images.
   2. I also like the car and trailer parking demo.

Lab Book: Record your results in your lab book.

11  Summary

We have tried to give a very quick idea of the facilities available in SCILAB and to the give an idea of using the SCILAB program in conjunction with a text editor to efficiently experiment with matrix problems. Also the demos available from the menu gives some good ideas. But now see if you can tackle the following problem.

12  Exercises

Exercises

Write a SCILAB program quadroots to compute and print the roots of a quadratic equation ax² + bx + c = 0 using the quadratic formula

x = −b ± √ b2 − 4ac 2a .

It should run with a command like

quadroots(1,3,2)

for the case of x² + 3x + 2 = 0. This program can be made to return both roots in a vector by calling them, for example, root1 and root2 in the function file and starting the function definition with the line

function [root1, root2] = quadroots(a,b,c)

These can then be stored into two variables such as r1 and r2 by running the function with the command

[r1 r2] = quadroots(1,3,2)

The program should avoid division by zero, printing an appropriate message in that case. Make your own decision about what to do in the case of complex roots, since SCILAB can compute and print complex values. Test quadroots on the following examples, checking the results in each case, and pay particular attention to the last example.

1. x² + 1 = 0
2. 0x² + 2x + 1 = 0
3. x² + 3x + 2 = 0
4. 4x² + 24x + 36 = 0
5. 10¹⁸ x² − x − 1 = 0
6. 10⁻¹⁸ x² − x + 1 = 0 (The roots are approximately 1+ 10⁻¹⁸ and 10¹⁸−1), which you can check by hand.)

The last example above illustrates the problem of catastrophic cancellation. Explain why it makes one root highly inaccurate but not the other, and devise a better algorithm for computing both roots accurately. You can use the fact that the roots x₁ and x₂ of ax² + bx+ c = 0 always satisfy x₁ x₂ = c/a. Write a modified program quadroots2 using this improved method, and test it on the same cases as above.

Unit Testing

It is very useful to write a test program which will test a program against known results. Here is a program test_quadroots together with some auxiliary functions which takes a quadroots function (quadroots or quadroots2) and tests the program with the examples given above. This is known as unit testing. Copy the following code into your SCILAB environment, and test your quadroots functions with the command test_quadroots(quadroots).

//--------------------------------------------------
//
// Test program for the quadroots function
//
// Test a range of inputs.
// Look at the expected results in the table
// to see how to code your function
//--------------------------------------------------
function test_quadroots(quadroots)


n=0
table = zeros(1,5);

//-------------------------------------
// Setup tests, a,b,c r1,r2 in a table
// r1 and r2 are expected roots. We are using
// %nan for non real roots
//-------------------------------------
// x^2 + 1 = 0
n=n+1; table(n,:) = [ 1.0, 0.0, 1.0 %nan, %nan]

// 0x^2 + 2x + 1 = 0
n=n+1; table(n,:) = [ 0.0, 2.0, 1.0, -0.5, %nan]

// x^2 + 3x + 2 = 0
n=n+1; table(n,:) = [ 1.0, 3.0, 2.0, -1.0, -2.0]

// 4x^2 + 24x + 36 = 0
n=n+1; table(n,:) = [ 4.0, 24.0, 36.0, -3.0, -3.0]

//10^18 x^2 - x - 1 = 0
n=n+1; table(n,:) = [ 1.0e18, -1.0, -1.0, 1.0000000005e-9, -9.999999995e-10]

//10^-18 x^2 - x + 1 = 0 (The roots are approximately 1 + 10e-18 and 10^18 - 1))
n=n+1; table(n,:) = [ 1.0e-18, -1.0, 1.0, 1.0e18, 1.0]


//--------------------------------------
// Test each problem.
//--------------------------------------
n = size(table,1)
success = 0.0
for i=1:n,
  success = success + ...
     print_test(table(i,1),table(i,2),table(i,3),table(i,4),table(i,5),quadroots)
end
printf('================================================\n')
printf('%g successful tests out of %g\n',success,n)

endfunction

//---------------------------------------------
// Print whether specific quadtest is correct.
// Returns 1.0 if correct, 0.0 otherwise
//---------------------------------------------
function success = print_test(a,b,c,r1,r2,quadroot)
success = 0.0

printf('================================================\n')
printf('TEST: %g x^2 + %g x + %g = 0\n',a,b,c)
printf('Expect   r1 = %20.15e, r2 = %20.15e \n',r1,r2)
[c1,c2]= quadroots(a,b,c)
printf('Obtained r1 = %20.15e, r2 = %20.15e\n',c1,c2)

if( (isnan(r1) & isnan(c1)) & (isnan(r2) & isnan(c1))) then
    printf('TEST PASSED\n')
    success = 1.0
    return
end

if( (isnan(r2) & isnan(c2)) & close(r1,c1) ) then
    printf('TEST PASSED\n')
    success = 1.0
    return
end

if (close(r1,c1)&close(r2,c2)) then
    printf('TEST PASSED\n')
    success = 1.0
    return
end

if (close(r1,c2)&close(r2,c1)) then
    printf('TEST PASSED\n')
    success = 1.0
    return
end

success = 0.0
printf('TEST FAILED\n')

endfunction

//------------------------------------
// Check whether two values are relatively
// close.
//------------------------------------
function r = close(x,y)

xx = max(abs(x),abs(y))

if xx ~= 0.0 then
    r = (abs(x-y)/xx)<%eps*2
else
    r = abs(x-y)<%eps*2
end
return


endfunction

Evolve your quadroots program to a stage in which it successfully passes all of the test_quadroots tests.

Lab Book: Add to the test_quadroots function a test for the quadratic 10⁻¹⁸ x² − x + 1 = 0 which should have roots approximately 1 + 10⁻¹⁸ and 10¹⁸−1. Copy your commented final code, with the functions quadroots, quadroots2 and your amended test_quadroots function into your lab book. Present the results of the test_quadroots function showing the success rate of your various implementations of quadroots.

---

Next Tutorial

---