	Расчет равновесий в аналитической химии

Главная » Программы » Обратная задача равновесия » DALSFEK

Программа DALSFEK для вычислений констант устойчивости

Программа впервые в отечественной литературе была описана в книге Хартли Ф., Бёргес К., Олкок Р.. Равновесия в растворах. – М.: Мир, 1980., в ней же приведены текст программы на языке Фортран, инструкция по подготовке данных, примеры использования.

Мы приводим выдержку из этой книги, которая касается программы DALSFEK (скачать в DJVU – 1.2МБ).

Программа DALSFEK для вычисления констант устойчивости по нелинейному методу наименьших квадратов с оценкой поправок

1. Введение

В этом приложении обсуждается программа DALSFEK вычислений по нелинейному методу наименьших квадратов [1] применительно к анализу данных, полученных с целью определения констант устойчивости [2, 3]. Описание этой программы полезно и для тех, кто хотел бы ее использовать, и для тех, кто создает свою программу. В приложении обсуждаются применяемые нами методы выполнения определенных задач, связанных с данной программой. Наша программа включает один из лучших алгоритмов, описанных в литературе [4] для обеспечения сходимости (см. разд. Б.7), однако не исключается возможность адаптации для использования других алгоритмов, поскольку программа написана таким образом, что в этом случае требуется лишь замена подпрограммы.

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

Блок-схема представлена на рис.; текст программы дается в конце приложения.

Начало

Шаг 1 Программа MASTER. Считывает управляющую переменную и определяет равновесия. Считывает начальные оценки параметров. Вызывает программу CYCLE

Шаг 2. Программа CYCLE Считывает массивы параметров констант и устанавливает массивы уточняемых параметров. Для каждой серии экспериментальных ванных считывает аналитические концентрации и вызывает каждый раз CONCENTRATION-программу

Шаг 3 Программа CONCENTRATION на основании определения равновесия и текущих значений параметров, находит концентрации частиц и хранит их. Вызывает программу OBSERVABLE

Шаг 4 Программа OBSERVABLE считывает экспериментальные значения зависимых переменных для ванного эксперимента и запоминает их. Находит расчетные значения концентраций частиц. Вычисляет разности (ошибки) и, следовательно, накапливает суммы квадратов ошибок. Возвращает управление программе CYCLE

Шаг 5 Программа CYCLE вычисляет и запоминает элементы матрицы В. Дает приращение для K и вызывает программу CONCENTRATION

Шаг 6 Программа CONCENTRATION находит концентрации частиц на основании измененных значений констант равновесия. Вызывает программу OBSERVABLE

Шаг 7 Программа OBSERVABLE вычисляет изменения рассчитанных значений зависимых переменных от концентраций частиц. Возвращает управление к программе CYCLE

Шаг 8 Программа CYCLE после повторения цикла для каждого эксперимента вычисляет и нормализует BWB. Вызывает программу LEAST-SQUARE

Шаг 9 Программа LEAST-SQUARE Вычисляет собстственные векторы и собственные значения BWB. Вычисляет поправки к параметрам. Пересчитывает сумму квадратов ошибок, вызывая программу CONCENTRATION

Шаг 10 Программа CONCENTRACION находит концентрации частиц, используя уточненные значения параметров. Вызывает программу OBSERVABLE

Шаг 11 Программа OBSERVABLE пересчитывает значения зависимых переменных, используя уточненные значения параметров. Пересчитывает разности и, следовательно, сумму квадратов ошибок. Возвращает управление программе LEAST-SQUARE

Шаг 12 Программа LEAST-SQUARE использует алгоритм Марквардта для решения вопроса, принять ли новые значения параметров либо следует внести новые поправки. Каждый раз при вычислении суммы квадратов ошибок вызывает программу CONCENTRATION. Передает управление программе CYCLE

Шаг 13 Программа CYCLE устанавливает, существует ли сходимость. Если существует, то вызывает программу ERROR В противном случае начинает новую итерацию с шага 2, используя новые значения параметров

Шаг 14 Программа ERROR вычисляет стандартные отклонения значений параметров. Вычисляет коэффициенты корреляции параметров

остановка

Рис. 1. Блок-схема программы DALSFEK для вычислений констант устойчивости по нелинейному методу наименьших квадратов.

2. Общее описание

2.1. Основная программа (MASTER)

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

2.2. Программа цикла (CYCLE)

Эта программа управляет применением дополнительной памяти в виде записи на магнитной ленте. В этом случае обеспечивается систематическое запоминание общих (аналитических) концентраций компонентов, измеренных параметров и концентрации каждой частицы по завершении цикла обработки каждой серии экспериментальных данных* (результатов одного потенциометрического титрования или одной серии измерения светопоглощения). Здесь вычисляется B-матрица дифференциалов (см. гл. 5), а эти вычисления необходимы для получения B^TWB. Следовало бы использовать эти элементы, полученные аналитически, что, однако, невозможно при определении констант устойчивости. DALSFEK вычисляет приближения разностей в частных производных согласно соотношению

∂o_i / ∂K_j = Δo_i / ΔK_j → (1)

---------
* Если при сравнении с числом значащих цифр уже получается ω_i-1<<1, то продолжайте процедуру по позициям (2) и (3), минуя сравнение с S(ω_i-i/V).

Хотя функция Δ, представляющая собой функцию приращений константы устойчивости, и является переменной величиной, однако нами показано, что ее значение 0,001 (0,1%) удовлетворительно. Перед вызовом итеративной программы вычисления наименьших квадратов нормализуют матрицу B^TWB к такому виду, чтобы диагональные элементы приняли значение единицы.

2.3. Программа для вычисления концентраций (CONCENTRATION)

В начале каждой итерации требуется определить концентрации частиц в системе. (Эта процедура необходима также при вычислении приближений разностей дифференциалов для констант устойчивости.) DALSFEK достигает этого при помощи метода итерации Гаусса — Ньютона на некоторых «предполагаемых» концентрациях, используя текущие значения констант устойчивости и данные по материальному балансу (общие концентрации). Подобные процедуры включены в ранее опубликованные программы для расчета равновесий в растворах [5,6].

Рассмотрим m равновесий в растворе, в которых участвуют n частиц. Тогда для этих равновесий можно написать m уравнений вида

ln K_i = a_i1 ln c₁ + ... + a_in ln c_n → (2)

где a_in > 0 для продуктов реакций и <0 для реагирующих веществ. В матричной записи

ln К= A ln С → (3)

Рассмотрим также уравнения материального баланса вида

t_j = q_j1c₁ + ... + q_jnC_n → (4)

где t_j — общая (аналитическая) концентрация компонента j. Величина q_ij принимает целочисленные значения и характеризует число атомов компонента j в частице i. В матричной записи

T = QC → (5)

Поскольку значения A и Q вводятся при определении равновесия, то, зная T и текущие значения K, можно уточнять набор предварительных оценок С до тех пор, пока не будут получены равенства (3) и (5).

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

2.4. Программа для обработки спектрофотометрических данных

Спектрофотометрия является одним из двух наиболее часто применяемых методов исследования равновесий в растворах, поэтому DALSFEK включает программу для обработки спектророфотометрических данных. На первой итерации считываются и запоминаются измеренные значения светопоглощения для каждого опыта, так же как и оценки молярных коэффициентов погашения каждой частицы при каждой длине волны (только при первом входе). На каждой из последующих итераций оценивают рассчитанные значения светопоглощения по концентрациям частиц, определенным в программе CONCENTRATION, предполагая соблюдение закона Бера

A^λ = ∑_{(i = 1 ÷ n)} ε_i^λ С_i l → (6)

где A^λ — светопоглощение при длине волны λ, ε_i^λ — молярный коэффициент погашения частицы i при длине волны λ, c_i — концентрация, l — толщина кюветы и n — общее число частиц в растворе. Такой подход позволяет определить расхождения между наблюдаемыми и рассчитанными значениями светопоглощения и, следовательно, сумму квадратов ошибок.

2.5. Программа для потенциометрических данных

Программа DALSFEK также включает подпрограмму для обработки данных, полученных потенциометрическим методом, вторым наиболее часто применяемым методом определения констант равновесия. Программа считывает и запоминает измеренные значения э.д. с., а также рассчитывает значения э. д. с. по уравнению Нернста на основании вычисленных концентраций. В программу для таких вычислений, конечно же, необходимо ввести информацию о том, какую из частиц характеризуют полученные потенциометрические данные.

2.6. Программы для оценки дисперсии (VARIANCE)

Эти программы могут быть включены, если нужно вычислить веса на основании оценок дисперсий разностей. Для проведения таких вычислений используют правила распространения ошибок [7], по подпрограммы требуют оценки ошибок экспериментальных данных: светопоглощения, э.д.с. и общих концентраций. Заметим, что дисперсии зависят от значений параметров, поэтому дисперсии необходимо вычислять на каждой итерации.

2.7. Программа алгоритма метода наименьших квадратов (LEAST-SQUARE)

Приведенная в этом приложении программа (DALSFEK) использует алгоритм Марквардта [4] для вычисления поправок к параметрам, применяемым в этой итерации. Этот алгоритм оказался очень удачным для решения наших задач, тем не менее можно было использовать и другие алгоритмы, например алгоритм Флетчера — Пауэлла. Целью последующей процедуры является получение собственных векторов и собственных значений B^TWB (квадратной симметричной матрицы) и вычисления, вклада каждого собственного значения в поправочный вектор параметров. Затем вводятся соответствующие поправки и вызываются подпрограммы для вычисления концентраций и зависимых переменных с целью пересчета суммы квадратов ошибок и последующего сравнения с ранее полученным значением. После этого применяют алгоритм Марквардта (разд 5.7) для решения вопроса, следует ли ослабить («damp») поправки или необходимо ввести новые значения параметров для выполнения следующей итерации. Ниже приведено краткое изложение применяемой стратегии.

Мы хотим, чтобы на i-й итерации значения ω_i (см. разд. 5.7) были такими, чтобы соблюдалось следующее условие: S_i+1 < S_i, где S — сумма квадратов ошибок. Допустим, ω принимает значение ω_i-1 из предыдущей итерации (первоначально задастся произвольное значение 0,01). Пусть V представляет собой некоторую величину >1,0 (выбор опять же произвольный). Вычислите S(ω_i-1) и S(ω_i-i/V)*, и тогда
1) если S(ω_i-i/V)<=S_i, то пусть ω_i = ω_i-1/V;
2) если S(ω_i-i/V)>S_i и S(ω_i-1)<=S_i,-, то пусть ω_i = ω_i-1;
3) если S(ω_i-i/V)>S_i и S(ω_i-1)>S_i, то необходимо увеличить ω последовательным умножением на V до тех пор, пока для некоторого наименьшего значения n не будет выполняться условие S(ω_i-i/V) <= Si, после чего пусть** ω_i-i/Vⁿ.
-----------
** В тех случаях, когда параметры плохо коррелируют, возможно увеличение ω_i до неразумно больших значений. Во избежание такой ситуации вводят масштабный множитель b_i, к поправочному вектору Δ, когда угол 0; между поправочным вектором Θ и самым крутым направлением спуска (которое определяется ω_i) меньше π/4.

Обычно можно начинать с достаточно больших значений ω_i, но часто по мере приближения к минимуму ω_i равномерно уменьшается.

2.8. Программа для вычисления ошибок

Для вычисления несмещенной оценки дисперсии с единичным весом, σ², (см. уравнение (4.34)) и стандартных отклонений параметров от диагональных элементов матрицы дисперсия-ковариация (уравнение (4.33)) на последней итерации осуществляют вход в эту программу. Кроме того, рассчитываются коэффициенты корреляции и выводятся на печать для того, чтобы можно было судить о состоянии проблемы и о воспроизводимости оценок параметров.

Комментарии к программе DALSFEK

Это сокращенный вариант программы расчета констант устойчивости по нелинейному методу наименьших квадратов при большом числе экспериментальных данных. Программа составлена Р. М. Олкоком и Д. Е. Роджерсом, сотрудниками химического отделения Саутгемптонского университета. Программа близка к алгоритму Марквардта, предназначенному для получения оптимальных поправок к параметрам на каждой итерации [Marquardt D. W., J. Soc. Ind. Арp'l. Math., 11, 431 (1963)1.

(1)
На одной карте отперфорировать следующие контрольные параметры:
«NOR» — число реакций (максимальное число 6). Формат I3.
«NOEK» — число констант равновесия (максимальное число 6). Формат I3.
«NOS» — число частиц, определяющих систему (максимальное число 8). Формат I3.
«NOW» — для спектрофотометрических данных это общее число длин волн, при которых измерялось светопоглощение (не более 20). Формат I3.
«NOCP»— число неизменяемых параметров (не более 50). Включает параметры, значения которых остаются неизменными при уточнении. Формат I3.
«NEXP» — число опытов (растворов), результаты которых используются для расчета. Формат I3.
«INCR»— величина, обратная дроби, добавленной к константам равновесия при численном дифференцировании. Например, если приращение должно быть равным 0,01, то INCR=100. Формат I4.
«NPMAX» — наибольшее число итераций. Формат I3.
«NOT» — число аналитических концентраций (основные случаи) (максимальное число 8). Формат I3.
«IFW» — параметр, определяющий выбор весов. Если IEW = 1, то используются веса, равные единице. При IFW = 0 веса рассчитывают из вычисленных ошибок экспериментальных величин. Формат I3.
«IFC» — возможность вычисления систематических ошибок в аналитических концентрациях. Если IFC = 1, ошибки вычисляют. Если IFC = 0, их не вычисляют. Формат I3.
NB (число уточняемых параметров; рассчитывается как NVAR) плюс: NEXP*NOT не должно превышать 80. (В варианте ICL 1907А-программы эта величина равна 60.)
«NDV» — общее число наблюдений, получаемых за один опыт. Для спектрофотомстрических данных NDV=NOW. Формат I3.
«NOP» — общее число параметров (фиксированных и уточняемых) на стадии уточнения. Для спектрофотометрических данных NOP = NOEK + NOS*NOW. Формат I3.
«IFRITE» — выбор вида печати выходных данных. Если: IFRITE = 0, то печатаются только текущие значения параметров, но если IFRITE = 1, то печатаются также поправки к параметрам на каждом смещении. Формат I3.

(2)
тперфорировать определение равновесия. Реагирующие вещества, константы равновесия и продукты реакции кодируются: в виде целых чисел и затем перфорируются для каждой реакции. Число карт равно NOR. Например, для равновесий 2А = В + 2С и B = D + C перфорируются две карты: 1 1 1 2 3 3 и 2 2 4 3. Первыми двумя цифрами обозначены первое и второе взаимодействующие вещества соответственно (во втором равновесии нет второго реагирующего вещества), третья цифра обозначает константу равновесия, а последние две — первый w второй продукты реакции. Формат 6Ш3.

(3)
Отперфорировать оценочные (предполагаемые) значения логарифмов (основание 10) констант равновесия. Формат 6F10.4

(4)
Отперфорировать целочисленные показатели степеней, в которые необходимо возвести концентрации частиц для получения констант равновесия. Число карт равно NOEK. Например, для указанных равновесий [см. (2)] перфорируются две карты: —2 1 2 0 и 0 —1 1 1. Формат 8I3.

(5)
Отперфорировать число неизменных параметров. Для спектро-фотометрического определения предположить, что константы равновесия дают параметры 1...NOEK, а параметры молярных коэффициентов погашения NOEK+1 ... NOP. Молярные коэффициенты погашения вводят в память в виде Е(1), Е(2), ... E(NOS) для длины волны 1, Е(1), Е(2) ... E(NOS) для длинны волны 2 и т. д. Например, определяя равновесия А = В при двух длинах волн, когда молярные коэффициенты А фиксированы, отперфорировать 2 4. Формат 26I3.

(6)
Отперфорировать используемые в эксперименте (растворе) 1 аналитические концентрации компонентов. Формат 7Е10.4.

(7)
Отперфорировать число атомов данного компонента в каждой частице. Число карт равно NOT. Каждая карта содержит NOS целых чисел. Например, для системы А+2В = АВ2 должны быть две карты: 1 0 1 и 0 1 2. Формат 813. Только для этого первого эксперимента отперфорировать набор «предполагаемых» равновесных концентраций (они не должны быть очень точными). Число равно NOS. Формат 8D7.3.

(8)
Отперфорировать «IFOBS» — целое число, указывающее, какую из имеющихся подпрограмм следует вводить. Например, если для каждого эксперимента необходимо ввести как спектрофотометрические, так и потенциометрические данные, то отперфорировать 12 или 21. В первом случае спектрофотометрические данные должны предшествовать потенциометрическим, во втором случае ввод данных должен быть обратным.

(9)
Если IFOBS = l, отперфорировать длину волны, а затем молярные коэффициенты погашения (число которых равно NOS). Всего NOW карт. Формат F5.1, 6Е10.4. (Это только для первого эксперимента.)
Отперфорировать NOW значений светопоглощения. Формат 16F5.3. (Десятичная точка включена в программу.)

(10)
Отперфорировать оцененную ошибку светопоглощения (абсолютное значение) и оцененные ошибки в аналитических концентрациях, выраженные в десятичных дробях, а не в процентах. Формат 10F6.3.

(11)
Если IFOBS = 2, отперфорировать нернстовский коэффициент, указав номер частицы (номер, который приписан этой частице при определении равновесия) и номер аналитического компонента, который входит в электрод сравнения*. Формат (F7.4, 213).
--------------
* Имеется в виду, видимо, индикаторный электрод. — Прим. перев.

(12)
Отперфорировать оцененное значение ошибки в э. д. с. (абсолютное значение) и оцененные значения ошибок в аналитических концентрациях, выраженные в десятичных дробях, а не в процентах. Формат 9F6.3.

(13)
Затем продолжать вводить массив данных в следующем порядке: аналитические концентрации, ошибки, экспериментальные данные (светопоглощение и/или э.д.с.); общее число карт для каждого эксперимента равно NEXP.

(14)
И наконец, отперфорировать следующую карту, в которой указаны величины S (контрольная переменная): CORREC — наибольшая разрешенная поправка, CONV — критерий сходимости и SCALE. Все эти величины выражены в долях параметров. Формат I3, 3Е12.4.

Литература

1. Alcock R. М., Hartley F. R., Rogers D. Е., J. Chem. Soc. (Dalton), 115 (1978).
2. Alcock R. М., Hartley F. R., Rogers D. E., Wagner J. L., J. Chem. Soc. (Dalton), 2189 (1975).
3. Alcock R. М., Hartley F. R., Rogers D. E„ Wagner J. L., J. Chem. Soc. (Dalton), 2194 (1975).
4. Marquardt D. WJ. Soc. Ind. Appl. Maths., 11, 431 (1963).
5. Perrin D. D., Sayce I. G., Talanta, 14, 833 (1967).
6. Ting-Po I., Nancollas G. HAnal. Chem., 44, 1940 (1972).
7. Lansbury R. C„ Price V. £., Smeeth A. G„ J. Chem. Soc. (A), 1896 (1965).