Численное моделирование

Поскольку Schubart Master Formula — это двойной интеграл по энергии и углу, его вычисление для любого конкретного случая требует численных методов. В этой главе описывается общая методология этого: как дискретизировать интеграл, как задавать входные данные и как убедиться, что вычисленный результат сошёлся. На всём протяжении акцент делается на методе, а не на конкретных числах; любые упоминаемые величины суть иллюстративные заполнители, обозначенные как таковые, а не измеренные результаты.

Дискретизация интеграла

Исходной точкой является расширенное физическое представление,

$$P(t) = \iint \Phi(E,\theta,t)\,\sigma(E)\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;dE\,d\theta .$$

Самый прямой численный подход заменяет непрерывные интегралы суммами по сетке. Энергетическая ось разбивается на $N_E$ интервалов, а угловая — на $N_\theta$ интервалов, и подынтегральное выражение вычисляется в представительных точках каждой ячейки. Записывая подынтегральное выражение как $f(E,\theta,t)$, простая квадратура выглядит так:

$$P(t) \approx \sum_{i=1}^{N_E}\sum_{j=1}^{N_\theta} w_{ij}\,f(E_i,\theta_j,t)\,\Delta E_i\,\Delta\theta_j ,$$

где $\Delta E_i$ и $\Delta\theta_j$ — ширины ячеек, а $w_{ij}$ — квадратурные веса. Для правила средней точки веса равны единице, а точки выборки лежат в центрах ячеек. Более точные схемы — трапеций, Симпсона или квадратура Гаусса–Лежандра — заменяют вычисление в средней точке взвешенными суммами, снижающими ошибку дискретизации при заданном числе точек. Выбор схемы — это стандартный компромисс между точностью и стоимостью вычисления, и ничто в концепции не отдаёт предпочтения какой-либо одной квадратуре перед другой.

Стратегия дискретизации и порядок сходимости

Поскольку область — произведение энергетической и угловой осей, двумерное правило проще всего строится как тензорное произведение одномерных правил, одного для $E$ и одного для $\theta$. Эта факторизация важна и для стоимости, и для точности. Со стороны стоимости тензорная сетка использует $N_E \times N_\theta$ вычислений подынтегрального выражения, так что совместное измельчение обеих осей квадратично дорого; когда вся структура сосредоточена на одной оси — скажем, острый резонанс в $\sigma(E)$ вдоль энергетической оси, — часто гораздо дешевле измельчать только эту ось. Со стороны точности каждое одномерное правило несёт собственный порядок сходимости, показатель $p$ в асимптотическом масштабировании ошибки

$$\left| P_{\text{exact}} - P_{N} \right| \;\sim\; C\,N^{-p} ,$$

где $N$ — число точек вдоль оси, а $C$ зависит от гладкости подынтегрального выражения. Для правила трапеций $p = 2$, так что уменьшение шага вдвое снижает ошибку примерно вчетверо; для правила Симпсона $p = 4$ на гладком подынтегральном выражении; а для квадратуры Гаусса–Лежандра сходимость ещё быстрее — фактически экспоненциальная для аналитического подынтегрального выражения без особенностей в области. Эти порядки суть в точности классические результаты, сведённые в таблицы в стандартных справочниках, таких как Press и др., Numerical Recipes, и они переносятся на двойной интеграл через тензорную конструкцию. Практическое следствие в том, что достижимую точность задаёт гладкость подынтегрального выражения, а не номинальный порядок схемы: правило четвёртого порядка, применённое через почти разрыв в $\sigma(E)$, локально деградирует до поведения первого порядка, и глобальная ошибка тогда определяется этой одной плохо разрешённой областью.

Задание входных данных

Моделирование настолько осмысленно, насколько осмысленны его входные данные. Нужно предоставить четыре вида входных данных. Во-первых, модель потока $\Phi(E,\theta,t)$: описание окружающего спектра и его углового распределения, либо измеренное, либо взятое из предполагаемого распределения. Во-вторых, сечение $\sigma(E)$: энергозависимая связь, обычно задаваемая как табулированная или параметризованная функция; её аналитическое поведение обсуждается в анализе сигмы. В-третьих, параметры передачи импульса, переноса и эффективности $\Delta p$, $v_{\text{ph}}$ и $\eta_{\text{conv}}$, которые могут быть константами или сами зависеть от $E$ и $\theta$. Определения и единицы всех их зафиксированы в определениях переменных. Поскольку несколько из этих входных данных несут неопределённость, ответственное моделирование распространяет эту неопределённость, а не сообщает единственное детерминированное число, — к этому мы обращаемся в валидации методом Монте-Карло.

Проверка сходимости

Центральный вопрос контроля качества для детерминированной квадратуры — сошёлся ли результат по разрешению сетки. Стандартная процедура — вычислять $P(t)$ на последовательности всё более мелких сеток, удваивая $N_E$ и $N_\theta$ на каждом шаге, и наблюдать, как меняется оценка. Когда последовательные измельчения изменяют результат менее чем на выбранный допуск, сетка считается достаточной. Более количественная версия этого теста использует предполагаемый порядок сходимости напрямую: если ошибка масштабируется как $N^{-p}$, то вычисление интеграла на трёх сетках с шагами $h$, $h/2$ и $h/4$ и сравнение последовательных разностей выявляет как наблюдаемый порядок, так и, посредством экстраполяции Ричардсона, улучшенную оценку. Записывая $P_h$, $P_{h/2}$, $P_{h/4}$ для трёх результатов, отношение

$$\frac{P_{h} - P_{h/2}}{P_{h/2} - P_{h/4}} \;\approx\; 2^{p}$$

должно восстановить ожидаемый порядок; значение, далёкое от $2^p$, — предупреждение о том, что подынтегральное выражение не так гладко, как предполагалось, и что номинальной оценке ошибки доверять нельзя. Для гладкого подынтегрального выражения удачно выбранная квадратура сходится быстро; острые особенности в $\sigma(E)$ или в потоке, такие как узкий резонанс, требуют более мелкой выборки в области особенности и являются частой причиной медленной сходимости. Адаптивные схемы, которые дробят ячейки везде, где локальная оценка ошибки превышает порог, и тем самым размещают больше точек там, где подынтегральное выражение быстро меняется, суть естественное средство.

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

Оценка и бюджетирование ошибок

Обоснованный численный результат сопровождается явным бюджетом ошибок, разделяющим отдельные вклады. Ошибка дискретизации — это член $C\,N^{-p}$ выше, оцениваемый на практике разностью последовательных измельчений $|P_{h/2} - P_{h/4}|$, которая для сходящейся последовательности ограничивает оставшуюся ошибку с точностью до малого множителя. Ошибка усечения возникает от обрезания энергетической области на конечном верхнем пределе; она ограничена интегралом по отброшенному хвосту, и область следует расширять, пока эта граница не упадёт ниже ошибки дискретизации. Входная ошибка, распространённая от неопределённости в $\Phi$, $\sigma(E)$, $\Delta p$, $v_{\text{ph}}$ и $\eta_{\text{conv}}$, обычно является наибольшей из трёх и обрабатывается стохастическими методами валидации методом Монте-Карло, а не измельчением сетки. Распространённая и поучительная ошибка — измельчать сетку до машинной точности, игнорируя входную неопределённость в десятки процентов, — полировать число, истинная ошибка которого определяется множителем, до которого квадратуре не дотянуться. Дисциплина состоит в том, чтобы сообщать все три вклада и позволить наибольшему из них управлять тем, сколько усилий заслуживает каждое дальнейшее измельчение.

Источники численной ошибки

Помимо ошибки дискретизации, внимательный практик следит за несколькими повторяющимися ловушками. Слишком раннее усечение энергетического интеграла отбрасывает поток в хвостах; область следует расширять, пока подынтегральное выражение не станет действительно пренебрежимо малым, как отмечено в математической архитектуре. Использование угловой сетки, игнорирующей правильную меру, — например, опущение множителя $\sin\theta$ там, где интегрирование по телесному углу его требует, — вносит систематическое смещение, которое никаким измельчением не устранить, поскольку это ошибка в подынтегральном выражении, а не в его выборке. А интерполяция грубо табулированного $\sigma(E)$ через острую особенность может недооценить или переоценить перекрытие с потоком. Каждой из этих ошибок можно избежать стандартной численной гигиеной, и каждую следует документировать вместе с любым сообщаемым результатом.

От детерминированного к статистическому вычислению

Детерминированная квадратура отвечает на вопрос «при данных точных входных данных чему равно $P(t)$?». Сама по себе она не отвечает на вопрос «при данной неопределённости этих входных данных, насколько мы должны быть уверены в $P(t)$?». Второй вопрос требует стохастических методов, в которых входные данные выбираются из их распределений, а интеграл оценивается статистически. Эти методы и порождаемые ими доверительные интервалы суть предмет валидации методом Монте-Карло и анализа сигмы. Сошедшаяся детерминированная квадратура выполняет там и вторую роль: запущенная на единственном фиксированном наборе входных данных, она является эталоном, с которым перекрёстно сверяется оценщик Монте-Карло, так что два метода скорее дополняют друг друга, чем являются альтернативами.

Дополнительное чтение

• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press — стандартный справочник по квадратурным схемам, порядку сходимости и адаптивному интегрированию.
• Stoer, J., & Bulirsch, R. (2002). Introduction to Numerical Analysis (3rd ed.). Springer — строгое изложение анализа ошибок квадратуры и экстраполяции Ричардсона.
• Davis, P. J., & Rabinowitz, P. (1984). Methods of Numerical Integration (2nd ed.). Academic Press — исчерпывающе о гауссовой квадратуре и многомерном интегрировании.