Алгоритмы

Эта страница описывает основные численные алгоритмы, используемые внутри PDIndexer. Это перенесённая и реорганизованная версия пояснительного PDF (PDIndexerAlgorithm.pdf), который раньше поставлялся вместе с дистрибутивом. Цель — передать что минимизируется и как это решается, а не дать полную математическую строгость.

Рассматриваются три темы:

1. Уточнение параметров решётки — линейный метод наименьших квадратов
2. Аппроксимация пиков — нелинейный метод наименьших квадратов по методу Марквардта и профильные функции
3. Вывод кубического сплайна — кривая фона

Теорию уравнений состояния (EOS) см. в разделе Уравнения состояния.

---

Уточнение параметров решётки

Обобщённый линейный метод наименьших квадратов

Пусть даны \( n \) наборов наблюдений \( (X^1_1, X^2_1, \dots, X^m_1, Z_1),\ (X^1_2, X^2_2, \dots, X^m_2, Z_2),\ \dots,\ (X^1_n, X^2_n, \dots, X^m_n, Z_n) \). Подгонка линейного уравнения наблюдений

\[Z = a^1 X^1 + a^2 X^2 + \dots + a^m X^m\]

по \( m \) параметрам \( (a^1, a^2, \dots, a^m) \) достигается минимизацией суммы квадратов остатков. В матричной форме

\[ \mathbf{a}=\begin{pmatrix}a^1\\ \vdots\\ a^m\end{pmatrix},\quad X=\begin{pmatrix}X^1_1 & X^2_1 & \cdots & X^m_1\\ X^1_2 & X^2_2 & \cdots & X^m_2\\ \vdots & & & \vdots\\ X^1_n & X^2_n & \cdots & X^m_n\end{pmatrix},\quad Y=\begin{pmatrix}Y_1\\ \vdots\\ Y_n\end{pmatrix},\quad W=\begin{pmatrix}W_1 & 0 & \cdots & 0\\ 0 & W_2 & & \\ \vdots & & \ddots & \\ 0 & & & W_n\end{pmatrix} \]

где \( W \) — диагональная матрица весов. Минимизация взвешенной суммы квадратов

\[\Delta^2 = (X\mathbf{a}-Y)^{\mathsf{T}}\,W\,(X\mathbf{a}-Y)\]

путём приравнивания к нулю её производной по \( \mathbf{a} \),

\[\delta\Delta^2 = 2\,\delta\mathbf{a}^{\mathsf{T}}\left(X^{\mathsf{T}}W X\,\mathbf{a} - X^{\mathsf{T}}W Y\right) = 0,\]

даёт решение

\[\mathbf{a} = (X^{\mathsf{T}}W X)^{-1}\,X^{\mathsf{T}}W Y.\]

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

Для уточнения параметров решётки уравнение наблюдений зависит от кристаллической системы, но в наиболее общем (триклинном) случае соотношение между межплоскостным расстоянием \( d \) и индексами \( (h,k,l) \),

\[\left(\frac{1}{d}\right)^2 = h^2 a^{*2} + k^2 b^{*2} + l^2 c^{*2} + 2kl\,b^*c^*\cos\alpha^* + 2lh\,c^*a^*\cos\beta^* + 2hk\,a^*b^*\cos\gamma^*,\]

рассматривается как линейная модель:

\[ \mathbf{a}=\begin{pmatrix}a^*a^*\\ b^*b^*\\ c^*c^*\\ 2b^*c^*\cos\alpha^*\\ 2c^*a^*\cos\beta^*\\ 2a^*b^*\cos\gamma^*\end{pmatrix},\quad X=\begin{pmatrix}h_1h_1 & k_1k_1 & l_1l_1 & k_1l_1 & l_1h_1 & h_1k_1\\ \vdots & & & & & \vdots\\ h_nh_n & k_nk_n & l_nl_n & k_nl_n & l_nh_n & h_nk_n\end{pmatrix},\quad Y=\begin{pmatrix}(1/d_1)^2\\ \vdots\\ (1/d_n)^2\end{pmatrix} \]

где \( a^*, b^*, \dots \) — обратные параметры решётки. Решение этой системы описанным выше линейным методом наименьших квадратов даёт компоненты обратного метрического тензора, из которых затем получаются параметры решётки.

Выбор весов

Вес зависит от ошибки. Если предположить, что ошибка присутствует только в угле дифракции \( 2\theta \), то отклик \( G = (1/d)^2 = 4\sin^2\theta/\lambda^2 \) на \( \theta \) равен

\[\frac{\partial G}{\partial\theta} = \frac{4\sin(2\theta)}{\lambda^2},\]

поэтому изменение \( \delta\theta \) сдвигает \( (1/d)^2 \) на \( \dfrac{4\sin(2\theta)}{\lambda^2}\delta\theta \). Следовательно, \( 1/\sin^2(2\theta) \) (обратная величина квадрата ошибки) является подходящим весом для \( (1/d)^2 \):

\[ W=\begin{pmatrix} 1/\sin^2(2\theta_1) & 0 & \cdots & 0\\ 0 & 1/\sin^2(2\theta_2) & & \\ \vdots & & \ddots & \\ 0 & & & 1/\sin^2(2\theta_n) \end{pmatrix}. \]

Здесь \( 1/\sin^2(2\theta) \) представляет лишь отношение обратных дисперсий точек, а не их абсолютное значение, однако оптимум всё равно находится верно: в выражении \( (X^{\mathsf{T}}W X)^{-1} X^{\mathsf{T}}W Y \) множитель \( W \) встречается дважды, поэтому абсолютный масштаб сокращается.

Ошибки параметров

Ошибки (дисперсии) \( \mathbf{a} \) получаются из диагонали \( (X^{\mathsf{T}}W X)^{-1} \), однако, поскольку \( W \) была задана лишь с точностью до отношения, абсолютный масштаб нужно определять отдельно. Используя определение дисперсии,

\[N - P = \sum_{i=1}^{n}\frac{\delta_i^2}{s_i}\]

(\( N \): число данных, \( P \): число параметров, \( \delta_i \): остаток \( i \)-го значения, \( s_i \): дисперсия \( i \)-го значения), масштаб дисперсии фиксируется по полученным параметрам как

\[s = \frac{\sum_{i=1}^{n}(\delta_i^2 w_i)}{N - P},\]

а её квадратный корень является ошибкой. Это ошибка обратных параметров решётки; для пересчёта в ошибку параметров решётки требуется дополнительно провести распространение ошибки, что принципиально несложно.

---

Аппроксимация пиков

Метод Марквардта

PDIndexer выполняет аппроксимацию пиков методом Марквардта (Levenberg–Marquardt), нелинейной итерационной схемой, похожей на метод Ньютона. Он сочетает быструю сходимость со стабильностью и находит оптимум с достаточной точностью.

Пусть аппроксимирующая функция — \( F = F(a_1, a_2, \dots, a_m, X) \), а остаток при начальных параметрах \( \mathbf{a}^0 \) равен

\[R = \sum_{i=1}^{n}\left[Y_i - F(\mathbf{a}^0, X_i)\right]^2.\]

Построим матрицу \( m\times m \) \( \alpha \) и вектор \( \beta \) размерности \( m \) следующим образом. Умножение только диагонали на \( (1+\lambda) \) — это ключевая идея метода Марквардта, где \( \lambda \) управляет стабильностью и скоростью сходимости:

\[\alpha_{jk} = \sum_{i=1}^{n} w_i\,\frac{\partial F}{\partial a_j}\frac{\partial F}{\partial a_k}\quad(j\neq k),\qquad \alpha_{jj} = (1+\lambda)\sum_{i=1}^{n} w_i\left(\frac{\partial F}{\partial a_j}\right)^2,\]

\[\beta_j = \sum_{i=1}^{n} w_i\left[Y_i - F(\mathbf{a}^0, X_i)\right]\frac{\partial F}{\partial a_j}.\]

Параметры обновляются по формуле

\[\mathbf{a}' = \mathbf{a}^0 + \alpha^{-1}\boldsymbol{\beta}.\]

Вычисляем новый остаток \( R' \) и:

• если \( R' < R \), обновление принимается, а \( \lambda \) уменьшается (в 0,1–0,5 раза);
• если \( R' > R \), обновление отклоняется, а \( \lambda \) увеличивается (в 2–10 раз).

Повторяем до тех пор, пока изменение \( R \) не станет достаточно малым. При \( \lambda \to 0 \) метод приближается к квадратично сходящемуся методу Гаусса–Ньютона; при большом \( \lambda \) он приближается к методу наискорейшего спуска вдоль градиента остатка \( \nabla R \). Непрерывное переключение между этими двумя режимами с помощью \( \lambda \) обеспечивает стабильную и быструю сходимость.

Профильные функции

PDIndexer предлагает функцию псевдо-Фойгта (смесь функций Гаусса и Лоренца), функцию Пирсона VII (функцию плотности вероятности), а также их асимметричные расширения Split Pseudo Voigt / Split Pearson VII. Для скорости и стабильности сходимости по умолчанию используется симметричная функция псевдо-Фойгта. Все функции нормированы на единичную площадь.

Symmetric Pseudo Voigt

\[f(x,\eta,H_k)=\frac{2}{H_k\pi}\left(\eta\,\frac{1}{1+4\left(\dfrac{x}{H_k}\right)^2}+(1-\eta)\sqrt{\pi\ln 2}\cdot 2^{-4\left(\frac{x}{H_k}\right)^2}\right)\]

Split Pseudo Voigt (Toraya 1990, modified)

\[f(x,\eta_l,\eta_h,A,H_k)=\frac{2}{H_k\pi+\dfrac{H_k\pi(\eta_h-\eta_l)(Z-1)}{(1+e^{A})(\eta_h+Z(1-\eta_h))}}\left(\eta_l\,\frac{1}{1+(1+e^{-A})\left(\dfrac{x}{H_k}\right)^2}+(1-\eta_l)Z\cdot 2^{-4(1+e^{-A})\left(\frac{x}{H_k}\right)^2}\right)\quad(x<0)\]

\[f(x,\eta_l,\eta_h,A,H_k)=\frac{2}{H_k\pi+\dfrac{H_k\pi(\eta_l-\eta_h)(Z-1)}{(1+e^{-A})(\eta_h+Z(1-\eta_h))}}\left(\eta_h\,\frac{1}{1+(1+e^{A})\left(\dfrac{x}{H_k}\right)^2}+(1-\eta_h)Z\cdot 2^{-4(1+e^{A})\left(\frac{x}{H_k}\right)^2}\right)\quad(x>0)\]

Symmetric Pearson VII

\[f(x,R,H_k)=\frac{2\sqrt{2^{1/R}-1}\,\Gamma(R)}{\sqrt{\pi}\,\Gamma(R-1/2)\,H_k}\left(1+4(2^{1/R}-1)\left(\frac{x}{H_k}\right)^2\right)^{-R}\]

Split Pearson VII (Toraya 1990, modified)

\[f(x,R_l,R_h,A,H_k)=\frac{2(1+e^{A})}{H_k\sqrt{\pi}}\left(e^{A}\frac{\Gamma(R_l-1/2)}{\Gamma(R_l)\sqrt{2^{1/R_l}-1}}+\frac{\Gamma(R_h-1/2)}{\Gamma(R_h)\sqrt{2^{1/R_h}-1}}\right)^{-1}\left(1+(1+e^{-A})^2(2^{1/R_l}-1)\left(\frac{x}{H_k}\right)^2\right)^{-R_l}\quad(x<0)\]

\[f(x,R_l,R_h,A,H_k)=\frac{2(1+e^{A})}{H_k\sqrt{\pi}}\left(e^{A}\frac{\Gamma(R_l-1/2)}{\Gamma(R_l)\sqrt{2^{1/R_l}-1}}+\frac{\Gamma(R_h-1/2)}{\Gamma(R_h)\sqrt{2^{1/R_h}-1}}\right)^{-1}\left(1+(1+e^{A})^2(2^{1/R_h}-1)\left(\frac{x}{H_k}\right)^2\right)^{-R_h}\quad(x>0)\]

Первые две функции симметричны относительно \( x=0 \), тогда как split-формы меняют форму в зависимости от знака \( x \), выражая асимметрию (например, хвост со стороны малых углов). В общем случае функция Пирсона VII обычно даёт лучшую аппроксимацию (меньший остаток), тогда как псевдо-Фойгта сходится более стабильно.

Обозначения

Символ	Значение
\( H_k \)	полная ширина на половине высоты (FWHM)
\( \pi \)	число «пи»
\( \eta \) (\( \eta_l, \eta_h \))	доля смешивания функций Лоренца/Гаусса (сторона малых углов / больших углов для split-форм)
\( \Gamma \)	гамма-функция
\( R \) (\( R_l, R_h \))	показатель степени Пирсона
\( A \)	параметр асимметрии
\( Z \)	нормировочная константа (\( \sqrt{\pi\ln 2} \))

Аппроксимирующая функция с фоном

На практике профильная функция \( f \) дополняется линейным фоном:

\[F(\theta,\Theta,B_1,B_2) = I\cdot f(\theta-\Theta) + B_1 + B_2(\theta-\Theta)\]

(\( I \): интегральная интенсивность, \( B_1, B_2 \): линейный фон, \( \Theta \): центр пика, \( \theta \): наблюдаемая позиция). В заданном диапазоне параметры варьируются методом Марквардта так, чтобы минимизировать \( R = \sum (Y - F)^2 \).

Частные производные каждой функции сложны; метод Марквардта использует эти аналитические градиенты. Ниже для справки приведены типичные выражения.

Частные производные симметричной функции псевдо-Фойгта

Обозначив \( u = \dfrac{\theta-\Theta}{H_k} \),

\[\frac{\partial F}{\partial I} = \frac{2}{H_k\pi}\left(\eta\,\frac{1}{1+4u^2}+(1-\eta)\sqrt{\pi\ln 2}\cdot e^{-4\ln 2\,u^2}\right)\]

\[\frac{\partial F}{\partial \eta} = \frac{2I}{H_k\pi}\left(\frac{1}{1+4u^2}-\sqrt{\pi\ln 2}\cdot e^{-4\ln 2\,u^2}\right)\]

\[\frac{\partial F}{\partial H_k} = -\frac{2I}{\pi}\left(\eta\,\frac{H_k^2-4(\theta-\Theta)^2}{\{H_k^2+4(\theta-\Theta)^2\}^2}+(1-\eta)\frac{\sqrt{\pi\ln 2}}{H_k^2}\,e^{-4\ln 2\,u^2}\right)\]

\[\frac{\partial F}{\partial \theta} = \frac{16(\theta-\Theta)I}{H_k^3\pi}\left(\eta\,\frac{1}{(1+4u^2)^2}+(1-\eta)\ln 2\sqrt{\pi\ln 2}\cdot e^{-4\ln 2\,u^2}\right)+B_2\]

\[\frac{\partial F}{\partial B_1} = 1,\qquad \frac{\partial F}{\partial B_2} = \theta-\Theta\]

Частные производные функции Пирсона VII

Простые производные по интенсивности и фону (\( \partial F/\partial I,\ \partial F/\partial B_1,\ \partial F/\partial B_2 \)) опущены. В оригинальном документе показатель степени Пирсона обозначается как \( R \), так и \( m \) (это одна и та же величина). Обозначив \( u = \dfrac{\theta-\Theta}{H_k} \),

\[\frac{\partial F}{\partial \Theta} = \frac{8mI(\theta-\Theta)}{H_k^2}(2^{1/R}-1)\left(1+4(2^{1/R}-1)u^2\right)^{-1-R}\]

\[\frac{\partial F}{\partial H_k} = \frac{8mI(\theta-\Theta)^2}{H_k^3}(2^{1/R}-1)\left(1+4(2^{1/R}-1)u^2\right)^{-1-R}\]

\[\frac{\partial F}{\partial m} = \left(1+4(2^{1/R}-1)u^2\right)^{-R}\left(\frac{4\cdot 2^{1/R}(\theta-\Theta)^2\ln 2}{m\left(H_k^2+4(2^{1/R}-1)(\theta-\Theta)^2\right)} - \ln\!\left(1+4(2^{1/R}-1)u^2\right)\right)\]

---

Вывод кубического сплайна

PDIndexer использует кривую кубического сплайна для построения фона. Истинную форму фона нельзя решить точно, но программа автоматически определяет области без пиков и соединяет обнаруженные точки сплайном, формируя кривую фона. Сплайн равномерно аппроксимирует данные, включая их производные, и приближение улучшается по мере уплотнения точек данных.

Пусть даны \( n \) точек \( (X_1,Y_1), (X_2,Y_2), \dots, (X_{n-1},Y_{n-1}), (X_n,Y_n) \); ищем кривую, которая кубична на каждом интервале и гладко соединяется так, что значение, наклон и кривизна совпадают в каждой точке (два крайних интервала \( \{-\infty, X_1\} \) и \( \{X_n, \infty\} \) считаются линейными).

Пусть функция на интервале \( \{X_{m-1}, X_m\} \) имеет вид

\[F_m = a_m X^3 + b_m X^2 + c_m X + d_m.\]

Внутренние точки (\( 2 \le m \le n-1 \)). Непрерывность значения, первой производной и второй производной даёт

\[F_m(X_m) = a_m X_m^3 + b_m X_m^2 + c_m X_m + d_m = Y_m\]

\[F_{m+1}(X_m) = a_{m+1} X_m^3 + b_{m+1} X_m^2 + c_{m+1} X_m + d_{m+1} = Y_m\]

\[F_m'(X_m) = 3a_m X_m^2 + 2b_m X_m + c_m = F_{m+1}'(X_m) = 3a_{m+1} X_m^2 + 2b_{m+1} X_m + c_{m+1}\]

\[F_m''(X_m) = 6a_m X_m + 2b_m = F_{m+1}''(X_m) = 6a_{m+1} X_m + 2b_{m+1}\]

— то есть \( 4n-8 \) условий.

Начало (\( m=1 \), левый крайний интервал линейный):

\[F_1(X_1) = c_1 X_1 + d_1 = Y_1\]

\[F_2(X_1) = a_2 X_1^3 + b_2 X_1^2 + c_2 X_1 + d_2 = Y_1\]

\[F_1'(X_1) = c_1 = F_2'(X_1) = 3a_2 X_1^2 + 2b_2 X_1 + c_2\]

\[F_1''(X_1) = F_2''(X_1) = 6a_2 X_1 + 2b_2 = 0\]

— 4 условия. Конец (\( m=n \)) аналогично даёт

\[F_n(X_n) = a_n X_n^3 + b_n X_n^2 + c_n X_n + d_n = Y_n\]

\[F_{n+1}(X_n) = c_{n+1} X_n + d_{n+1} = Y_n\]

\[F_n'(X_n) = 3a_n X_n^2 + 2b_n X_n + c_n = F_{n+1}'(X_n) = c_{n+1}\]

\[F_n''(X_n) = 6a_n X_n + 2b_n = F_{n+1}''(X_n) = 0\]

— ещё 4 условия.

В сумме \( 4n \) условий определяют \( 4n \) неизвестных, сводя задачу к системе линейных уравнений. Записав её в матричном виде и обратив матрицу, задачу легко решить.

---

Связанные страницы

• 6. Аппроксимация дифракционных пиков — практическое применение
• Уравнения состояния — теория EOS, включая уравнения Birch–Murnaghan и Mie–Grüneisen