Числовой инструментарий


Автор: © Christoph Spiel
Перевод: © Роман Шумихин.


Некоторые думают, что GNU/Linux хорошая операционная система, но для неё нет достаточного количества приложений, чтобы преуспеть на рынке. Это может быть и верно для desktop систем, но совершенно неверно по отношению к числовым инструментальным средствам. В этой области у пользователей GNU/Linux огромный и превосходный) выбор, даже слишком огромный чтобы описать все приложения. Поэтому данная серия статей расскажет вам о трёх выдающихся приложениях:

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

Чтобы узнать больше о числовых приложениях смотрите http://sal.kachinatech.com/A/2/

Введение

Что могут делать все эти программы? Разве не достаточно бумаги и карандаша или электронной таблицы?

Основные задачи числовых приложений:

Как бы то ни было, все они предоставляют пользователю полные языки программирования, которые более того могут быть легко расширены, поэтому количество числовых задач, которые они могут решать почти безгранично.

Числовая математика

И зачем же нам нужна эта числовая математика? Числовая математика -- это раздел математики, который разрабатывает, анализирует и применяет методы для выполнения вычислений с конечной точностью. Например, компьютеры как раз используют числовую математику.

Почему компьютеры работают с числами конечной точности? Почему никто не разработал схему, которая позволяет хранить точные числа?

  1. Решение задач с числами с конечной точности быстрее -- намного, намного быстрее. Возьмём для примера сумму всех квадратных корней от 1 до миллиона. На моём компьютере, выполняя точное вычисление в программе MuPAD-1.4.2, доступной по адресу http://math-www.uni-paderborn.de/MuPAD/index.html,

        time(sum(sqrt(i), i = 1..10^6));
    

    занимает около 40 секунд, тогда как вычисление приближённого результата в программе Tela-1.32

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

    занимает 0.31 секунды, то-то и оно, ответ конечной точности получен более чем в 100 раз быстрее! Другими словами, мы можем обработать в сто раз больше данных ограниченной точности за один отрезок времени.

  2. Если использовать хорошие алгоритмы, те которые разработали числовые математики, и быть достаточно осторожным, мы можем получать удивительно точные ответы даже с числами конечной точности.
  3. Признайтесь, большую часть времени вам не нужны точные результаты! Хорошего приближения с "достаточным количеством правильных цифр" будет достаточно.

Организация статьи

В этой серии статей мы рассмотрим общие черты трёх приложений, которые мы собирались обсудить. Мы будем использовать GNU/Octave в большинстве примеров. Там, где есть важные различия, на которые нужно обязательно обратить внимание, мы добавим параграф Особенности в конце раздела.

Технические детали для очень любопытных будут помещены в параграф Детальная информация.

Запуск и выход из приложений

Чтобы дать вам практический опыт, давайте запустим каждое приложение, запросим справку по функции и выйдем из программы.

GNU/Octave
    [email protected]:~/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
[email protected]:~/articles/numerical-workbenches $

Вы можете использовать команду exit или нажать C-d, чтобы выйти из программы GNU/Octave.

GNU/Octave предлагает пользователю автоматическое дополнение имён функций, т.е. если введена часть имени функции, при нажатии пользователем на клавишу Tab, имя функции дополняется насколько это возможно. Повторное нажатие клавиши Tab покажет список всех оставшихся вариантов дополнения команды.

Scilab
После такого запуска Scilab:
    [email protected]:~/articles/numerical-workbenches $ scilab

мы получаем новое X-window, в котором выполняется интерпретатор Scilab. Запрос о помощи открывает окно xless(1x). (Обе эти ссылки - на скрин-шоты, нажмите на них.)

Чтобы выйти из Scilab, введите quit или exit.

Scilab также может быть запущена не в оконном режиме, с помощью ключа -nw в командной строке:

    [email protected]:~/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 довольно краток, однако, справочная система объёмна на сколько это необходимо. Tela также предлагает автоматическое дополнение имён функций, как и в GNU/Octave.
    [email protected]:~/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.
[email protected]:~/articles/numerical-workbenches $

Выйти из Tela можно нажатием C-d.

Лучше, чем карманный калькулятор!

Теперь, когда мы знаем, как запускать и закрывать программы, давайте посмотрим на них в действии.

Простые выражения

Мы хотим увидеть:

  1. Можем ли мы записывать математические выражения таким же способом, как мы привыкли в школе. Да!
    1 + 2 * 3 ^ 4 следует понимать как 1 + (2 * (3 ^ 4)), получаем ответ 163. Не следует воспринимать это выражение как ((1 + 2) * 3) ^ 4, что равно 6561,
  2. Сколько битов нужно чтобы сохранить число 10^6, и
  3. Насколько крута наша дорога? (измеряя в градусах) Наш гараж в 7 метрах от улицы и на полметра выше её.

Поехали.

    [email protected]:~/articles/numerics $ octave

Все три программы управляются из консоли.Т.е. пользователь получает подсказку, когда приложение готово принять ввод. Мы введём наш первый вопрос так, как мы записали его на бумаге. Нажатие клавиши Enter завершает строку, программа оценивает её и возвращает результат в переменной ans (более подробная информация о переменных будет дальше).

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

Ага, очевидно GNU/Octave знает элементарную школьную математику!

Наш второй вопрос требует логарифмической функции log, которая возвращает натуральный логарифм своего аргумента; это логарифм по основанию e.

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

Мы заключаем, что для хранения 1,000,000 нужно 20 бит.

Наконец, насколько крут наш выезд из гаража? Что нам здесь нужно - это угловая функция, именно арк-тангенс, записывается это так atan(argument).

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

Хмм, не кажется ли вам, что получилось слишком плоско? Из давно забытых уроков математики, мы помним, что арк-тангенс 1 равен 45 градусам. Давайте проверим это!

    octave:4> atan(1)
ans = 0.78540

Опс, ответ получился в 57 раз меньше! Что же нам теперь выбрасывать программу? Подождите, 57 при умножении на число пи равняется почти 180. Это означает, что GNU/Octave возвратила результат в радианах, а не в градусах. Все угловые функции работают с радианами, угол в 360 градусов равен 2 пи радиан.

Попробуем снова, но на этот раз преобразуем результат из радиан в градусы:

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

Около 4-х градусов, выглядит неплохо. Наш гараж точно не будет затоплен при сильном ливне.

Детальная информация

Особенности

Переменные

В предыдущей секции мы не получили ничего нового по сравнению с карманным калькулятором, но та черта, ,благодаря которой наши программы "сделают" карманные калькуляторы и записные таблицы -- это имена, которые мы можем давать параметрам или результатам вычислений, они называются переменными.

Положим, наша лучшая половина хочет построить сад во дворе, но мы хотим смотреть баскетбол. Поэтому нам быстро нужна какая-нибудь отмазка, например, что у нас недостаточно компоста на желаемый размер сада. Ха -- замечательная идея!

[ПХЯСМНЙ: ОКЮМ ДКЪ ЙКСЛАШ]

Из нашего маленького плана мы берём следующие длины в футах (1 фут = 12 дюймов = 30,48 см):

    houseside_length = 10 (длина со стороны дома)
creekside_length = 6  (длина со стороны ручья)
width = 2             (ширина)

Наша лучшая половина также сказала, что слой новой почвы должен быть, по крайней мере, 5 дюймов (12,7 см), поэтому

    height = 5 / 12 (высота)

Призываем на помощь GNU/Octave!

    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

Ящик с компостом 6' x 4' (в дюймах) и сейчас содержит 8 дюймов полезного компоста.

    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    (объём компоста)

О нет, мы только что вырыли себе могилу. Мы получили достаточно компоста! Может быть записать матч на видео?

Детальная информация

Структурированные данные

До настоящего времени мы не использовали компьютер там, где он действительно хорош: периодически повторяющаяся работа.

Векторы

Скажем, мы получили длинный счёт из бакалейной лавки. [Здесь должна быть ваша реклама!] Как мы можем посчитать НДС в долларах на каждый предмет, давая общее количество, и ставку НДС в процентах? Формула

            процент_ндс / 100
ндс = --------------------- * общее_количество
1 + процент_ндс / 100

тривиальна, но мы не хотим набирать её снова и снова.

Список общего количества в счёте формирует то, что числовые программы называют вектором. Векторы строятся из значений, заключая их в квадратные скобки и разделяя запятыми между собой, вот так:

    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

Вектор строится слева направо, используя предоставленные числа в таком же порядке, в каком они были введены.

Разве не будет чудесно, если мы просто напишем: gross *(vat_percent/100) / (1 + vat_percent/100) и получим НДС для каждого из них? Это действительно просто.

    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    (gross - общее количество)
vat 
      0.097477  0.325794  0.051682  0.359159  0.062804  0.062804  0.062804  0.062804

Ура -- работает! Впервые мы получили удобство и выразимость: один знак умножения выполняет подряд 8 умножений.

Что произошло? vat_percent -- это единственное значение, которое называется в теории чисел скаляром, чтобы отличить его от векторов. Хорошо, если vat_percent скаляр, тогда vat_percent/100, 1 + vat_percent/100 и a -- тоже скаляры. Окончательно, скаляр a должен быть умножен на вектор gross. Что мы хотели и получили -- это a, умноженная на каждый компонент вектора gross. Это касается любого оператора, не только умножения! В основном,

     вектор оп  скаляр

и

     скаляр оп  вектор

применяют скаляр к каждому компоненту вектора в соответствии с операцией оп. В нашем примере, как если бы мы написали следующее

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

где мы представили новый кусок синтаксиса: индексы векторов. Каждый компонент (скаляр) вектора может быть получен по его индексу, который является номером его места в векторе. Например, чтобы получить второй элемент gross, мы пишем

    octave:5> gross(2)
ans = 4.9800

Присваивать компоненты вектора можно таким же образом. Просто поместите проиндексированный вектор с левой стороны от знака присваивания, например, gross(2)= 5.12.

О чём ещё можно думать как о векторе чисел, кроме нашего чека? Любые последовательности значений! Большую часть времени значения будут относительными, как температура, измеряемая каждый день в 8 часов утра, как диаметры стопки железных прутьев, скорости всех движущихся в западном направлении через Buffalo Street, на углу West Clinton Street в среду, 18-го Апреля  2001-го года автомобилей. Поскольку мы живём в Цифровую Эпоху, многие последовательности данных, можно рассматривать как векторы: каждый кусок музыки на CD - это вектор звуковых амплитуд и индексные метки дискретных моментов во времени.

Детальная информация

Особенности

В следующем месяце

Christoph Spiel

Chris управляет консалтинговой компанией по открытому программному обеспечению в верхней Баварии, Германия. По своей специальности - физике -- у него есть докторская степень от Munich University of Technology -- его основные интересы лежат в области теории чисел, гетерогенных программных сред и программного инжиниринга. С ним можно связаться по адресу [email protected].
Copyright © 2001, Cristoph Spiel.
Published in Issue 69 of Linux Gazette, August 2001