Обзор свободного математического ПО

Это конспект моего доклада на семинаре, организованном нашей LUG совместно с университетом. Соответственно, я не мог охватить всё - у меня на доклад было где-то 15 минут.

Вступление

Известные пакеты - это гиганты всё-в-одном

Когда мы говорим о математическом ПО, на ум приходят такие гиганты, как Maple, Mathematica, MatLAB… У них есть одно общее свойство: они пытаются охватить всё. Конечно, Mathematica известна прежде всего как система для символьных вычислений, а Matlab - для численных, но одновременно в Mathematica есть мощные алгоритмы для вычислений с плавающей точкой, а в Matlab - пакет для символьных вычислений. Причём эти второстепенные функции в программах по сравнению с программами, для этого предназначенными, выглядят убого и смешно. А небезызвестный MathCAD пытается включить в себя всё, при этом всё реализовано так себе. Причина проста: нельзя объять необъятное.

Свободные программы - делают одно дело хорошо

В противоположность этому, большинство свободных программ следует философии UNIX, гласящей: программа должна делать одно дело, но делать его хорошо. Свободного математического ПО очень много, при этом бóльшая часть их предназначена для какой-нибудь одной задачи. Например, есть программы, которые только и умеют, что строить сетку для метода конечных разностей. Или программа, которая предназначена для вычисления цифр числа Пи. Или программа, которая умеет только строить графики, но зато очень хорошо.

Однако, есть и программы, в той или иной степени являющиеся аналогами известных пакетов. Я расскажу о трёх.

Символьные вычисления: Maxima

История проекта

Начну я с истории этого проекта.

Сначала я напомню, что компьютеры - это, вообще-то, Электронные Вычислительные Машины, они создавались для вычислений над числами. Однако уже в конце 50-х появилась идея, что можно заставить компьютер работать не только с числами, но и с алгебраическими выражениями. В начале 60-х начали появляться первые системы компьютерной алгебры. И, конечно, такая система нужна была одному мирному американскому ведомству (департаменту энергетики, это практически подразделение Пентагона). Был объявлен тендер, и его выиграл проект под названием Macsyma (пишется через CS). В течение многих лет DOE Macsyma развивалась как коммерческий проект, финансируемый правительством. В 1982-м году Уильям Шелтер создал форк Macsyma, называемый Maxima. В начале 90-х распался СССР, кончилась холодная война, и косвенным следствием этого стало практически полное прекращение финансирования DOE Macsyma. К концу 90-х проект практически загнулся. Исходники Macsyma по кусочкам распродали, и они оказались в Maple и Mathematica. В 1998-м Уильям Шелтер добился от DOE разрешения на публикацию исходных текстов Maxima под лицензией GPL. Maxima стала свободной программой. В 2001-м Шелтер скончался, но к этому моменту над Maxima работало уже довольно много людей, и они подхватили проект.

Интерфейс: командная строка или wxMaxima

Maxima имеет традиционный для UNIX интерфейс командной строки, однако также умеет слушать сетевой порт, работая как сервер. Этот факт используют различные оболочки (фронтенды), предоставляющие графический интерфейс. Наиболее распространены TeXmacs и wxMaxima. TeXmacs - это научный текстовый редактор, в котором можно в документ вставить сессию Maxima. wxMaxima выглядит примерно так:

Последняя версия, 0.8.0, стала больше походить на Mathematica и Maple: раньше командная строка для ввода была отдельно, внизу.

Lisp-подобный язык

Язык Maxima берёт основные идеи из Lisp, так как Maxima написана на Lisp-e. При этом он похож одновременно на языки Mathematica и Maple, так как эти программы позаимствовали многие идеи и часть кода из Macsyma. Чтобы избежать долгого и нудного перечисления возможностей, я приведу пример решения типичных задач с первого курса.

Пример

Пусть дана функция

maxima>> f(x) := x*tanh(x) + x + 1/x + 2;

Проверим, не является ли она чётной или нечётной:

maxima>> f(-x);

Как видим, функция не является ни чётной, ни нечётной. Найдём пределы функции на плюс-минус бесконечности:

maxima>> limit(f(x),x,-inf);

maxima>> limit(f(x),x,inf);

Итак, на плюс бесконечности функция уходит в бесконечность. Нет ли у неё наклонной асимптоты?

maxima>> limit(f(x)/x, x,inf);

Наклонная асимптота есть - y=kx+b, причём k=2. Найдём b:

maxima>> limit(f(x)-2*x, x,inf);

Наконец, построим график:

maxima>> plot2d(f(x), [x,-5,5], [y,-10,10]);

Найдём производную нашей функции:

maxima>> diff(f(x),x);

И заодно - неопределённый интеграл:

maxima>> integrate(f(x), x);

Интеграл до конца "не взялся". Можно показать, что этот интеграл в элементарных функциях и не берётся. Однако Maxima умеет брать некоторые из таких интегралов, используя специальные функции:

maxima>> part: risch(x/(exp(2*x)+1), x);

(здесь я присваиваю результат интегрирования переменной part). Таким образом, интеграл f(x) будет равен

maxima>> ir: -2*part + log(x) + x^2 + 2*x;

Что-то ужасное. Раскроем скобки:

maxima>> expand(ir);

Дифференциальные уравнения

Или вот пример более сложных вычислений. Пусть надо решить дифференциальное уравнение:

maxima>> eq: 'diff(y,x) + x*y = 1-x^2;

Знак апострофа здесь используется, чтобы указать, что не надо сейчас вычислять производную, а сохранить обозначение.

maxima>> solution: ode2(eq,y,x);

Вот и решение. erf здесь - это специальная функция, известная как функция ошибки Лапласа. После раскрытия скобок получим вот что:

maxima>> expand(solution);

По Maxima есть некоторое количество русскоязычных руководств, которые можно найти в интернете. На мой взгляд, самое удачное введение с обзором возможностей содержится в цикле статей Тихона Тарнавского в журнале LinuxFormat. Сейчас эти статьи выложены в открытый доступ, в том числе на русском сайте Maxima. Документация по продвинутым возможностям maxima существует, к сожалению, только на английском языке. Официальная документация составляет 712 страниц.

Численные вычисления: Scilab

Scilab совместим с MatLAB-ом

Наиболее известный пакет для численных расчётов - это MatLAB. Scilab создавался как конкурент matlab-а, более скромный по ценовой политике. Однако коммерчески проект себя не оправдал, и исходные коды были открыты под лицензией, похожей на GNU GPL. Язык scilab сделан по возможности совместимым с матлабом, так что большинство ваших наработок из matlab заработают в scilab. Только вот, как известно, основная мощь matlab-a сосредоточена в его тулбоксах - отдельно поставляемых модулях. Модули для scilab-а тоже есть, однако их сильно меньше.

Octave - это GPL-аналог Matlab

Позже появился проект GNU Octave, нацеленный на создание аналога matlab-a, распространяемого по GNU GPL без всяких заморочек. Язык тоже практически совместим с матлабом, но здесь нет аналога Simulink - средства моделирования и симулирования динамических систем.

Зато Octave имеет чисто консольный интерфейс (конечно, графические фронтенды тоже есть, самый развитый - QtOctave), что позволяет использовать его в скриптах, для автоматизации расчётов, и упрощает встраивание в сложные программные комплексы. Для Octave написаны десятки пакетов расширений.

По Scilab есть статьи на русском языке, кроме того, не так давно в издательстве AltLinux вышла книга `Scilab: Решение инженерных и математических задач'. Книгу можно приобрести в интернет-магазине, кроме того, её электронная версия свободно доступна на сайте AltLinux.

Обработка данных: GNU R

Обзор

Формально, средства обработки данных относятся к программам для численных расчётов, ибо всё что они делают - это вычисления над числами. Однако, как известно, специализированный инструмент всегда лучше универсального. Под словами обработка данных скрывается довольно много различных видов деятельности: статистический анализ, статистическое моделирование, выборка только нужных данных, преобразование данных, построение различных графиков и гистограмм.

Программы для обработки данных можно разделить по типичному размеру выборки, для которого они предназначены. Для небольших выборок подойдёт, например, Statistica. Для средних по размеру выборок хорошо подходит GNU R (она хранит все данные в оперативной памяти, так что на типичном PC получим ограничение в 1-2-4 гигабайта). Для больших и очень больших объёмов данных (от сотен гигабайт до сотен терабайт) предназначены разработанные в CERN свободные системы PAW и ROOT.

GNU R - это интерпретируемый язык программироваммирования, предназначенный для статистического анализа и моделирования. R - это свободная реализация давно существующего языка S. Язык этот весьма эклектичен, он местами похож на C, местами - на Python, местами - на Haskell. Для GNU R существует почти полторы тысячи пакетов расширений (написанных на самом R, на C или Fortran), собранных в репозитории CRAN (Comprehensive R Archive Network).

Типы данных - числа, строки, факторы, векторы, списки и таблицы данных

Основные типы данных в языке - это числа, строки, факторы, векторы, списки и таблицы данных (data frames). Фактор - это данные, которые могут принимать одно из нескольких значений (пол; сорт дерева; логический тип и др). Векторы являются аналогами массивов - это набор из нескольких значений одного типа, размер вектора меняться не может. Тут же надо заметить, что в R нету скаляров; например, число - это, с точки зрения R, вектор из одного элемента. Списки - это обобщение векторов, они могут содержать объекты разных типов, и длина их может меняться. Кроме того, отдельным элементам списка можно присвоить имена, и обращаться к элементам не по номерам, а по именам. Пример:

R>> lst <- list(1,2,3)

(присваивание в R обозначается обычно знаком ←, хотя можно использовать и более привычное =; кроме того, есть форма value → variable). Для обращения к элементам списка по номеру используются двойные квадратные скобки:

R>> lst[[2]]

[1] 2

Назначим имена элементам списка:

R>> names(lst) <- c('first','second','third')

(функция c создаёт векторы). Теперь к элементам списка можно обращаться по именам:

R>> lst$third

[1] 3

Таблица данных (фрейм данных) в R - это список, состоящий из векторов. Создаются таблицы данных чаще всего загрузкой из внешнего файла.

Пример

Скажем, в файле airquality.dat находятся данные замеров качества воздуха:

"Ozone" "Solar.R" "Wind" "Temp" "Month" "Day"
"1" 41 190 7.4 67 5 1
"2" 36 118 8 72 5 2
"3" 12 149 12.6 74 5 3
"4" 18 313 11.5 62 5 4
"5" NA NA 14.3 56 5 5
"6" 28 NA 14.9 66 5 6
"7" 23 299 8.6 65 5 7
"8" 19 99 13.8 59 5 8
"9" 8 19 20.1 61 5 9
"10" NA 194 8.6 69 5 10
.......................

В первой строке - названия полей, дальше идут сами данные. Пропущенные (неизвестные) данные обозначены как NA. Загрузим эти данные в R:

R>> air <- read.table('airquality.dat', sep=' ', header=TRUE)

Здесь мы указываем имя файла, разделитель (пробел), а также указываем, что в первой строке записаны имена полей. К полям таблицы мы можем теперь обращаться как к элементам списка - например, air$Ozone. Посмотрим, что R знает о структуре наших данных:

R>> str(air)

'data.frame':       153 obs. of  6 variables:
$ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
$ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
$ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
$ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
$ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
$ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

Теперь мы можем, например, посмотреть описательную статистику по всем полям таблицы:

R>> summary(air)

    Ozone           Solar.R           Wind             Temp
Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00
1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00
Median : 31.50   Median :205.0   Median : 9.700   Median :79.00
Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88
3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00
Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00
NA's   : 37.00   NA's   :  7.0
Month            Day
Min.   :5.000   Min.   : 1.00
1st Qu.:6.000   1st Qu.: 8.00
Median :7.000   Median :16.00
Mean   :6.993   Mean   :15.80
3rd Qu.:8.000   3rd Qu.:23.00
Max.   :9.000   Max.   :31.00

Для каждого поля показаны минимум, максимум, медиана и две квартили, среднее значение и количество пропущенных данных. Осталось только среднеквадратичное отклонение:

R>> sd(air)

Ozone  Solar.R     Wind     Temp    Month      Day
NA       NA 3.523001 9.465270 1.416522 8.864520

Как видим, R считает среднеквадратичное отклонение для полей Ozone и Solar.R неизвестным - из-за того, что в этих полях есть пропущенные данные. Мы можем явно указать, что на пропущенные данные не надо обращать внимание:

R>> sd(air, na.rm=TRUE)

    Ozone   Solar.R      Wind      Temp     Month       Day
32.987885 90.058422  3.523001  9.465270  1.416522  8.864520

Построим простейшую линейную модель - исследуем зависимость концентрации озона от температуры:

R>> ot <- lm(Ozone ~ Temp, data=air)

R>> ot

Call:
lm(formula = Ozone ~ Temp, data = air)

Coefficients:
(Intercept)         Temp
-146.995        2.429

То есть, если приближать зависимость линейной Ozone = k*Temp + b, то k=2.429, а b=-146.995, при увеличении температуры концентрация озона в среднем растёт.

По GNU R есть довольно много материалов на русском, в частности, методические рекомендации по лабораторным работам для вузов. Также есть хорошее введение в R, содержащееся в цикле статей А.Б. Шипунова и Е.М.Балдина в журнале LinuxFormat, сейчас эти статьи есть в открытом доступе. Продвинутая документация, к сожалению, только на английском, зато её много, включая толстые книги. Официальное руководство к R занимает 2541 страницу.