1168 lines
35 KiB
HTML
1168 lines
35 KiB
HTML
<!--startcut ==============================================-->
|
|
<!-- *** BEGIN HTML header *** -->
|
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
|
|
<HTML><HEAD>
|
|
<title>Numerical Workbenches, part II LG #70</title>
|
|
</HEAD>
|
|
<BODY BGCOLOR="#FFFFFF" TEXT="#000000" LINK="#0000FF" VLINK="#0000AF"
|
|
ALINK="#FF0000">
|
|
<!-- *** END HTML header *** -->
|
|
|
|
<CENTER>
|
|
<A HREF="http://www.linuxgazette.com/">
|
|
<IMG ALT="LINUX GAZETTE" SRC="../gx/lglogo.png"
|
|
WIDTH="600" HEIGHT="124" border="0"></A>
|
|
<BR>
|
|
|
|
<!-- *** BEGIN navbar *** -->
|
|
<IMG ALT="" SRC="../gx/navbar/left.jpg" WIDTH="14" HEIGHT="45" BORDER="0" ALIGN="bottom"><A HREF="mcgucken.html"><IMG ALT="[ Prev ]" SRC="../gx/navbar/prev.jpg" WIDTH="16" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="index.html"><IMG ALT="[ Table of Contents ]" SRC="../gx/navbar/toc.jpg" WIDTH="220" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><A HREF="../index.html"><IMG ALT="[ Front Page ]" SRC="../gx/navbar/frontpage.jpg" WIDTH="137" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="http://www.linuxgazette.com/cgi-bin/talkback/all.py?site=LG&article=http://www.linuxgazette.com/issue70/spiel.html"><IMG ALT="[ Talkback ]" SRC="../gx/navbar/talkback.jpg" WIDTH="121" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><A HREF="../faq/index.html"><IMG ALT="[ FAQ ]" SRC="./../gx/navbar/faq.jpg"WIDTH="62" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="tranter.html"><IMG ALT="[ Next ]" SRC="../gx/navbar/next.jpg" WIDTH="15" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><IMG ALT="" SRC="../gx/navbar/right.jpg" WIDTH="15" HEIGHT="45" ALIGN="bottom">
|
|
<!-- *** END navbar *** -->
|
|
<P>
|
|
</CENTER>
|
|
|
|
<!--endcut ============================================================-->
|
|
|
|
<H4 ALIGN="center">
|
|
"Linux Gazette...<I>making Linux just a little more fun!</I>"
|
|
</H4>
|
|
|
|
<P> <HR> <P>
|
|
<!--===================================================================-->
|
|
|
|
<center>
|
|
<H1><font color="maroon">Numerical Workbenches, part II</font></H1>
|
|
<H4>By <a href="mailto:cspiel@hammersmith-consulting.com">Christoph Spiel</a></H4>
|
|
</center>
|
|
<P> <HR> <P>
|
|
|
|
<!-- END header -->
|
|
|
|
|
|
|
|
|
|
<p>In
|
|
<A HREF="../issue69/spiel.html">Part I</A> ,
|
|
we looked at the most basic operations of the numerical
|
|
workbenches GNU/Octave 2.1.34, Scilab 2.6, and Tela 1.32. This
|
|
time we will talk about matrices, have a look at some of the predefined
|
|
functions, learn how to write our own functions, and introduce flow control
|
|
statements. The article closes with a brief discussion of the applications'
|
|
input and output facilities.</p>
|
|
|
|
<h2><a name="matrices">Matrices</a></h2>
|
|
|
|
<p>Vectors help a lot if data depend on a single parameter. The different
|
|
parameter values are reflected by different index values. If data depend on
|
|
two parameters, vectors are a clumsy container and a more general structure,
|
|
which allows for two independent indices is needed. Such a structure is called
|
|
a matrix. Matrices are packed like a fresh six-pack: they are rectangular
|
|
storage containers and no bottle -- oops -- element is missing.</p>
|
|
|
|
<p>Matrices are, for example, built from scalars as the next transcript of a
|
|
GNU/Octave session demonstrates.</p>
|
|
|
|
<pre>
|
|
octave:1> # temperature rain sunshine
|
|
octave:1> # degF inches hours
|
|
octave:1> weather_data = [ 73.4, 0.0, 10.8; ...
|
|
> 70.7, 0.0, 8.5; ...
|
|
> 65.2, 1.3, 0.7; ...
|
|
> 68.2, 0.2, 4.1]
|
|
weather_data =
|
|
</pre>
|
|
|
|
<pre>
|
|
73.40000 0.00000 10.80000
|
|
70.70000 0.00000 8.50000
|
|
65.20000 1.30000 0.70000
|
|
68.20000 0.20000 4.10000
|
|
</pre>
|
|
|
|
<p>Three new ideas appear in the example. First, we have introduced some
|
|
comments to label the columns of our matrix. A comment starts with a pound
|
|
sign ``<code>#</code>'' and extends until the end of the line. Second,
|
|
the rows of a matrix are separated by semi-colons ``<code>;</code>'', and
|
|
third, if an expression stretches across two or more lines, the unfinished
|
|
lines must end with the line-continuation operator ``<code>...</code>''.</p>
|
|
|
|
<p>Similarly to vectors, matrices can not only be constructed from scalars, but
|
|
from vectors or other matrices. If we had some variables holding the weather
|
|
data of each day, like</p>
|
|
|
|
<pre>
|
|
weather_mon = [73.4, 0.0, 10.8]
|
|
weather_tue = [70.7, 0.0, 8.5]
|
|
weather_wed = [65.2, 1.3, 0.7]
|
|
weather_thu = [68.2, 0.2, 4.1]
|
|
</pre>
|
|
|
|
<p>we would have defined <code>weather_data</code> with</p>
|
|
|
|
<pre>
|
|
weather_data = [weather_mon; weather_tue; weather_wed; weather_thu]
|
|
</pre>
|
|
|
|
<p>or, on the other hand, if we had the data from the various instruments
|
|
as</p>
|
|
|
|
<pre>
|
|
temperature = [73.4; 70.7; 65.2; 68.2]
|
|
rain = [0.0; 0.0; 1.3; 0.2]
|
|
sunshine = [10.8; 8.5; 0.7; 4.1]
|
|
</pre>
|
|
|
|
<p>we would have defined <code>weather_data</code> with</p>
|
|
|
|
<pre>
|
|
weather_data = [temperature, rain, sunshine]
|
|
</pre>
|
|
|
|
<p>The fundamental rule is: <em>Commas separate columns, semi-colons separate
|
|
rows.</em></p>
|
|
|
|
<p>The scalars living in matrix <code>m</code> are accessed by applying
|
|
two indices: <code>m(row, column)</code>, where <em>row</em> is the row-index,
|
|
and <em>column</em> is the column index. Thus, the amount of rain fallen on
|
|
Wednesday is fetched with the expression</p>
|
|
|
|
<pre>
|
|
octave:10> weather_data(3, 2)
|
|
ans = 1.3000
|
|
</pre>
|
|
|
|
<p>Entries are changed by assigning to them:</p>
|
|
|
|
<pre>
|
|
octave:11> weather_data(3, 2) = 1.1
|
|
weather_data =
|
|
</pre>
|
|
|
|
<pre>
|
|
73.40000 0.00000 10.80000
|
|
70.70000 0.00000 8.50000
|
|
65.20000 1.10000 0.70000
|
|
68.20000 0.20000 4.10000
|
|
</pre>
|
|
|
|
<p>Now that we have defined <code>weather_data</code> we want to work with it.
|
|
We can apply all binary operations that we have seen in last month's article
|
|
on vectors. However, for this particular example, computing</p>
|
|
|
|
<pre>
|
|
rain_forest_weather_data = weather_data + 2.1
|
|
siberian_summer_weather_data = weather_data / 3.8
|
|
</pre>
|
|
|
|
<p>does not make much sense, though the computer will not complain at all. In
|
|
the first example it would dutifully add <code>2.1</code> to every element of
|
|
<code>weather_data</code>, in the second it would -- obedient like a sheepdog
|
|
-- divide each element by <code>3.8</code>.</p>
|
|
|
|
<p>Say we want to do something meaningful to <code>weather_data</code> and
|
|
convert all temperatures from degrees Fahrenheit to degrees Celsius. To that
|
|
end, we need to access all elements in the first column. The vector of
|
|
interest is</p>
|
|
|
|
<pre>
|
|
octave:16> temp = [weather_data(1, 1); ...
|
|
> weather_data(2, 1); ...
|
|
> weather_data(3, 1); ...
|
|
> weather_data(4, 1)]
|
|
temp =
|
|
</pre>
|
|
|
|
<pre>
|
|
73.400
|
|
70.700
|
|
65.200
|
|
68.200
|
|
</pre>
|
|
|
|
<p>Obviously, the row-indices <code>[1, 2, 3, 4]</code>
|
|
form a vector themselves. We can use a shortcut and write this vector of
|
|
indices as</p>
|
|
|
|
<pre>
|
|
temp = weather_data([1, 2, 3, 4], 1)
|
|
</pre>
|
|
|
|
<p>In general, any vector may be used as index vector. Just watch out that no
|
|
index is out of range. Ordering of the indices does matter (for example <code>
|
|
weather_data([2, 1, 4, 3], 1)</code> puts Tuesday's temperature in front) and
|
|
repeated indices are permitted (for example <code>weather_data([3, 3, 3, 3, 3,
|
|
3, 3], 1)</code> holds Wednesday's temperature seven times).</p>
|
|
|
|
<p>In our example, the index-vector can be generated by a special built-in,
|
|
the range generation operator ``<code>:</code>''. To make a vector that
|
|
starts at <em>low</em> and contains all integers from <em>low</em> to <em>
|
|
high</em>, we say</p>
|
|
|
|
<pre>
|
|
low:high
|
|
</pre>
|
|
|
|
<p>For example</p>
|
|
|
|
<pre>
|
|
octave:1> -5:2
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
-5 -4 -3 -2 -1 0 1 2
|
|
</pre>
|
|
|
|
<p>Our weather data example now simplifies to</p>
|
|
|
|
<pre>
|
|
temp = weather_data(1:4, 1)
|
|
</pre>
|
|
|
|
<p>Accessing a complete column or row is so common that further shortcuts
|
|
exist. If we drop both, <em>low</em> and <em>high</em> from the
|
|
colon-operator, it will generate all valid column indices for us. Therefore,
|
|
we reach at the shortest form to get all elements in the first column.</p>
|
|
|
|
<pre>
|
|
octave:17> temp = weather_data(:, 1)
|
|
temp =
|
|
</pre>
|
|
|
|
<pre>
|
|
73.400
|
|
70.700
|
|
65.200
|
|
68.200
|
|
</pre>
|
|
|
|
<p>With our new knowledge, we extract the sunshine hours on Tuesday,
|
|
Wednesday, and Thursday</p>
|
|
|
|
<pre>
|
|
octave:19> sunnyhours = weather_data(2:4, 3)
|
|
sunnyhours =
|
|
</pre>
|
|
|
|
<pre>
|
|
8.50000
|
|
0.70000
|
|
4.10000
|
|
</pre>
|
|
|
|
<p>and Tuesday's weather record</p>
|
|
|
|
<pre>
|
|
octave:20> tue_all = weather_data(2, :)
|
|
tue_all =
|
|
</pre>
|
|
|
|
<pre>
|
|
70.70000 0.00000 8.50000
|
|
</pre>
|
|
|
|
<p>Now it is trivial to convert the data on the rain from inches to
|
|
millimeters: Multiply the second column of <code>weather_data</code> by 25.4
|
|
(Millimeters per Inch) to get the amount of rain in metric units:</p>
|
|
|
|
<pre>
|
|
octave:21> rain_in_mm = 25.4 * weather_data(:, 2)
|
|
rain_in_mm =
|
|
</pre>
|
|
|
|
<pre>
|
|
0.00000
|
|
0.00000
|
|
27.94000
|
|
5.08000
|
|
</pre>
|
|
|
|
<p>We have already seen that vectors are compatible with scalars</p>
|
|
|
|
<pre>
|
|
1.25 + [0.5, 0.75, 1.0]
|
|
</pre>
|
|
|
|
<p>or</p>
|
|
|
|
<pre>
|
|
[-4.49, -4.32, 1.76] * 2
|
|
</pre>
|
|
|
|
<p>Scalars are also compatible with matrices.</p>
|
|
|
|
<pre>
|
|
octave:1> 1.25 + [ 0.5, 0.75, 1.0; ...
|
|
> -0.75, 0.5, 1.25; ...
|
|
> -1.0, -1.25, 0.5]
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
1.75000 2.00000 2.25000
|
|
0.50000 1.75000 2.50000
|
|
0.25000 0.00000 1.75000
|
|
</pre>
|
|
|
|
<pre>
|
|
octave:2> [-4.49, -4.32, 1.76; ...
|
|
> 9.17, 6.35, 3.27] * 2
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
-8.9800 -8.6400 3.5200
|
|
18.3400 12.7000 6.5400
|
|
</pre>
|
|
|
|
<p>In each case the result is the scalar applied to every element in the
|
|
vector or matrix.</p>
|
|
|
|
<p>How about vectors and matrices? Obviously, an expressions like</p>
|
|
|
|
<pre>
|
|
[7, 4, 9] + [3, 2, 7, 6, 6]
|
|
[2, 4; 1, 6] - [1, 1, 9, 4]
|
|
</pre>
|
|
|
|
<p>do not make any sense. In the first line the vectors disagree in size (3
|
|
vs. 5 elements), in the second line they have different shapes
|
|
(2 columns and 2 rows vs. 4 columns and 1 row). To make
|
|
sense, vectors or matrices that are used in an addition or subtraction must
|
|
have the same shape, which means the same number of rows and the same number
|
|
of columns. The technical term for ``shape'' in this context is dimension. We
|
|
can query the dimension of anything with the built-in function <code>
|
|
size()</code>.</p>
|
|
|
|
<pre>
|
|
octave:22> size(weather_data)
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
4 3
|
|
</pre>
|
|
|
|
<pre>
|
|
octave:23> size(sunnyhours)
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
3 1
|
|
</pre>
|
|
|
|
<p>The answer is a vector whose first element is the number of rows, and whose
|
|
second element is the number of columns of the argument.</p>
|
|
|
|
<p>Multiplication and division of matrices can be defined in two flavors, both
|
|
of which are implemented in the numerical workbenches.</p>
|
|
|
|
<ul>
|
|
<li>Element by element multiplication or division of two vectors or matrices
|
|
of same dimensions: The number in the first row and first column of the first
|
|
matrix is multiplied by the number in the first row and first column of the
|
|
second matrix and so on for every element.
|
|
|
|
<pre>
|
|
a = [3, 3; ...
|
|
6, 4; ...
|
|
6, 3]
|
|
b = [9, 3; ...
|
|
8, 2; ...
|
|
0, 3]
|
|
</pre>
|
|
|
|
<pre>
|
|
octave:1> a .* b
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
27 9
|
|
48 8
|
|
0 9
|
|
</pre>
|
|
|
|
<p>The element-by-element operators are preceded by a dot: element-by-element
|
|
multiplication ``<code>.*</code>'' and element-by-element division
|
|
``<code>./</code>''.</p>
|
|
</li>
|
|
|
|
<li>Matrix multiplication as known from Linear Algebra: <em>
|
|
c</em> = <em>a</em> * <em>b</em>, where <em>a</em> is a <em>
|
|
p</em>-times-<em>q</em> matrix and <em>b</em> is a <em>
|
|
q</em>-times-<em>r</em> matrix. The result <em>c</em> is a <em>
|
|
p</em>-times-<em>r</em> matrix.
|
|
|
|
<p>Example:</p>
|
|
|
|
<pre>
|
|
a = [3, 3; ...
|
|
6, 4; ...
|
|
6, 3]
|
|
</pre>
|
|
|
|
<pre>
|
|
b = [-4, 0, 1, -4; ...
|
|
-1, -3, 2, 0]
|
|
</pre>
|
|
|
|
<pre>
|
|
octave:1> a * b
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
-15 -9 9 -12
|
|
-28 -12 14 -24
|
|
-27 -9 12 -24
|
|
</pre>
|
|
|
|
<p>Although we have not seen <code>for</code>-loops yet (they will be discussed
|
|
<a href="#flow control statements">farther down</a>), I would like to write
|
|
the code behind the matrix multiplication operator ``<code>*</code>'' to
|
|
give the reader an impression of the operations involved.</p>
|
|
|
|
<pre>
|
|
for i = 1:p
|
|
for j = 1:r
|
|
sum = 0
|
|
for k = 1:q
|
|
sum = sum + a(i, k)*b(k, j)
|
|
end
|
|
c(i, j) = sum
|
|
end
|
|
end
|
|
</pre>
|
|
|
|
<p>Compare these triply nested <code>for</code>-loops with the simple
|
|
expression <code>c = a * b</code>.</p>
|
|
</li>
|
|
|
|
<li>Matrix division? You cannot divide by a matrix! However, operator
|
|
``<code>/</code>'' is defined for vectors and matrices. But writing <em>
|
|
x</em> = <em>b</em> / <em>a</em>, where <em>a</em> and <em>
|
|
b</em> are matrices or vectors has nothing to do with division at all! It
|
|
means: please solve the system of linear equations
|
|
|
|
<pre>
|
|
x * a = b
|
|
</pre>
|
|
|
|
<p>for <em>x</em>, given matrix <em>a</em> and the
|
|
right-hand-side(s) <em>b</em>. Here ``<code>*</code>'' denotes matrix
|
|
multiplication as defined in the previous item, and the same rules for
|
|
compatible dimensions of <em>a</em> and <em>b</em> apply.</p>
|
|
|
|
<pre>
|
|
a = [-2, 3, 1; ...
|
|
7, 8, 6; ...
|
|
2, 0, -1]
|
|
</pre>
|
|
|
|
<pre>
|
|
b = [-26, 5, -6; ...
|
|
24, 53, 26]
|
|
</pre>
|
|
|
|
<pre>
|
|
octave:1> x = b / a
|
|
x =
|
|
</pre>
|
|
|
|
<pre>
|
|
7.00000 -2.00000 1.00000
|
|
7.00000 4.00000 5.00000
|
|
</pre>
|
|
|
|
<p>Isn't that an easy way to solve a system of linear equations? Imagine you
|
|
had to write the code which does exactly that.</p>
|
|
|
|
<p>Finally, let us verify the result by multiplying with <em>a</em> again</p>
|
|
|
|
<pre>
|
|
octave:2> x*a
|
|
ans =
|
|
</pre>
|
|
|
|
<pre>
|
|
-26.0000 5.0000 -6.0000
|
|
24.0000 53.0000 26.0000
|
|
</pre>
|
|
|
|
<p>which, as expected, recovers <em>b</em>.</p>
|
|
</li>
|
|
</ul>
|
|
|
|
<p><strong>Details</strong></p>
|
|
|
|
<ul>
|
|
<li>For convenience GNU/Octave and Scilab define an alternative matrix
|
|
division operator ``<code>\</code>''. <em>x</em> = <em>
|
|
a</em> \ <em>b</em> solves the linear system of equations
|
|
|
|
<pre>
|
|
a * x = b
|
|
</pre>
|
|
|
|
<p>for <em>x</em>, given matrix <em>a</em> and the
|
|
right-hand-side(s) <em>b</em>. This is the form most users prefer,
|
|
because here <em>x</em> is a column vector, whereas
|
|
operator ``<code>/</code>'' returns <em>x</em> as row-vector.</p>
|
|
</li>
|
|
|
|
<li>operator ``<code>\</code>'' has the dotted cousin ``<code>.\</code>''
|
|
and the relation <em>a</em> ./ <em>b</em> == <em>
|
|
b</em> .\ <em>a</em> holds.</li>
|
|
</ul>
|
|
|
|
<p><strong>Differences</strong></p>
|
|
|
|
<ul>
|
|
<li>Scilab and Tela use C++-like comments
|
|
|
|
<pre>
|
|
// This is a Scilab or a Tela comment
|
|
</pre>
|
|
</li>
|
|
|
|
<li>Tela does not need or understand the line continuation operator
|
|
``<code>...</code>''
|
|
|
|
<pre>
|
|
weather_data = #(73.4, 0.0, 10.8;
|
|
70.7, 0.0, 8.5;
|
|
65.2, 1.3, 0.7;
|
|
68.2, 0.2, 4.1)
|
|
</pre>
|
|
|
|
<p>In interactive mode, Tela does not handle multi-line expressions as the
|
|
above. Multi-line expressions must be read from a file (with <code>
|
|
source("filename.t")</code>).</p>
|
|
</li>
|
|
|
|
<li>In Tela the operators ``<code>*</code>'' and ``<code>/</code>'' work
|
|
element by element, this is, they work like ``<code>.*</code>'' and
|
|
``<code>./</code>'' do in GNU/Octave and Scilab. Matrix multiplication
|
|
(<em>a</em> * <em>b</em> in GNU/Octave or Scilab) is written as
|
|
|
|
<pre>
|
|
a ** b
|
|
</pre>
|
|
|
|
<p>or</p>
|
|
|
|
<pre>
|
|
matmul(a, b)
|
|
</pre>
|
|
|
|
<p>solving systems of linear equations (<em>b</em> / <em>a</em> in
|
|
Octave or Scilab) as</p>
|
|
|
|
<pre>
|
|
linsolve(a, b)
|
|
</pre>
|
|
</li>
|
|
</ul>
|
|
|
|
<h2><a name="builtin matrix functions">Built-In Matrix Functions</a></h2>
|
|
|
|
<p>Ugh -- far too many to mention! The workbenches supply dozens of predefined
|
|
functions. Here I can only wet the reader's
|
|
appetite.</p>
|
|
|
|
<dl>
|
|
<dt><strong><a name="item_Generating_Special_Matrices">Generating Special
|
|
Matrices</a></strong><br>
|
|
</dt>
|
|
|
|
<dd>Several matrices occur often enough in computations that they have been
|
|
given their own generating functions. These are for example, <em>
|
|
m</em>-times-<em>n</em> matrices filled with zeros: <code>zeros(m, n)</code>
|
|
or ones: <code>ones(m, n)</code>, or <em>n</em>-times-<em>n</em> diagonal
|
|
matrices, where the diagonal consists entirely of ones: <code>eye(n)</code> or
|
|
the diagonal is set to numbers supplied in a vector: <code>diag([a1, a2, ...,
|
|
I<an>])</code>.</dd>
|
|
|
|
<dt><strong><a name="item_Analyzing_Matrices">Analyzing
|
|
Matrices</a></strong><br>
|
|
</dt>
|
|
|
|
<dd>Getting the smallest or largest element in matrix <em>a</em>: <code>
|
|
min(a)</code>, <code>max(a)</code>, or totaling matrix <em>a</em>: <code>
|
|
sum(a)</code>.
|
|
|
|
<p><strong>Differences:</strong> GNU/Octave's <code>min(a)</code>, <code>
|
|
max(a)</code>, and <code>sum(a)</code> return the column-wise result as a row
|
|
vector. To get the minimum, maximum, and sum of all elements in
|
|
matrix <em>a</em>, use <code>min(min(a))</code>, <code>
|
|
max(max(a))</code>, <code>sum(sum(a))</code>.</p>
|
|
</dd>
|
|
|
|
<dt><strong><a name="item_Linear_Algebra">Linear Algebra</a></strong><br>
|
|
</dt>
|
|
|
|
<dd>We mentioned that systems of linear equations, like <em>
|
|
x</em> * <em>a</em> = <em>b</em>, are solved for <em>x</em>
|
|
with the slash operator ``<code>/</code>''. But many more linear algebra
|
|
functions exist, for example singular value decomposition: <code>
|
|
svd(a)</code>, or eigenvalue computation: <code>eig(a)</code>.
|
|
|
|
<p><strong>Differences:</strong> In Tela uses <code>SVD(a)</code> instead of
|
|
<code>svd(a)</code>, and instead of <code>eig(a)</code>, Scilab uses <code>
|
|
spec(a)</code> to compute the eigenvalue spectrum.</p>
|
|
</dd>
|
|
</dl>
|
|
|
|
<p>One note on performance: basically, all three applications are
|
|
interpreters. This means that each expression is first parsed, then the
|
|
interpreter performs desired computations, finally calling the functions
|
|
inside of the expressions -- all in all a relatively slow process in
|
|
comparison to a compiled program. However, functions like those shown above
|
|
are used in their compiled form! They execute almost at top speed. What the
|
|
interpreter does in these cases is to hand over the complete matrix to a
|
|
compiled Fortran, C, or C++ function, let it do all the work, and then
|
|
pick up the result.</p>
|
|
|
|
<p>Thus we deduce one of the fundamental rules for successful work with
|
|
numerical workbenches: prefer compiled functions over interpreted code.
|
|
It makes a tremendous difference in execution speed.</p>
|
|
|
|
<h2><a name="user defined functions">User Defined Functions</a></h2>
|
|
|
|
<p>No matter how many functions a program may provide its users, they are never
|
|
enough. Users always need specialized functions to deal with their problems,
|
|
or they simply want to group repeated, yet predefined operations. In other
|
|
words, there always is a need for user-defined functions.</p>
|
|
|
|
<p>User functions are best defined in files, so that they can be used again in
|
|
later sessions. For GNU/Octave, functions files end in <em>.m</em>, and are
|
|
loaded either <a href="#automagical_explanation">automagically</a> or with
|
|
<code>source("<em>filename.m</em>")</code>. Scilab calls its
|
|
function files <em>.sci</em>, and requires them to be loaded with <code>
|
|
getf("<em>filename.sci</em>")</code>. Tela functions are stored
|
|
in <em>.t</em>-files and loaded with <code>
|
|
source("<em>filename.t</em>")</code>. As big as the differences
|
|
are in loading functions, all workbenches use quite similar syntax for the
|
|
definition of functions.</p>
|
|
|
|
<p>GNU/Octave and Scilab</p>
|
|
|
|
<pre>
|
|
function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN)
|
|
# function body
|
|
endfunction
|
|
</pre>
|
|
|
|
<p>Tela</p>
|
|
|
|
<pre>
|
|
function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN)
|
|
{
|
|
// function body
|
|
};
|
|
</pre>
|
|
|
|
<p>where <em>arg1</em> to <em>argN</em> are the functions' arguments (also
|
|
known as parameters), and <em>res1</em> to <em>resN</em> are the return
|
|
values. Yes, trust your eyes, multiple return values are permitted, what might
|
|
come as a surprise to most readers who are acquainted with popular programming
|
|
languages. However, this is a necessity, as no function is allowed to change
|
|
any of its input arguments.</p>
|
|
|
|
<p>Enough theory! let us write a function that takes a matrix as input and
|
|
returns a matrix of the same dimensions, with the entries rescaled to lie in
|
|
the interval (0, 1).</p>
|
|
|
|
<pre>
|
|
### Octave
|
|
</pre>
|
|
|
|
<pre>
|
|
function y = normalize(x)
|
|
## Return matrix X rescaled to the interval (0, 1).
|
|
</pre>
|
|
|
|
<pre>
|
|
minval = min(min(x))
|
|
maxval = max(max(x))
|
|
</pre>
|
|
|
|
<pre>
|
|
y = (x - minval) / (maxval - minval)
|
|
endfunction
|
|
</pre>
|
|
|
|
<p>Now define a Scilab function that returns the spectral radius on a matrix.
|
|
We use <code>abs()</code> which returns the magnitude of its (possibly
|
|
complex) argument.</p>
|
|
|
|
<pre>
|
|
// Scilab
|
|
</pre>
|
|
|
|
<pre>
|
|
function r = spectral_radius(m)
|
|
// Return the spectral radius R of matrix M.
|
|
</pre>
|
|
|
|
<pre>
|
|
r = max(abs(spec(m)))
|
|
endfunction
|
|
</pre>
|
|
|
|
<p>Finally, we write a Tela function which computes the Frobenius norm of a
|
|
matrix.</p>
|
|
|
|
<pre>
|
|
// Tela
|
|
</pre>
|
|
|
|
<pre>
|
|
function x = frobenius(m)
|
|
// Return the Frobenius norm X of matrix M.
|
|
{
|
|
x = sqrt(sum(abs(m)^2))
|
|
};
|
|
</pre>
|
|
|
|
<p><strong>Details:</strong></p>
|
|
|
|
<a name="automagical_explanation"></a>
|
|
<p>GNU/Octave's ``automagical'' function file loading works the following way:
|
|
if Octave runs into an undefined function name it searches the list of
|
|
directories specified by the built-in variable <code>LOADPATH</code> for files
|
|
ending in .m that have the same base name as the undefined function; for
|
|
example, <code>x = my_square_root(2.0)</code> looks for the file <em>
|
|
my_square_root.m</em> in the directories listed in <code>LOADPATH</code>.</p>
|
|
|
|
<h2><a name="flow control statements">Flow Control Statements</a></h2>
|
|
|
|
<p>All code we have written thus far executes strictly top-to-bottom, we have
|
|
not used any flow control statements such as conditionals or loops.</p>
|
|
|
|
<p>Before we manipulate the flow of control, we should look at logical
|
|
expressions because the conditions used in conditionals and loops depend on
|
|
them. Logical expressions are formed from (1.) numbers,
|
|
(2.) comparisons, and (3.) logical expressions catenated with
|
|
logical operators.</p>
|
|
|
|
<ol>
|
|
<li>Zero means logically false, any number not equal to zero means logically
|
|
true, hence C-programmers should feel at home.</li>
|
|
|
|
<li>The usual gang of comparison operators exist: less-than
|
|
``<code><</code>'', less-or-equal ``<code><=</code>'', greater-than
|
|
``<code>></code>'', greater-or-equal ``<code>>=</code>'', and equal
|
|
``<code>==</code>''.
|
|
|
|
<p><strong>Differences:</strong> The inequality operator varies quite a bit
|
|
among the programs. (Octave cannot decide whether it feels like C, Smalltalk,
|
|
or Pascal. Scilab wants to be Smalltalk and Pascal at the same time. :-)</p>
|
|
|
|
<pre>
|
|
!= ~= <> # Octave
|
|
~= <> // Scilab
|
|
!= // Tela
|
|
</pre>
|
|
</li>
|
|
|
|
<li>Complex logical expressions are formed with logical operators ``and'',
|
|
``or'' and ``not'' whose syntax is borrowed from C. However, each program uses
|
|
its own set of operators. Thus, we have to list some
|
|
|
|
<p><strong>Differences:</strong></p>
|
|
|
|
<pre>
|
|
and or not
|
|
---- ---- ----
|
|
& | ! ~ # Octave
|
|
& | ~ // Scilab
|
|
&& || ! // Tela
|
|
</pre>
|
|
</li>
|
|
</ol>
|
|
|
|
<p>We are all set now for the first conditional, the <code>
|
|
if</code>-statement. Note that the parenthesis around the conditions are
|
|
mandatory (as they are in C). The <code>else</code>-branches are optional in
|
|
any case.</p>
|
|
|
|
<pre>
|
|
# Octave // Scilab // Tela
|
|
</pre>
|
|
|
|
<pre>
|
|
if (cond) if cond then if (cond) {
|
|
# then-body // then-body // then-body
|
|
else else } else {
|
|
# else-body // else-body // else-body
|
|
endif end };
|
|
</pre>
|
|
|
|
<p><em>cond</em> is a logical expression as described above.</p>
|
|
|
|
<p><code>while</code>-statements:</p>
|
|
|
|
<pre>
|
|
# Octave // Scilab // Tela
|
|
</pre>
|
|
|
|
<pre>
|
|
while (cond) while cond while (cond) {
|
|
# body // body // body
|
|
endwhile end };
|
|
</pre>
|
|
|
|
<p>Again, <em>cond</em> is a logical expression.</p>
|
|
|
|
<p><code>for</code>-statements in Octave and Scilab walk through the columns
|
|
of <em>expr</em> one by one. Most often <em>expr</em> will be a vector
|
|
generated with the range operator ``<code>:</code>'', like <code>for i =
|
|
1:10</code>. Tela's <code>for</code>-statement is the same as C's.</p>
|
|
|
|
<pre>
|
|
# Octave // Scilab // Tela
|
|
</pre>
|
|
|
|
<pre>
|
|
for var = expr for var = expr for (init; cond; step) {
|
|
# body // body // body
|
|
endfor end };
|
|
</pre>
|
|
|
|
<p>Here come some examples which only show things we have discussed so
|
|
far.</p>
|
|
|
|
<p>Octave</p>
|
|
|
|
<pre>
|
|
function n = catch22(x0)
|
|
## The famous catch-22 function: it is
|
|
## impossible to compute that it will
|
|
## stop for a specific input. Returns
|
|
## the number of loops.
|
|
</pre>
|
|
|
|
<pre>
|
|
n = 0
|
|
x = x0
|
|
while (x != 1)
|
|
if (x - floor(x/2)*2 == 0)
|
|
x = x / 2
|
|
else
|
|
x = 3*x + 1
|
|
endif
|
|
n = n + 1
|
|
endwhile
|
|
endfunction
|
|
</pre>
|
|
|
|
<p>Scilab</p>
|
|
|
|
<pre>
|
|
function m = vandermonde(v)
|
|
// Return the Vandermonde matrix M based on
|
|
// vector V.
|
|
</pre>
|
|
|
|
<pre>
|
|
[rows, cols] = size(v)
|
|
m = [] // empty matrix
|
|
if rows < cols then
|
|
for i = 0 : (cols-1)
|
|
m = [m; v^i]
|
|
end
|
|
else
|
|
for i = 0 : (rows-1)
|
|
m = [m, v^i]
|
|
end
|
|
end
|
|
endfunction
|
|
</pre>
|
|
|
|
<p>Tela</p>
|
|
|
|
<pre>
|
|
function vp = sieve(n)
|
|
// Sieve of Erathostenes; returns vector of
|
|
// all primes VP that are strictly less than
|
|
// 2*N. 1 is not considered to be a prime
|
|
// number in sieve().
|
|
{
|
|
vp = #(); // empty vector
|
|
if (n <= 2) { return };
|
|
</pre>
|
|
|
|
<pre>
|
|
vp = #(2);
|
|
flags = ones(1, n + 1);
|
|
for (i = 0; i <= n - 2; i = i + 1)
|
|
{
|
|
if (flags[i + 1])
|
|
{
|
|
p = i + i + 3;
|
|
vp = #(vp, p);
|
|
for (j = p + i; j <= n; j = j + p)
|
|
{
|
|
flags[j + 1] = 0
|
|
}
|
|
}
|
|
}
|
|
};
|
|
</pre>
|
|
|
|
<h2><a name="input/output">Input/Output</a></h2>
|
|
|
|
<p>We have been using with the workbenches a lot. At some point we would like
|
|
to call it a day, but we do not want to lose all of our work. Our functions
|
|
are already stored in files. It is time to see how to make our data
|
|
persist.</p>
|
|
|
|
<h3><a name="simple input and output">Simple Input and Output</a></h3>
|
|
|
|
<p>All three applications at least have one input/output (I/O) model that
|
|
borrows heavily from the C programming language. This model allows close
|
|
control of the items read or written. Often though, it is unnecessary to take
|
|
direct control over the file format written. If variables must be saved just
|
|
to be restored later, simplified I/O commands will do.</p>
|
|
|
|
<ul>
|
|
<li>Octave offers the most flexible solution with the <code>
|
|
save</code>/<code>load</code> command pair.
|
|
|
|
<pre>
|
|
save filename varname1 varname2 ... varnameN
|
|
</pre>
|
|
|
|
<p>saves the variables named <em>varname1</em>, <em>varname2</em>, ..., <em>
|
|
varnameN</em> in file <em>filename</em>. The complementary</p>
|
|
|
|
<pre>
|
|
load filename varname1 varname2 ... varnameN
|
|
</pre>
|
|
|
|
<p>command restores them from <em>filename</em>. If <code>load</code> is given
|
|
no variable names, all variables form <em>filename</em> are loaded. Handing
|
|
over names to <code>load</code> selects only the named variables for
|
|
loading.</p>
|
|
|
|
<p>Note that the <code>save</code> and <code>load</code> commands do not have
|
|
parenthesis and their arguments are separated by spaces not commas. Filename
|
|
and variable names are strings.</p>
|
|
|
|
<pre>
|
|
save "model.oct-data" "prantl" "reynolds" "grashoff"
|
|
load "model.oct-data" "reynolds"
|
|
</pre>
|
|
|
|
<p>By default <code>load</code> does not overwrite existing variables, but
|
|
complain with an error if the user tries to do so. When it is save to discard
|
|
of the values of existing variables, add option ``<code>-force</code>''
|
|
to <code>load</code>, like</p>
|
|
|
|
<pre>
|
|
load -force "model.oct-data" "reynolds"
|
|
</pre>
|
|
|
|
<p>and variable <code>reynolds</code> will be loaded from file <em>
|
|
model.oct-data</em> no matter whether it has existed before or not.</p>
|
|
</li>
|
|
|
|
<li>Scilab's simple I/O parallels that of GNU/Octave:
|
|
|
|
<pre>
|
|
save(filename, var1, var2, ..., varN)
|
|
</pre>
|
|
|
|
<p>However, the variables <em>var1</em>, ..., <em>varN</em> are not strings,
|
|
but appear literally. This means that the name of a variable is not stored in
|
|
the file. The association between name and contents is lost!</p>
|
|
|
|
<p>The complementary function</p>
|
|
|
|
<pre>
|
|
load(filename, varname1, varname2, ..., varnameN)
|
|
</pre>
|
|
|
|
<p>restores the contents of <em>filename</em> in the variables named <em>
|
|
varname1</em>, <em>varname2</em>, ... <em>varnameN</em>.</p>
|
|
</li>
|
|
|
|
<li>Tela lets the users save her variables with the
|
|
|
|
<pre>
|
|
save(filename, varname1, varname2, ..., varnameN)
|
|
</pre>
|
|
|
|
<p>function, preserving the association between variable name and variable
|
|
contents. The complementary</p>
|
|
|
|
<pre>
|
|
load(filename)
|
|
</pre>
|
|
|
|
<p>function loads all variables stored in <em>filename</em>. It is not
|
|
possible to select specific variables.</p>
|
|
</li>
|
|
</ul>
|
|
|
|
<h3><a name="matrix oriented i/o">Matrix Oriented I/O</a></h3>
|
|
|
|
<p>As we use matrices so often, specialized functions exist to load and save
|
|
whole matrices. Especially loading a matrix with a single command is
|
|
convenient and efficient to read data from experiments or other programs.</p>
|
|
|
|
<p>Let us assume, we have the ASCII file <em>datafile.ascii</em> which
|
|
contains the lines</p>
|
|
|
|
<pre>
|
|
# run 271
|
|
# 2000-4-27
|
|
#
|
|
# P/bar T/K R/Ohm
|
|
# ====== ====== ======
|
|
19.6 0.118352 0.893906e4
|
|
15.9846 0.1 0.253311e5
|
|
39.66 0.378377 0.678877e4
|
|
13.6 0.752707 0.00622945e4
|
|
12.4877 0.126462 0.61755e5
|
|
</pre>
|
|
|
|
<p>and sits in the current working directory. The file's five leading lines
|
|
are non-numeric. They are skipped by the workbenches, but possibly aid the
|
|
user in identifying her data. I have intentionally taken a data set which is
|
|
not neatly formatted, as are most data files. Matrix-loading functions split
|
|
the input at whitespace not at a specific column, thus they are happy with
|
|
<em>datafile.ascii</em>.</p>
|
|
|
|
<p>We load the data into GNU/Octave with</p>
|
|
|
|
<pre>
|
|
octave:1> data = load("datafile.ascii")
|
|
data =
|
|
</pre>
|
|
|
|
<pre>
|
|
1.9600e+01 1.1835e-01 8.9391e+03
|
|
1.5985e+01 1.0000e-01 2.5331e+04
|
|
3.9660e+01 3.7838e-01 6.7888e+03
|
|
1.3600e+01 7.5271e-01 6.2294e+01
|
|
1.2488e+01 1.2646e-01 6.1755e+04
|
|
</pre>
|
|
|
|
<p>or into Scilab</p>
|
|
|
|
<pre>
|
|
-->data = fscanfMat("datafile.ascii")
|
|
data =
|
|
</pre>
|
|
|
|
<pre>
|
|
! 19.6 0.118352 8939.06 !
|
|
! 15.9846 0.1 25331.1 !
|
|
! 39.66 0.378377 6788.77 !
|
|
! 13.6 0.752707 62.2945 !
|
|
! 12.4877 0.126462 61755. !
|
|
</pre>
|
|
|
|
<p>or into Tela</p>
|
|
|
|
<pre>
|
|
>data1 = import1("datafile.ascii")
|
|
>data1
|
|
#( 19.6, 0.118352, 8939.06;
|
|
15.9846, 0.1, 25331.1;
|
|
39.66, 0.378377, 6788.77;
|
|
13.6, 0.752707, 62.2945;
|
|
12.4877, 0.126462, 61755)
|
|
</pre>
|
|
|
|
<p>In all three examples data will contain a 5-times-3 matrix with all the
|
|
values from <em>datafile.ascii</em>.</p>
|
|
|
|
<p>The complementary commands for saving a single matrix in ASCII format
|
|
are</p>
|
|
|
|
<pre>
|
|
save("data.ascii", "data") # GNU/Octave
|
|
fprintfMat("data.ascii", data, "%12.6g") // Scilab
|
|
export_ASCII("data.ascii", data) // Tela
|
|
</pre>
|
|
|
|
<p>Note that Scilab's <code>fprintfMat()</code> requires a third parameter
|
|
that defines the output format with a C-style template string.</p>
|
|
|
|
<p>Of course none of the above save commands writes the original header, the
|
|
lines starting with hash-symbols, of <em>datafile.ascii</em>. To write these,
|
|
we need the ``low-level'', C-like input/output functions, which featured in
|
|
each of the three workbenches.</p>
|
|
|
|
<h3><a name="clike input/output">C-like Input/Output</a></h3>
|
|
|
|
<p>For a precise control of the input and the output, C-like I/O models are
|
|
offered. All three applications implement function</p>
|
|
|
|
<pre>
|
|
printf(format, ...)
|
|
</pre>
|
|
|
|
<p>Moreover, GNU/Octave and Tela follow the C naming scheme with their C-style
|
|
file I/O:</p>
|
|
|
|
<pre>
|
|
handle = fopen(filename)
|
|
fprintf(handle, format, ...)
|
|
fclose(handle)
|
|
</pre>
|
|
|
|
<p>whereas Scilab prefixes these functions with an ``<code>m</code>'' instead
|
|
of an ``<code>f</code>''</p>
|
|
|
|
<pre>
|
|
handle = mopen(filename)
|
|
mprintf(handle, format, ...)
|
|
mclose(handle)
|
|
</pre>
|
|
|
|
<p>Whether the function is called <code>fprintf()</code> or <code>
|
|
mprintf()</code>, they work the same way.</p>
|
|
|
|
<p><EM>Next Month: Graphics, function plotting and data plotting.</EM></p>
|
|
|
|
|
|
|
|
|
|
|
|
<!-- *** BEGIN bio *** -->
|
|
<SPACER TYPE="vertical" SIZE="30">
|
|
<P>
|
|
<H4><IMG ALIGN=BOTTOM ALT="" SRC="../gx/note.gif">Christoph Spiel</H4>
|
|
<EM>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
|
|
<A HREF="mailto:cspiel@hammersmith-consulting.com">cspiel@hammersmith-consulting.com</A>.</EM>
|
|
|
|
<!-- *** END bio *** -->
|
|
|
|
<!-- *** BEGIN copyright *** -->
|
|
<P> <hr> <!-- P -->
|
|
<H5 ALIGN=center>
|
|
|
|
Copyright © 2001, Christoph Spiel.<BR>
|
|
Copying license <A HREF="../copying.html">http://www.linuxgazette.com/copying.html</A><BR>
|
|
Published in Issue 70 of <i>Linux Gazette</i>, September 2001</H5>
|
|
<!-- *** END copyright *** -->
|
|
|
|
<!--startcut ==========================================================-->
|
|
<HR><P>
|
|
<CENTER>
|
|
<!-- *** BEGIN navbar *** -->
|
|
<IMG ALT="" SRC="../gx/navbar/left.jpg" WIDTH="14" HEIGHT="45" BORDER="0" ALIGN="bottom"><A HREF="mcgucken.html"><IMG ALT="[ Prev ]" SRC="../gx/navbar/prev.jpg" WIDTH="16" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="index.html"><IMG ALT="[ Table of Contents ]" SRC="../gx/navbar/toc.jpg" WIDTH="220" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><A HREF="../index.html"><IMG ALT="[ Front Page ]" SRC="../gx/navbar/frontpage.jpg" WIDTH="137" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="http://www.linuxgazette.com/cgi-bin/talkback/all.py?site=LG&article=http://www.linuxgazette.com/issue70/spiel.html"><IMG ALT="[ Talkback ]" SRC="../gx/navbar/talkback.jpg" WIDTH="121" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><A HREF="../faq/index.html"><IMG ALT="[ FAQ ]" SRC="./../gx/navbar/faq.jpg"WIDTH="62" HEIGHT="45" BORDER="0" ALIGN="bottom"></A><A HREF="tranter.html"><IMG ALT="[ Next ]" SRC="../gx/navbar/next.jpg" WIDTH="15" HEIGHT="45" BORDER="0" ALIGN="bottom" ></A><IMG ALT="" SRC="../gx/navbar/right.jpg" WIDTH="15" HEIGHT="45" ALIGN="bottom">
|
|
<!-- *** END navbar *** -->
|
|
</CENTER>
|
|
</BODY></HTML>
|
|
<!--endcut ============================================================-->
|