Программа для расчета сложных равновесий при "плохих" приближениях Текст научной статьи по специальности «Математика»

С. М. Кузнецов, В. И. Голованов Программа для расчета сложных равновесий при "плохих" приближениях

Методами компьютерной графики построены и изучены, на ряде примеров, поверхности функций МНК-рассогласования SS(0), где 0 = || ln K1, ln K2, ..., ln Ks || - вектор логарифмов констант образования. Обнаружены общие закономерности в топологии поверхности SS(0). Показано, что топологические особенности связаны со свойствами собственных векторов(СВ) и собственных чисел(СЧ) матрицы || д2SS / dlnKj dlnKj ||. Предложен алгоритм корректировки движения к глобальному минимуму SS(0) в дополнение к той или иной известной процедуре отыскания вектора 0 = arg min SS(0). Суть корректировки заключается в том, что в промежутках между шагами итерирования вычисляются и анализируются СЧ и СВ. Если обнаруживается пологий участок поверхности по признаку СЧ, ~ 0, то текущий вектор констант 0k получает приращение в направлении СВ, до тех пор, пока не будет найден выход к глобальному минимуму. Это позволяет задавать в качестве первых приближений практически любые значения 0О, что глобально решает проблему "плохих" приближений при параметрической идентификации модели сложного равновесия. Алгоритм реализован как структурный элемент, дополняющий итерационную процедуру Марквардта, в доступной через Internet программе.

Центральной проблемой моделирования сложных равновесий на ЭВМ с целью ее параметрической идентификации является задание начального набора параметров 0 = || 0Ь 02, ..., 0s ||. Стартовые значения параметров затем уточняют последовательным приближением к искомой оценке 0 = || 01 , 02 , ..., 0s ||, исходя из выбранного критерия аппроксимации. Чаще других применяется принцип метода наименьших квадратов. МНК-оценкой параметров является вектор

0* = arg min SS(0). (1)

Сумму квадратов отклонений SS, иначе функцию рассогласования, записывают матричным произведением

SS = OeT WOe, (2)

где O0 - вектор рассогласования (невязок) в точке 0 параметрического пространства, W - весовая матрица. В свою очередь O0 задан разностью между экспериментально наблюдаемым (у) и вычисленным по уравнению функции отклика (п0) векторами:

O = у - П0. (3)

Так как в (2) является сложной нелинейной функцией относительно 0, для поиска минимума применяют тот или иной метод оптимизаций 1]. Обычно используют метод Ньютона или его модификации (метод Гаусса-Ньютона, метод Левенберга-Марквардта и др.), основанные на итерационном выражении[2]:

0к+1 = 0к - - И^О]"1 ^^О или (4а)

0к+1 = 0к - (И1 + И2)-1 ц. (4б)

В рекуррентных формулах (4): ^ - матрица первых частных производных вектора-функции п по вектору 0; Ип - матрица вторых производных; О -вектор рассогласования (см. (3)); И1 - матрица сумм кросс-произведений первых производных; И2 = - Ипт^О; ц - градиент 88; к - номер итерации (все матрицы рассчитываются в точке 0к). Матрицу Гессе функции рассогласования находят как

И = И1 + И2. (5)

Проблема «плохих» приближений связана с плохой обусловленностью матрицы И в некоторых далеко отстоящих от глобального минимума областях параметрического пространства. Вырождение матрицы Гессе приводит к невозможности ее обращения, и вектор поправок Д0 становится неопределенным, т.е. бесконечно большим. Вырождение И связывают со сложной топологией поверхности 88(0). Областями плохой обусловленности являются, прежде всего, «овраги». Предполагают [3, с.88] сложную искривленную форму контуров оврагов. Только при попадании первых приближений в область глобального минимума сходимость итерационных процедур гарантирована. Попадание, на пути к минимуму, в узкий, а главное, длинный овраг может потребовать очень большого числа итераций или даже привести к остановке процесса.

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

Моделирование функций рассогласования

-г, N х 1

В эксперименте элементы вектора у получают при реализации некоторого плана эксперимента Е№ в котором, по крайней мере, одна из независимых переменных (факторов) принимает N уровней значений. В качестве факторов, для простоты, выберем только общие концентрации независимых компонентов (ключевых веществ, базисных частиц). Дополнительно потребуем выполнения принципа постоянства свойств среды (постоянства ионной

силы, в частности), на фоне которой протекают изучаемые реакции. Вектор-функция откликов у полностью задана, если:

1 ) задана структура, а также параметры 0 модели равновесия;

2) задан план EN ;

3) определена с точностью до неизвестных параметров зависимость свойства системы от ее равновесного состава.

Физико-химическую модель равновесия запишем, как это принято при описании сложного равновесия [4], матрицей состава (табл. 1).

В компактной матричной форме система уравнений реакций образования всех частиц равновесной системы запишем как:

УБТ = А или |№|| ||Б^Г = ||А1||, 1 = Т7 ] = Т~т , (6)

где V - матрица п х т стехиометрических коэффициентов в табл.1,

Б - вектор-строка базисных частиц порядка т, А - вектор-столбец порядка п всех химических форм, отвечающих за изменение термодинамического потенциала О системы при ее переходе от начального к равновесному состоянию.

Таблица 1

Обобщенная матрица состава системы со сложным равновесием.

1п К Б: Компоненты (Б^ Б2 ... Бт Частицы (А1)

0 1 0 0 0 А1 = Б1

0 0 1 0 0 А2 = Б2

0 0 0 0 1 Ат = Бт

Ьп Кт+1 'Ут+и Vm+1,2 Vm+1,m Ат+1

Ьп Кт+2 Vm+2,1 Vm+2,2 Vm+2,m Ат+2

1п Кп Vn,1 Vn,2 Ап

п - т = г > 0

Систему уравнений материального баланса запишем в виде: vTexp а = с, (7)

где ехр а = || ехр а1 || - вектор-столбец показательных функций, причем элементы вектора а являются логарифмами равновесных концентраций частиц а1 = 1п [А1]. В свою очередь, обобщенный закон действия масс для сложного равновесия удобно представить как

а = р + vb . (8)

1 Здесь рассмотрим только однооткликовую модель.

132

Компоненты векторов в и Ь - это ^ = 1п К! и Ь = 1п [Б^, соответственно. Вектор с в уравнении (7) состоит из общих концентраций базисных частиц. Набор из N векторов общих концентраций образует спектр матрицы плана .

Решение нелинейной системы уравнений (7), (8) относительно Ь получают методом Ньютона-Гаусса, который опишем рекуррентной формулой [5]

/ т Ч и Г деп х (с - V ехр ар) = Ьр -

Ьр+1 = Ьр ■

--Ч, т .-1-1 Г-ч -1-1

Э(с - V ехр ар)

дЬр

1Е-

дЬр

ер, (9)

где ер - текущий вектор невязок общих концентраций. Каждый (к, 1)-й элемент матрицы шхш первых производных Эер/ЭЬр (матрицы Якоби) в уравнении (9) находят по формуле

Эе п

~ЭЬ = ^ v1l ехр аг . (10)

7=1

Найденные для точки плана Едг равновесные концентрации позволяют рассчитать П х вектор значений функций отклика.

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

При моделировании поверхности 88(0) рассмотрим два наиболее распространенных варианта формулирования отклика:

п

П = 1п [А] и п = Ъц [А] . (11)

7=1

Первый относится, например, к ионометрии, а второй - к методу растворимости. Также первая функция может быть косвенным результатом по-тенциометрии, а вторая - спектрофотометрии.

Исходя из выражений (8) и (11), очевидно, что П является функцией набора уточняемых параметров 0 4 х 1 е в пх \ т.е. П = П(0). Тогда уравнение (2) перепишем, заменив для простоты матрицу W единичной2:

88(0) = (у - П(0))т(у - П(0)). (12)

Вектор у в уравнении (12) будем вычислять как

у = п(0ист) + е, (13)

где 0ист - истинное значение вектора || 1п Кь 1п К2, ..., 1п К ||, заданное в вычислительном эксперименте, е - вектор погрешностей, который получали с использованием генератора нормально распределенных случайных чисел. Уравнения (9) - (13) позволяют генерировать значения функции рассогласо-

2 Это равносильно распространенному допущению о равной точности измерений откликов в разных точках плана эксперимента.

Рис. 1. Функция рассогласования при моделирования титрования двухосновной кислоты щелочью: В = \\Н+, В^Ц; п = рН; в = || ^КИЕ ^КН2В \ \ = \ \ 6, 12 \ \. А - 3Б диаграмма; Б - диаграмма линий постоянных уровней

вания 88(в) как результат решения прямой задачи равновесия на множестве значений в е Я Для этого точки параметрического пространства в, на котором отображается задавали в виде многомерной сетки, покрывающей некоторую область (подпространство), в центре которой находится точка вист.

На рис. 1А приведена функция рассогласования для простейшего из рассмотренных нами примеров. Рис. 1Б представляет «вид сверху» этой функции с линиями постоянных уровней 88. При рассмотрении примера на рис. 1 и других примеров был установлен ряд общих закономерностей в топологии 88(в).

1. Поверхность 88(в) характеризуется тремя основными элементами: плато на периферии объемлющего пространства, оврагами, котловиной глобального минимума, к которой сходятся овраги.

2. Величина перепада "овраг - глобальный минимум" зависит от вариации погрешностей £, т.е. от качества экспериментальных данных. Максимальный перепад наблюдали, когда вектор ошибок задавали равным нулю. При этом, разумеется, 88(вист)=0. Когда погрешности довольно велики, котловина сливается с оврагом самого нижнего уровня.

3. Число оврагов равно числу оцениваемых параметров, увеличенному на единицу3, т.е. 5+1. Причем 5 оврагов проходят вдоль каждой из 5 коор-

3 При условии, что все частицы, для которых оцениваются параметры, значимы.

динатных осей параметрического пространства. Дополнительный (4+1)-й овраг ориентирован по «диагонали» подпространства. Угол между направлением дополнительного оврага и координатными осями изменяется при изменении матрицы стехиометрических коэффициентов. Число плато (или подпространств с постоянным значением) зависит от параметров системы и вида функции п, но всегда не превышает 4+1.

4. При движении от периферии к центру, т.е. к глобальному минимуму, плато обрываются каскадами оврагов, которые «втекают» в котловину. Таким образом, движение по поверхности функции рассогласования представляет собой преодоление череды пологих и крутых участков. Перечисленные выше правила, по-видимому, носят общий характер и справедливы не только для Я3, но и для пространств высших порядков. Отмеченные наблюдения легли в основу нового алгоритма движения к минимуму, который мы обсудим позже. Особо обратим внимание на правило п. 3 перечня, согласно которому овраги линейные вдоль главной оси. Возможно, это обусловлено спецификой задач химического равновесия. Линейность оврагов существенно упрощает расчеты движения к минимуму и исключает вмешательство исследователя в расчетные процедуры. Отметим также, что плато -особый овраг с многими главными осями.

Тактика движения по поверхности должна быть различной в зависимости от характера ее рельефа на том или ином шаге поиска оптимума. Алгоритмы поиска, основанные на расчете градиентных характеристик, не имеют эффекта на «плоских» участках 88. На них нет градиента. Отсюда, при попадании в овраг следует отказаться от использования процедур (4), а использовать процедуру движения только вдоль большей оси оврага.

Для записи алгоритма в терминах программирования нужно установить аналитические признаки оврагов и признаки их главных осей. Эти признаки удается обнаружить при объяснении топологии функции рассогласования.

Объяснение топологии функции рассогласования

Топологию функции рассогласования объясним, исходя из ее локального описания квадратичной формой. Разложим функцию 88(0) в ряд Тейлора по степеням приращений Д0 в окрестности некоторой точки 0р и ограничимся членами второго порядка:

88(0р+Д0) - 88(0р) + £(0/ Д0 + 1/2 Д0т Н Д0, (14)

где £(0^) - вектор-градиент 88, Н - гессиан 88.

В стационарных точках 0р поверхности 88(0), т.е. в точках экстремумов, перегибов или вырожденностей, градиент равен нулю. Тогда вместо уравнения (1 4) имеем

Д88 - V2 Д0 т Н Д0. (15)

Уравнение (15) является квадратичной формой, описывающей эллиптический параболоид, центрированный в точке 0Р. Заметим, что дифференцирование

(1 4) по в и приравнивание результата нулю приводит, в частности, к общим выражениям (4) для расчета центра эллиптического параболоида, заданного его точкой вр4 Исходя из геометрической трактовки квадратичной формы, контуры равных уровней 8 8 (в), представляют собой различные типы эллипсов (в общем случае - это эллипсоиды). Из теории квадратичной формы следует, что большей оси эллипсоида соответствует меньшее по величине собственное число X матрицы Гессе. Когда какое-либо Хг- ~ 0, гессиан становится плохо обусловленным, а соответствующая ему квадратичная поверхность приобретает элементы вырожденности в виде оврагов, плато, подпространств постоянных уровней и т. п. Таким образом, обнаружение собственных чисел, близких нулю, равносильно обнаружению оврагов (или других вырожденностей поверхности).

Ориентация в параметрическом пространстве эллипсоидов равных уровней характеризуется собственными векторами матрицы Н. Отсюда собственный вектор, соответствующий СЧ = 0, указывает направление вдоль оврага. Движение в одном из двух направлений вектора приводит к глобальному минимуму. Поскольку заранее неизвестно, какое из двух направлений ведет к цели, приходится алгоритмически исследовать оба варианта (если, конечно, первый из выбранных не привел к положительному результату).

Здесь следует отметить два важных момента. Аналитическим признаком плато является гессиан, все элементы которого равны нулю (нулевой гессиан). Существует несколько равноправных направлений движения по его поверхности. Однако это не является непреодолимым препятствием. Оказывается, что выбрав любое направление из точки Р на рис. 1Б, мы попадаем в один из оврагов, который в конце концов выводит нас к глобальному минимуму.

Второе замечание относится к определенности гессиана. Гессиан Н квадратичной формы (15) не всегда является положительно полуопределенным, т.е. таким, что все X! > 0. Он может быть и отрицательно полуопределенным и неопределенным вовсе. Это, между прочим, означает, что формула (4) описывает движение к стационарной точке с g(в) = 0, которой может быть и максимум. Отсюда процесс поиска (4) может увеличивать критериальную функцию 88. Гессиан отрицательно полуопределен, т.е. все X! < 0, в промежутке между плато, являющимся максимумом, и дном оврага, являющимся минимумом, до точки перегиба (на верхней стенке оврага). Следствием этого является то, что необходимо применять алгоритмы, модифицирующие, при необходимости, матрицу Н так, чтобы обеспечить ее положительную определенность.

Одним из методов, устраняющих проблему положительной неопределенности гессиана (5), является метод Марквардта, приспособленный авторами программы БЛЬ8РЕК[3, 7] для решения задач моделирования

4 Вернее, точкой вр с ее градиентными характеристиками.

136

равновесий. Здесь, в отличии от (4), текущий вектор параметров рассчитывается по формуле

0к+1 = 0к - VT(Л+ц ^^ 8(0р)т, (16)

где V - квадратная матрица порядка 8, столбцы которой - это собственные вектора ^ Л - диагональная матрица с собственными векторами H по диагонали; I - единичная матрица; ц - такой скаляр, чтобы все (А,; + ц) > 0. Варьирование ц позволяет изменять вектор поправки от направления, рассчитанного в предположении квадратичной аппроксимации 88, до направления антиградиента - направления наискорейшего спуска. Другие методы по другому решают эту проблему. Например, метод Гаусса-Ньютона пренебрегает поправками на вторые производные в матрице ^ (4Б). Тогда гессиан равен матрице ^ = JTWJ, которая по определению положительно полуопределена. Применение того или иного метода, однако, не решает проблему движения к минимуму в условиях «овражности» минимизируемой функции. Поправки, полученные по выбранному алгоритму, не приводят к уменьшению 88 в овраге или на плато, что приводит к прерыванию вычислений или зацикливанию. В любом случае «стандартные» методы требуют задания достаточно точных начальных приближений. Другими словами, стартовая точка параметрического пространства 0О должна лежать в области котловины. Следует отметить, что в этой области все методы так или иначе справляются с поставленной задачей. Нами проверялись следующие методы: Девидона-Флетчера-Пауэлла (редакция программы AUTOEQUIL[8]), деформируемого многогранника^], Марквардта(программа DALSFEK[3, 7]), Ньютона-Гаусса (редакция программы ^ЮТ[9]). Все они показали примерно одинаковую ско-

5

рость расчета и сходимость .

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

Алгоритм решения обратной задачи равновесия и его исследование

(рис. 2)

1. Выбрать начальное приближение (для универсальности можно рекомендовать выбирать значения констант заведомо заниженными, так, как будто соответствующие им равновесия не вносят вклада в функцию отклика).

2. «Стандартным» методом нелинейного МНК решить обратную задачу равновесия. При выборе метода руководствоваться рассуждениями из предыдущей главы. Запомнить текущее значение 88.

5 При достаточно высокой мощности вычислительной техники.

137

3. По формуле (5) найти гессиан Н и его спектральное разложение (Собственные числа и собственные векторы находят стандартной процедурой спектрального разложения H = V ЛV [11]. В случае нулевого гессиана применяется тривиальное соотношение V = I, Х=0, где I - единичная, 0 -нулевая матрицы). Зафиксировать все собственные числа, равные нулю.

4. Если нет СЧ, равных нулю, идти на п. 9.

5. Выбрать очередное СЧ (равное нулю) и присвоить вектору поправок значение его собственного вектора.

6. Изменять вектор параметров на величину вектора поправок (возможно с некоторым множителем), пока не уменьшится значение 88 или значение какой-либо константы не превысит критического значения.

7. Если значение 88 уменьшилось, запомнить текущие значения параметров и перейти на п. 2.

8. Иначе, если направление вектора поправок не изменялось, изменить направление на обратное и перейти на п. 6. В противном случае - на п. 4.

9. Рассчитать статистические параметры и закончить работу алгоритма. Наш алгоритм существенно отличается от давно известного «метода оврагов» [6, с. 209], в котором также используется овражность рельефа целевой функции. Алгоритм метода оврагов рассчитан на то, чтобы пройти вдоль оврага и выйти в котловину около минимума. Этим методом удается находить минимумы функций 5-10 переменных. Отмечают [6], что «метод довольно капризен». Для каждой функции нужно подбирать свой овражный шаг, интерактивно наблюдать за ходом расчета и вносить коррективы.

Описанный выше алгоритм программно реализован для персональных компьютеров с операционной системой Б08^т32. В качестве «стандартного» метода нелинейного МНК использовался метод Марквардта. Демоверсию программы с примерами и файлами справки можно загрузить с WWW-страницы http://sinisha.chat.ru/products/newdalsfek/. Для тестирования и отладки программы испытывали теоретические модели равновесий различной сложности6. Испытание алгоритма на модельных системах подтвердило его работоспособность при варьировании стартовых точек параметрического пространства, уровней погрешностей в данных, числа откликов. Положительные во всех случаях, результаты испытаний указывают на общность алгоритма, общность лежащих в его основе правил и теоретических предпосылок, которые мы обсуждали выше.

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

6 Более двух десятков примеров.

Таблица 2

Результаты самосогласования констант для систем Н+-еп-М2+

(1М КС1)[10]

(Первые приближения общих констант равновесия равны -10.0)

Частицы Рассчитанные общие константы Ступенчатые константы

рассчитанные [10]

Н еп+ 10.070 ± 0.005 10.070 ± 0.005 10.05

Н2еп2+ 17.391 ± 0.006 7.321 ± 0.008 7.31

№(еп)2+ 7.69 ± 0.04 7.69 ± 0.04 7.66

№(еп)22+ 14.10 ± 0.04 6.41 ± 0.05 6.40

№(еп)з2+ 18.72 ± 0.05 4.63 ± 0.06 4.55

Мп(еп)2+ 2.86 ± 0.04 2.86 ± 0.04 2.73

Мп(еп)22+ 4.82 ± 0.05 1.96 ± 0.06 2.06

Мп(еп)з2+ 5.90 ± 0.07 1.08 ± 0.09 0.88

Бе(еп)2+ 4.30 ± 0.05 4.30 ± 0.05 4.28

Бе(еп)22+ 7.61 ± 0.05 3.31 ± 0.07 3.25

Бе(еп)з2+ 9.64 ± 0.05 2.02 ± 0.08 1.99

Со(еп)2+ 5.92 ± 0.03 5.92 ± 0.03 5.89

Со(еп)22+ 10.79 ± 0.03 4.87 ± 0.04 4.83

Со(еп)32+ 13.93 ± 0.03 3.14 ± 0.04 3.10

тальных данных Я.Бъеррума, представленных объединенной матрицей плана, с использованием также объединенной матрицы состава. Таким образом, матрица состава содержала 14 уточняемых констант образования: двух констант протонирования еп и 12 констант ступенчатого образования комплексов. Первые приближения логарифмов всех общих констант образования были заданы равными = -10. Разумеется не потому, что посчитали их лучшими, а для иллюстрации возможностей программы. Найденные при этом самосогласованные вектора общих (полных) и частных ступенчатых констант можно сравнить с константами Я. Бъеррума в табл. 2. Значения ступенчатых констант Я. Бъеррума "накрываются" 2а пределами найденных нами оценок этих констант практически во всех случаях, кроме, пожалуй, системы с Мп2+. Тем не менее, найденные нами константы более предпочтительны для использования при решении прямых задач моделирования, поскольку они прошли самосогласование и поэтому их значения являются лучшими точечными оценками действительных значений.

1. Задать началь ные приближения 6 о

2. Решить обратную зад 0 = 31 Положить 88 ачу "стандартным" методом g min SS(6). >1 = SS(6 ); 6тек = 6

3. Вычислить У, Л матрицы Гессе Н=УТЛ У. Создать массив Л с номерами СЧ = 0. Положить: i = 1.

5. Присвоить вектору поправок (Э) СВ, соответствующий СЧ с индексом ^ Э = Уи .

Установить направление = 1

9. Рассчитать статистические параметры

СТОП

втек = 0 тек +

k - корректирующий множитель

882 = 88(0тек)

6тек = 6 ;

dir = -1 (меняем направление)

Рис. 2. Блок-схема алгоритма решения обратной задачи моделирования равновесий

Список литературы

1. Химельблау Д. Прикладное нелинейное программирование. М.: Мир, 1975. 534 с.

2. Демиденко Е.З. Линейная и нелинейная регрессии. М.: Финансы и статистика, 1981. 302 с.

3. Хартли Р., Бергес К., Олкок Р. Равновесия в растворах. М.: Мир, 1983. 360с.

4. Бек М., Надьпал И. Исследование комплексообразования новейшими методами. М.: Мир, 1989. 413 с.

5. Васильев В.П., Бородин В. А., Козловский Е.В. Применение ЭВМ в химико-аналитических расчетах. М.: Высш. шк., 1993. 112 с.

6. Калинкин Н.Н. Численные методы. М.: Наука, 1978. 512 с.

7. Alcock R.M., Hartley F.R., Rogers D.E. // J. Chem. Soc (Dalton). 1978. V. 2. P. 115.

8. Кирьянов Ю.А., Николаева Л.С., Евсеев А.М. Автоматизированная система математического моделирования химических равновесий с учетом кинетики баланса масс (AUTOEQUIL) // Мат. моделирование хим. равновесий. М.: Изд-во МГУ, 1988. С. 146-186.

9. Бугаевский А.А., Холин Ю.В., Коняев Д.С. //Журн. физ. химии. 1993. Т. 38, № 2. С. 350.

10. Бъеррум Я. Образование аминов металлов в водном растворе. М.: Ивдстр. лит., 1961. 308 с.

11. Сборник научных программ на ФОРТРАНЕ: Пер. с англ.; Под ред. С.Я. Виленкина. М.: Статистика, 1974. Вып. 2. 224 с.