CC BY
29УДК 004.67
ТЕХНОЛОГИЯ МНОГОФАКТОРНОГО РЕГРЕССИОННОГО АНАЛИЗА ДАННЫХ ФИЗИЧЕСКОГО И МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ НА ОСНОВЕ
ПРОГРАММЫ GNUPLOT Хан Полина Вениаминовна
к.ф.-м.н, научный сотрудник, Институт систем энергетики им. Л.А. Мелентьева СО РАН, 664130 г. Иркутск, ул. Лермонтова 130, e-mail: polinakhan@gmail. com
Таиров Эмир Асгадович д.т.н., главный научный сотрудник, Институт систем энергетики им. Л.А. Мелентьева СО РАН, 664130 г. Иркутск, ул. Лермонтова 130, e-mail: tairov@isem.irk.ru
Аннотация. В статье представлено описание процесса обработки данных эксперимента по критическому истечению пароводяной смеси через неподвижную засыпку твёрдых частиц. Приводится обзор программных средств регрессионного анализа. На основе существующих математических моделей и экспериментальных данных составлено описание зависимостей коэффициента скольжения фаз и показателя политропы от давления и массового паросодержания на входе в засыпку, и предложены функции, обладающие необходимыми свойствами. Параметры этих функций определены при помощи итерационного алгоритма, в основу которого положен многофакторный регрессионный анализ, доступный в программе gnuplot. Рассмотрены особенности интеграции скрипта gnuplot в скрипт командного интерпретатора bash. Полученные замыкающие соотношения позволили составить предсказательную модель для вычисления критической массовой скорости пароводяной смеси в засыпке из частиц диаметром 2 и 4 мм при начальном давлении от 0.6 до 1.2 МПа и начальном расходном массовом паросодержании от 0.016 до 0.178. Разработанная технология может быть применена для реализации гибкой и эффективной системы обработки любых экспериментальных данных, включающей предварительную подготовку, регрессионный анализ и построение графиков. Ключевые слова: многофакторная нелинейная регрессия, программные средства, критическое истечение, пароводяная смесь, зернистый слой.
Введение. К настоящему времени сложилась вполне эффективная методология построения математических моделей сложных физических и технических систем при существенной степени неизученности протекающих физических процессов.
На первом этапе исходная система идеализируется для исключения эффектов, которые кажутся несущественными или малоизученными. Следующим этапом является построение математической модели на основе применения фундаментальных законов сохранения массы, импульса, энергии в рамках принятых при идеализации исходной системы физических допущений. Как правило, при этом применяется процедура усреднения во времени и пространстве. Далее утраченная при упрощении информация восполняется
вводом системы замыкающих соотношений. На последующем этапе анализа достоверности модели производится сопоставление предсказаний модели с систематическими экспериментальными данными. При этом модель уточняется путём надлежащего подбора содержащихся в замыкающих соотношениях коэффициентов.
Данная работа посвящена проблеме автоматизации процесса регрессионного анализа экспериментальных данных при построении модели критического истечения пароводяной смеси через засыпку из шаровых частиц. Приводятся описание экспериментальной установки, теоретической модели и постановка задачи регрессии. Представлен обзор компьютерных систем, применимых для решения данной задачи. Обоснован выбор программы gnuplot и разобраны ключевые моменты выполнения регрессии и интеграции gnuplot со средой командного интерпретатора при реализации итерационного алгоритма.
1. Эксперимент, теоретическая модель и постановка задачи регрессии. Экспериментальные данные были получены на установке «Высокотемпературный контур», схема которой представлена на рис. 1. Данная установка позволяет моделировать и изучать процессы истечения, теплообмена и распространения звука в многофазных средах, характерные для энергетических установок.
Рис. 1. Схема экспериментальной установки: 1- бак-аккумулятор, 2 - фильтр, 3 - насос, 4 -подогревательный участок, 5 - рабочий участок, 6 - конденсатор, 7 - автоматический измеритель расхода, 8 - дроссельная шайба, 9 - кольцевые отборы давления, 10, 11 -регулирующие вентили, 12 - электромагнитный клапан. Р}, Р2, Р3, Р4 - измерения давления,
Т, Т2, Т3, Т4 - измерения температуры.
Пароводяная смесь, полученная путём дросселирования, проходила снизу-вверх через вертикальный рабочий участок, представляющий собой полый цилиндр со сменными вставками засыпок из жестко закрепленных стальных шариков одинаковых размеров. Входные значения давления смеси и расходного массового паросодержания поддерживались средствами ручного и автоматического регулирования. Массовый расход измерялся объемным способом после конденсации пароводяной смеси. Снижение выходного давления до критических значений обеспечивалось при помощи вентиля выхода в атмосферу.
Путём измерения расхода при постепенном снижении выходного давления были получены значения критического расхода для следующих значений входных параметров: диаметр шаровых частиц засыпки, ё - 2 и 4 мм, высота засыпки, Н - 50, 100, 250 и 355 мм, входное давление, Р1 - 0.6, 0.9 и 1.2 МПа, входное расходное массовое паросодержание, х1 -0.016, 0.022, 0.033, 0.055, 0.096 и 0.178 [7].
Для описания опытных данных по критическому истечению парожидкостной смеси через неподвижный слой шаровых частиц использована модель газовой динамики зернистого слоя [4]. Согласно этой модели, при й << Я, критическая массовая скорость определяется следующей формулой:
где р - плотность, ^ - среднеобъемная скорость среды в пористом пространстве, т — пористость среды, = 0.508 — 0. 56 ( 1 — т) — относительное минимальное проходное сечение, а п - показатель политропы. Переменные с индексом 1 относятся ко входному сечению, а без индекса - к произвольному сечению, либо засыпке в целом.
В ранее опубликованной работе [7] было установлено соответствие результатов эксперимента формуле (1) в части пропорциональности критической массовой скорости фактору ( ) . В
данной работе формула (1) была адаптирована для пароводяной смеси путём определения входной плотности и показателя политропы. Плотность равновесной пароводяной смеси определяется следующими соотношениями:
[. Р" 1-х~1~1 ...
1 + ] , (2)
р = р ' 'ф + р '( 1 —<р ), (3)
где ф - истинное объемное паросодержание, а 5 = и' '/и/' - коэффициент скольжения фаз. Переменные с одним штрихом относятся к насыщенной жидкости, переменные с двумя штрихами — к насыщенному пару, а без штриха — к смеси. В данной работе была решена задача определения замыкающих соотношений для коэффициента скольжения ^ и показателя политропы пароводяной смеси п.
Показатель политропы при заданном коэффициенте скольжения определялся следующим образом. Процесс расширения смеси при прохождении через засыпку представляет собой многократное дросселирование, а потому протекает при постоянной энтальпии. Это условие позволяет определить расходное массовое паросодержание в произвольном сечении:
к' 1 ( 1 - )+к' ' 1 - к' , X = —-Ц-—-. (4)
На основе формул (2-4) и имеющегося коэффициента скольжения в каждом сечении, определяемого текущим давлением Р от входного до атмосферного, была вычислена плотность смеси □. Полученный таким образом набор пар (РьР{) для заданных входных условий аппроксимировался уравнением политропы
(Р/Р1) = (р/р1 )". (5)
Задача регрессии формулируется следующим образом: найти соотношения для коэффициента скольжения и показателя политропы на основе данных по величинам критических массовых скоростей при различных значениях параметров засыпки и состояниях смеси на входе с учётом уравнений (1-5).
2. Обзор программных средств регрессионного анализа. Регрессионный метод идентификации моделей широко используется в самых различных областях знания, и современному исследователю доступны как алгоритмы регрессии, так и многочисленные программные средства, их реализующие. Регрессия делится на линейную и нелинейную, а также на однофакторную и многофакторную. Линейная регрессия позволяет вычислить неизвестные коэффициенты в явном виде. Нелинейная требует применения итерационных оптимизационных подходов. Однофакторная позволяет обойтись скалярными вычислениями, тогда как многофакторная требует выполнения векторных операций. В рассматриваемой задаче критическая массовая скорость зависит от четырёх параметров, а искомые зависимости для коэффициента скольжения и показателя политропы — от двух. При этом все зависимости являются нелинейными. Поэтому необходимо использовать программные средства, решающие задачу многофакторной нелинейной регрессии.
Программные средства, включающие методы регрессионного анализа, можно условно подразделить на электронные таблицы, системы компьютерной математики, среды статистической обработки данных и библиотеки процедур для языков программирования.
Электронные таблицы являются самым популярным средством регрессионного анализа данных, поскольку они установлены практически на каждом компьютере: Excel [18] входит в состав платного пакета Microsoft Office, a Calc [12, 16], является частью бесплатных пакетов ПО LibreOffice и OpenOffice.org. Excel и Calc предлагают два способа получения коэффициентов регрессии: добавление линии тренда с ее уравнением на график и функции LINEST и LOGEST для прямого вычисления коэффициентов регрессии. В обоих способах отсутствует возможность аппроксимации данных произвольно заданной нелинейной функцией.
Среди систем компьютерной математики (СМК), Mathcad [14], Mathematica [15], Octave [2], Maxima [8] поддерживают все виды регрессии, вплоть до многофакторной нелинейной. В Matlab [13] многофакторная регрессия ограничена обобщенной линейной, а нелинейная регрессия общего вида применяется только для функций одной переменной. Все СМК позволяют создавать программы, выполняемые в их среде, и любые операции, включая регрессию, осуществлять в пакетном режиме. Основным недостатком СМК для решения поставленной задачи является их избыточность. Это очень мощные и многофункциональные программы, и устанавливать их только ради регрессии нерационально. Особенно это относится к платным Mathcad, Mathematica и Matlab, лицензионные версии которых стоят много больше лицензионных версий Microsoft Office.
К средам статистической обработки данных можно отнести R[21], SPSS[6], Stata[20], SAS[19] и Statistica[3]. R - это бесплатная среда анализа данных, основанная на языке программирования S. R имеет все основные управляющие конструкции, скалярные и векторные типы данных, средства работы с файлами, графику и все популярные статистические методы, включая регрессию всех видов. Коммерческие системы статистической обработки SPSS, Stata, SAS и Statistica, а также построенные на основе R коммерческие программные продукты Revolution R, RODM, кроме этих основных возможностей, включают возможность эффективной работы с очень большими массивами данных с использованием параллельных вычислений.
Особняком стоит gnuplot [11]. Это консольное приложение со скриптовым языком имеет широчайшие возможности по графическому представлению данных, включая
контурные диаграммы на произвольных сетках, двумерные и трехмерные диаграммы, диаграммы с погрешностями и многое другое. Кроме построения графиков, gnuplot умеет находить коэффициенты множественной нелинейной регрессии. Язык gnuplot позволяет задать настройки графиков, определения пользовательских функций и значений переменных, и все это сохранять в файлы и загружать из файлов. Чтение файлов данных происходит непосредственно в ходе построения графика, либо выполнения регрессии без использования векторных переменных. Отсутствуют какие-либо управляющие конструкции, кроме цикла со счётчиком внутри функций построения диаграмм. Настройки позволяют строить диаграммы в отдельном окне или в файле любого графического типа, либо в текстовом файле табличного типа.
Также для языков программирования высокого уровня, таких, как Fortran, C++, C#, Java, Python и других, существуют свободно распространяемые и коммерческие библиотеки математических функций, включающих все виды регрессии [9, 17]. И, наконец, достаточно эффективные для большинства приложений алгоритмы известны уже не одно десятилетие, и могут быть реализованы самостоятельно.
В результате обзора программного инструментария для анализа данных эксперимента методом нелинейной регрессии выработаны следующие рекомендации.
1. Пользователям систем компьютерной математики или систем статистического анализа данных стоит воспользоваться функциями регрессии этих систем. Это будет наиболее эффективный метод, как по затратам рабочего времени на программирование, так и по скорости обработки данных.
2. Тем, у кого установлена система программирования, и написание небольших программ не вызывает затруднений, стоит воспользоваться функциями регрессии из хорошо зарекомендовавших себя математических библиотек.
3. В иных случаях стоит установить бесплатные R или gnuplot. При этом изучение R займет больше времени, но позволит полностью автоматизировать процесс обработки данных в пределах одной среды программирования. Изучить gnuplot можно за один день. Но для реализации итерационного процесса придётся задействовать командный интерпретатор, как описано ниже для интерпретатора bash.
3. Реализация и результаты. Как указано выше, процесс расширения смеси является изоэнтальпийным и при этом, для применения теории газовой динамики зернистого слоя, он должен аппроксимироваться политропным. Поэтому функции s(P, х) и n(f^1,x1) должны быть согласованы, и алгоритм определения их коэффициентов становится итерационным.
Общая схема алгоритма приводится на рисунке 1.
Полученные значения n для всех входных данных аппроксимировались формулой вида
n(Xi) = an + (1.0 - aj(l - exp(-xi/bn)) (6)
Выбор формулы обусловлен расположением точек и тем фактом, что, при стремлении паросодержания к единице любая характеристика пароводяной смеси стремится к ее значению для насыщенного пара. Изоэнтальпийный процесс для пара - это изотермический процесс, для которого п = 1. На рис. 3 показаны точки и кривая для показателя политропы после последней итерации при значении коэффициентов регрессии
an = 0.481, Ьп = 0.145. (7)
Рис. 2. Общая схема алгоритма
Отметим, что точки для различных входных значений давления на рис. 3 ложатся на одну кривую. Коэффициент скольжения фаз был определен в виде:
Р, х) = 1 + а5(Р)ехр ( - + С5( р}х)) (8)
где а5, Ь5, с5 - квадратичные функции от давления. Данный выбор обусловлен тем, что при х стремящимся к нулю мелкие пузырьки движутся в потоке жидкости, а при х, стремящимся к единице, мелкие капли движутся в потоке пара [10], и в обоих этих крайних случаях коэффициент скольжения равен единице. В промежуточных состояниях пар движется быстрее жидкости, и 5 > 1. Результат после последней итерации приводится на рисунке 4 при значении коэффициентов регрессии
а5 = 7.3 3 - 8.2 7 Р + 4 . 3 5 Р2,
Ь5 = 0.0 592 -0. 124 Р + 0.0 746 Р2, (9)
с5 = 2 . 7 1 - 1 .9 0 Р + 1 . 3 9 Р2.
Начальное приближение функции Р , х) было выбрано на основе формулы для узких каналов [5], применяемой также для течения в засыпках [1].
3 = (РСг-/Р)° ■ 3 8 (10)
1 г
0.9 -
0.8 -0.7
0.6 -
0.5 }
0.4 -
0 0.2 0.4 0.6 0.8 1 Л1
Рис. 3. Показатель политропы. Точки -опорные значения, линия - уравнение (6).
si
О 0.2 0.4 0.6 0.8 1 Х1
Рис. 4. Коэффициент скольжения. Точки -опорные значения, линии - уравнение (8).
3.5
2.5
1.5
1 Р.-0.6 МРа
Р -0.9 МРа -
f: у Р -1.2 МРа - ............ ■
f * k-. Ч Р =0.6 МРа +
Р =0.9 МРа X
, Р = 1.2 МРа ж
"L......
Необходимость итерационного подбора коэффициентов двух аппроксимирующих функций (6) и (8), обусловила выбор технологии, позволяющей писать скрипты и осуществлять пакетную обработку. Выбор пал на программу gnuplot, как минимальный по сложности бесплатный пакет, включающий многофакторную нелинейную регрессию общего вида. Поскольку gnuplot не содержит управляющих конструкций циклов и ветвления, для реализации описанного выше алгоритма было написано семейство скриптов на языках командного интерпретатора bash, потокового текстового процессора awk и gnuplot. Далее описаны основные возникшие при программировании задачи и методы их решения.
1. Нахождение коэффициентов регрессии по табличным данным осуществляется вызовом функции fit программы gnuplot. Параметры данной функции позволяют указать любую скалярную аппроксимирующую функцию, с одной или несколькими переменными и любым количеством коэффициентов, не превышающим количества точек в данных, имя файла данных, а также любое преобразование столбцов данных.
2. Возвращение коэффициентов регрессии из gnuplot в среду bash. Gnuplot использует стандартный консольный ввод, но не использует консольный вывод. Поэтому использован следующий приём: gnuplot сохраняет коэффициенты регрессии в текстовом файле, а awk извлекает оттуда только числовое значение коэффициентов.
3. Вызов параметризованного скрипта gnuplot из скрипта bash. Поскольку gnuplot не принимает значения переменных командной строки, была использована технология here-document, позволяющая генерировать скрипт для gnuplot прямо в точке вызова с подстановкой значений переменных bash в текст скрипта.
Заключение. С использованием описанного подхода построена предсказательная модель для критического истечения пароводяной смеси через засыпку из шаровых частиц. Во всем рассмотренном диапазоне параметров относительная разница между измеренными и вычисленными значениями составляет не более 7% (рис. 5).
(Pw)cr kg/(m" s) 1000
750
500
250
; p =0.3 MPa +
p =0.6 MPa X
p =0.9 MPa ж
- \ n p = 1.2 MPa □
X П \
■ "------- .....в......... -—ж........
0.00 0.05 0.10 0.15 0.20 шл<> о.ОО 0.05 0.10 0.15 0.20 х'
а) б)
Рис. 5. Критическая массовая скорость: а) в зависимости от фактора (й/Я)05 при постоянном входном давлении; б) в зависимости от входного паросодержания для одной засыпки. Точки - экспериментальные данные, линии - расчет.
Разработанная технология обработки результатов эксперимента на основе gnuplot
позволяет осуществлять регрессию самого общего вида, в том числе итерационно или в
сочетании с другими типами преобразований.
Работа выполнена в ЦКП «Высокотемпературный контур» в Институте систем
энергетики им. Л.А. Мелентьева СО РАН при поддержке РФФИ (проект №17-08-00709а)/
СПИСОК ЛИТЕРАТУРЫ
1. Авдеев А.А., Балунов Б.Ф., Рыбин Р.А., Созиев Р.И., Филиппов Г.А. Гидродинамическое сопротивление при течении двухфазной смеси в шаровой засыпке // Теплофизика высоких температур. 2003. Т. 41, №3. С. 432 - 438. DOI: 10.1023/A:1024246728189.
2. Алексеев Е.Р., Чеснокова О.В. Введение в Octave для инженеров и математиков. ALT Linux, 2012.
3. Боровиков В. STATISTICA: искусство анализа данных на компьютере. С-Пб.: Питер, 2003.
4. Гольдштик М.А. Процессы переноса в зернистом слое. Новосибирск: Изд. Института теплофизики СО РАН, 2005. 358 с.
5. Кириллов П.Л., Бобков В.П., Жуков А.В., Юрьев Ю.С. Теплогидравлические процессы в ЯЭУ: Справочник по теплогидравлическим расчетам в ядерной энергетике / под ред. д.т.н., проф. П Л. Кириллова. Т. 1. М.: ИздАТ, 2010.
6. Руководство пользователя по базовой системе Statistics 20. USA: IBM Corporation, 2011.
7. Таиров Э.А, Покусаев Б.Г., Быкова С.М. Критическое истечение парожидкостного потока через слой шаровых частиц // Теплофизика высоких температур. 2016. Т. 54, №2. С. 277-286.
8. Чичкарев Е.А. Компьютерная математика с Maxima: Руководство для школьников и студентов. М.: ALT Linux, 2012. 384 c.
9. Anderson E., Bai Z., Bischof C., Blackford S., Demmel J., Dongarra J., Croz J. Du, Greenbaum A., Hammarling S., McKenney A., Sorensen D. LAPACK Users' Guide. Third Edition. Society for Industrial and Applied Mathematics, 1999.
10. Delhaye J.M., Giot M., Riethmuller M.L. Thermohydraulics of two-phase systems for industrial design and nuclear engineering. Hemisphere Publishing Corporation, 1981.
11. Janert P.K. Gnuplot in Action. Understanding Data with Graphs. USA: Manning Publications Co., 2010
12. LibreOffice 3. Calc Guide. Using Spreadsheets in LibreOffice 2. GNU Public License, 2011.
13. Martinez W.L., Martinez A.R. Computational Statistics Handbook in MATLAB, 2nd ed. Chapman and Hall. CRC, 2007.
14. Mathsoft Mathcad 11 User Guide. Cambridge, USA: Mathsoft Engineering & Education Inc., 2002.
15. Mure§an M. Introduction to Mathematica with Applications. Springer, 2017.
16. OpenOffice.org3. Calc Guide. Using Spreadsheets in OpenOffice.org. GNU Public License, 2010.
17. Pemstein D., Quinn K. M., Martin A. D., The Scythe Statistical Library: An Open Source C++ Library for Statistical Computation // Journal of Statistical Software. 2011. Vol. 42. Issue 12.
18. Remenyi D., Onofrei G., English J. An Introduction to Statistics using Microsoft Excel. Curtis Farm, UK: Academic Publishing Ltd., 2009.
19. Ron Cody. An Introduction to SAS. University Edition. USA: SAS Institute, 2015.
20. STATA. User's Guide Release 13. StataCorp, 2013.
21. Venables W. N., Smith D. M. and the R Core Team. An Introduction to R. Notes on R: A Programming Environment for Data Analysis and Graphics. Version 3.4.0, 2017.
UDK 004.67
TECHNOLOGY OF MULTIPLE REGRESSION ANALYSIS OF PHYSICAL AND MATHEMATICAL MODELING DATA BASED ON GNUPLOT PROGRAM
Polina V. Khan
Ph.D., Researcher, Melentiev Energy Systems Institute Siberian Branch of the Russian Academy of Sciences 130, Lermontov Str., 664033, Irkutsk, Russia, e-mail: polinakhan@gmail.com
Emir A. Tairov D.Sc., Chief Researcher, Melentiev Energy Systems Institute Siberian Branch of the Russian Academy of Sciences 130, Lermontov Str., 664033, Irkutsk, Russia, e-mail: tairov@isem.irk.ru
Abstract. In the present work, the data processing for the experiment on the critical flow of a vapor-liquid mixture through the densely-packed beds of spherical particles is described. The overview of existing software for regression analysis is given. The effect of the input flow quality on the phase slip coefficient and polytropic coefficient is described and appropriate approximating functions are proposed. The functions parameters are identified by using an iterative program, based on the multiple regression
analysis, available in gnuplot software. The technique of data exchange between command interpreter bash and gnuplot is considered. The closing relationships obtained together with the granular layer theory constitute the critical mass flow rate model for the input pressure from 0.6 to 1.2 MPa and input flow quality from 0.016 to 0.178. The developed approach of bash-gnuplot integration can be used for any experiment data processing and plotting.
Keywords: multiple nonlinear regression, computer programs, critical flow, vapor-liquid mixture, granular layer.
References
1. Avdeev A.A., Balunov B.F., Rybin R.A., Soziev R.I., Filippov G.A. Gidrodinamicheskoe soprotivlenie pri techenii dvuhfaznoj smesi v sharovoj zasypke [Hydrodynamic drag of a flow of two phase mixture in a pebble bed] // High Temperature. 2003, vol. 41, iss. 3, pp. 377-383. DOI: 10.1023/A:1024246728189. (in Russian)
2. Alekseev E.R., Chesnokova O.V. Vvedenie v Octave dlja inzhenerov i matematikov [Introduction to Octave for Engineers and Mathematicians]. Moscow, ALT Linux Publ., 2012. (in Russian)
3. Borovikov V. STATISTICA: iskusstvo analiza dannyh na komp'jutere [Art of data analysis on computer]. Saint-Petersburg, Piter Publ., 2003. (in Russian)
4. Goldstick M.A. Transfer processes in granular layer. Novosibirsk, Kutateladze Institute of Thermophysics Russian Academy of Sciences Publ., 2005. (in Russian)
5. Kirillov P.L., Bobkov V.P., Zhukov A.V., Jur'ev Ju.S. Teplogidravlicheskie processy v JaJeU: Spravochnik po teplogidravlicheskim raschetam v jadernoj jenergetike [Heat transfer and hydraulic processes in Nuclear Power Plants. Reference book on thermohydraulic calculations in nuclear power industry]. / under ed. Dr. Prof. Kirillov P.L., vol. 1, Moscow, IzdAT Publ., 2010. (in Russian)
6. User Guide. Basic Statistics 20. USA, IBM Corporation, 2011.
7. Tairov E.A., Pokusaev B.G., Bykova S.M., Vapor-liquid critical flow through a layer of spherical particles // High Temperature, 2016, vol. 54, iss. 2, pp. 277-286.
8. Chichkarev E.A. Komp'juternaja matematika s Maxima: Rukovodstvo dlja shkol'nikov i studentov [Computer mathematics with Maxima. Guide for students]. Moscow, ALT Linux Publ., 2012. 384 p. (in Russian).
9. Anderson E., Bai Z., Bischof C., Blackford S., Demmel J., Dongarra J., Croz J. Du, Greenbaum A., Hammarling S., McKenney A., Sorensen D. LAPACK Users' Guide. Third Edition. Society for Industrial and Applied Mathematics, 1999.
10. Delhaye J.M., Giot M., Riethmuller M.L. Thermohydraulics of two-phase systems for industrial design and nuclear engineering. Hemisphere Publishing Corporation, 1981.
11. Janert P.K. Gnuplot in Action. Understanding Data with Graphs. USA, Manning Publications Co., 2010
12. LibreOffice 3. Calc Guide. Using Spreadsheets in LibreOffice 2. GNU Public License, 2011.
13. Martinez W.L., Martinez A.R. Computational Statistics Handbook in MATLAB, 2nd ed. Chapman and Hall. CRC, 2007.
14. Mathsoft Mathcad 11 User Guide. Cambridge, USA, Mathsoft Engineering & Education Inc., 2002.
15. Mure§an M. Introduction to Mathematica with Applications. Springer, 2017.
16. OpenOffice.org3. Calc Guide. Using Spreadsheets in OpenOffice.org. GNU Public License, 2010.
17. Pemstein D., Quinn K.M., Martin A.D., The Scythe Statistical Library: An Open Source C++ Library for Statistical Computation // Journal of Statistical Software, 2011, vol. 42. iss. 12.
18. Remenyi D., Onofrei G., English J. An Introduction to Statistics using Microsoft Excel. Curtis Farm, UK, Academic Publishing Ltd., 2009.
19. Ron Cody. An Introduction to SAS. University Edition. USA, SAS Institute, 2015.
20. STATA. User's Guide Release 13. StataCorp, 2013.
21. Venables W.N., Smith D.M. and the R Core Team. An Introduction to R. Notes on R: A Programming Environment for Data Analysis and Graphics. Version 3.4.0, 2017.
