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

Автор: © 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/

Введение

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

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

• Предварительная обработка (да -- между тем нам нужно чтобы компьютеры разговаривали с компьютерами) и заключительная обработка данных в случае необходимости "склеивания" двух числовых приложений вместе.
• Числовые оптимизации (обе, и линейные и нелинейные),
• Самостоятельное моделирование,
• Визуализация данных,
• Интеллектуальные "карманные" калькуляторы,
• Быстрое прототипирование специальных числовых приложений, которые будут затем написаны, например, на C++ или Fortran-90.

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

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

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

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

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

    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 $

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

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

Scilab

После такого запуска Scilab:

    cspiel@hydra:~/articles/numerical-workbenches $ scilab

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

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

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

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

    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 можно нажатием C-d.

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

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

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

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

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

Поехали.

    cspiel@orion:~/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-х градусов, выглядит неплохо. Наш гараж точно не будет затоплен при сильном ливне.

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

• Числа могут быть вещественными, либо комплексными. Элементарные операции +'', -'', *'', and /'', вместе с потенцированием ^'', работают, как и ожидается, на вещественных и комплексных числах.
• Часто используемые основные функции - это нахождение абсолютного значения (модуль числа) abs(arg), функция, определяющая знак числа, sign(arg), и квадратный корень sqrt(arg).
• Есть две логарифмические функции, одна с основанием e: log(arg), и с основанием 10: log10(arg). Экспоненциальная функция exp(arg) - это обратная функция к log(arg).
• Все приложения предлагают богатый набор угловых и гиперболических функций: sin(arg), cos(arg), tan(arg), sec(arg), csc(arg); asin(arg), acos(arg), atan(arg), acsc(arg); sinh(arg), cosh(arg), tanh(arg), sech(arg), csch(arg) asinh(arg), acosh(arg), atanh(arg), acsch(arg).

Особенности

• В Tela нет удобной переменной ans.
• В GNU/Octave и Tela мнимые буквенные числа записываются с добавлением "i", "j", "I", или "J" к числу. Например, 1i, -8.99I, 324J. Scilab определяет специальную константу для мнимой единицы sqrt(-1), которая записывается как <%i>. Поэтому, в Scilab мнимые числа выглядят как произведения: -8.99*%i, %i*324.

Переменные

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

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

Из нашего маленького плана мы берём следующие длины в футах (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 - это вектор звуковых амплитуд и индексные метки дискретных моментов во времени.

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

• Говоря математическим языком, все векторы здесь - кортежи.
• Многие встроенные функции, например sin, могут быть использованы на векторах.

      v = [0.12, 0.89, 0.78, 0.10]
  sin(v)
  

  или

      sin([0.12, 0.89, 0.78, 0.10])
  

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

Особенности

• Tela использует другой синтаксис для построения векторов, который напоминает Scheme или Smalltalk: Вектор окружается круглыми скобками, предваряемые знаком решётки, например: gross = #(1.49, 4.98, ..., 0.96).
• Tela использует другой синтаксис для индексирования компонентов вектора, который напоминает Pascal. Индекс заключается в квадратные скобки; например: gross[2]. Предупреждение для программистов на Си: хотя квадратные скобки выглядят как в Си, нижний индекс всегда 1.

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

• Матрицы
• Функции, определённые пользователем
• Операции управления потоками
• Ввод и вывод