Журнал вычислительной математики и математической физики, 2020, T. 60, № 7, стр. 1095-1110

1. ВВЕДЕНИЕ

Использование квадратурных формул на сфере, позволяющих получить необходимую точность расчета полей излучения при минимальном числе угловых направлений, является важным компонентом численного решения уравнения переноса методом дискретных ординат. Обычно используемые для решения уравнения переноса квадратуры на сфере: ${{S}_{n}}$ [2] и $E{{S}_{n}}$ [3] квадратуры Карлсона порядка $n$, $n = 2,4,6, \ldots $, содержащие $n(n + 2)$ узлов на сфере в $x,\;y,\;z$ геометрии и $n(n + 4)$ узлов в $r,\;\vartheta ,\;z$ геометрии, из которых $2n$ узлов являются вспомогательными, имеющими нулевые веса. $E{{S}_{n}}$ квадратура имеет одинаковые веса. Их размер ${{w}_{0}}$ зависит только от порядка квадратуры $n$: ${{w}_{0}} = {{4\pi } \mathord{\left/ {\vphantom {{4\pi } {(n(n + 2))}}} \right. \kern-0em} {(n(n + 2))}}$. ${{S}_{n}}$ и $E{{S}_{n}}$ квадратуры имеют относительно низкую алгебраическую точность: они точны для всех сферических гармоник $Y_{l}^{m}(\theta ,\varphi )$ с $l \leqslant n - 1$ [4] и $l \leqslant 3$ соответственно. Квадратура Лежандра–Чебышёва $L{{C}_{n}}$ типа произведения порядка $n$, $n = 2,4,6, \ldots $, на сфере [5]:

где $\left\{ {{{w}_{l}},{{\mu }_{l}}} \right\}$ – набор весов и узлов гауссовой квадратуры на интервале $\left[ { - 1,1} \right]$ порядка $n$, точна для всех сферических гармоник $Y_{l}^{m}(\theta ,\varphi )$ с $l \leqslant 2n - 1$, однако использует $2{{n}^{2}}$ узлов в $x,\;y,\;z$ геометрии и $2n\left( {n + 1} \right)$ узлов в $r,\;\vartheta ,\;z$ геометрии, что значительно больше, чем у $E{{S}_{n}}$ и ${{S}_{n}}$ квадратур этого же порядка.

Вместо использования квадратуры Чебышёва порядка $2n$ на каждом слое $l$ можно предложить вариант квадратуры Лежандра–Чебышёва [2], [6], [7], в которой порядок квадратуры Чебышёва зависит от номера слоя $l$, возрастая от минимального значения 2 у полюсов при $l = 1$ и $l = n$ до максимального значения $2n$ вблизи экватора при $l = {n \mathord{\left/ {\vphantom {n 2}} \right. \kern-0em} 2} - 1$ и $l = {n \mathord{\left/ {\vphantom {n 2}} \right. \kern-0em} 2} + 1$:

В [6] квадратура (1.2) обозначена как ${{P}_{n}} - {{T}_{n}}$ квадратура. Другое название этой квадратуры, которого мы будем придерживаться: Legendre–Chebyshev Triangular ($LC{{T}_{n}}$) квадратура [7]. Расположение узлов на сфере у этой квадратуры аналогично принятому у $E{{S}_{n}}$ и ${{S}_{n}}$ квадратур. Коэффициент эффективности квадратуры на сфере [8]: равен отношению числа правильно интегрируемых этой квадратурой сферических гармоник $Y_{l}^{m}(\theta ,\varphi )$ с $l \leqslant L$к числу параметров $3N$ квадратурной формулы с числом узлов $N$. Для $L{{C}_{n}},$ ${{S}_{n}}$ и $E{{S}_{n}}$ квадратуры в $x,\;y,\;z$ геометрии он равен соответственно ${2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-0em} 3}$, ${n \mathord{\left/ {\vphantom {n {3\left( {n + 2} \right)}}} \right. \kern-0em} {3\left( {n + 2} \right)}}$ и ${{16} \mathord{\left/ {\vphantom {{16} {3n\left( {n + 2} \right)}}} \right. \kern-0em} {3n\left( {n + 2} \right)}}$. Порядок ${{S}_{n}}$ квадратур ограничен, так как при $n > 20$ у квадратуры появляются отрицательные веса. $L{{C}_{n}}$, $LC{{T}_{n}}$ и $E{{S}_{n}}$ квадратуры не имеют подобных ограничений. Отметим также, что $L{{C}_{n}}$, $LC{{T}_{n}}$, ${{S}_{n}}$ и   $E{{S}_{n}}$ квадратуры обладают симметрией относительно отражений: $\theta \to - \theta $, $\varphi \to - \varphi $, ${\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2} - \varphi \to {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2} + \varphi $.

Целью данной работы является разработка алгоритма для построения квадратурных формул, обладающих симметрией правильной призмы (диэдра) и точных для всех сферических многочленов со степенью, не превышающей некоторого максимального значения $L$. Узлы у построенных квадратурных формул должны располагаться на так называемых “орбитах” (множестве точек с одним и тем же полярным углом). Данная работа является развитием работы А.Н. Казакова и В.И. Лебедева [1]. Построенная квадратура, в отличие от [1], не содержит узлов при $\varphi = 0,{\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2},\pi ,{{3\pi } \mathord{\left/ {\vphantom {{3\pi } 2}} \right. \kern-0em} 2}$, на полюсах $\theta = \pm {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2}$ и экваторе $\theta = 0$ сферы.

На возможность использования квадратуры [1] в $r,\;\vartheta ,\;z$ геометрии внимание одного из авторов этой работы (А.М. Волощенко) обратил В.И. Лебедев. Поиск программы QR, написанной А.Н. Казаковым для расчета квадратуры [1], к сожалению, оказался безрезультативным. Это вызвало необходимость разработки новой программы, оформленной в виде подпрограммы momscall, обращение к которой позволяет получить набор направляющих косинусов и весов квадратуры заданного порядка точности. Одновременно было выполнено уточнение вида требуемой квадратуры и алгоритма ее расчета. Подпрограмма momscall, написанная А.А. Руссковым, затем была интегрирована в 2D и 3D ${{S}_{N}}$ коды КАСКАД-С и КАТРИН из пакета CNCSN [9], позволяющие проводить расчеты радиационной защиты реакторных установок в 2D $x,\;z$, $r,\;z$ и $r,\;\vartheta $ и 3D $x,\;y,\;z$ и $r,\;\vartheta ,\;z$ геометриях.

В статье представлены результаты сопоставления точности квадратуры Казакова–Лебедева с вышеописанными ограничениями, обладающей симметрией группы диэдра, с обычно используемыми для решения уравнения переноса в 3D $x,\;y,\;z$ и $r,\;\vartheta ,\;z$ геометриях $E{{S}_{n}}$, ${{S}_{n}}$, $L{{C}_{n}}$ и $LC{{T}_{n}}$ квадратурами, в том числе, при решении практической задачи радиационной защиты.

2. РАСШИРЕННАЯ ГРУППА ДИЭДРА И ОБЩИЙ ВИД КВАДРАТУРНОЙ ФОРМУЛЫ С СООТВЕТСТВУЮЩЕЙ СИММЕТРИЕЙ

Пусть

есть единичная сфера в пространстве ${{R}^{3}}$. Рассмотрим поверхностный интеграл по сфере вида где $w(s) \geqslant 0$ – интегрируемая на сфере $S$ весовая функция, такая что $I[1] = 1$.

Введем полярные координаты $r,\varphi ,\theta $: $x = r\sin \theta \cos \varphi $, $y = r\sin \theta \sin \varphi $, $z = r\cos \theta $. Тогда $ds = \sin \vartheta d\varphi d\vartheta $. В наиболее общем виде квадратурная формула для сферы выглядит следующим образом:

Для хороших квадратурных формул типа Гаусса значение коэффициента эффективности квадратуры $\eta $ (1.3) достаточно близко к единице $\eta \approx 1$. Рассмотрим задачу о построении квадратурной формулы типа Гаусса для поверхностного интеграла (2.1). Алгебраическая степень точности для квадратурной формулы в случае поверхностного интеграла определяется аналогично одномерному случаю как максимальная степень многочлена из некоторого класса, точно интегрируемого по поверхности сферы.

Следуя [1], введем понятие расширенной группы диэдра $\overline {{{D}_{m}}} $. Пусть $m \geqslant 1$ – натуральное число; $U$, $V$ и $T$ – следующие преобразования пространства ${{R}^{3}}$:

т.е. поворот вокруг оси $Oz$ на угол ${{2\pi } \mathord{\left/ {\vphantom {{2\pi } m}} \right. \kern-0em} m}$ и зеркальная симметрия относительно плоскостей $xOz$ и $xOy$ соответственно. Эти преобразования порождают группу, которую, следуя терминологии Ф. Клейна, назовем расширенной группой диэдра $\overline {{{D}_{m}}} = \left\langle {U,V,T} \right\rangle $. Порядок группы $\overline {\left| {{{D}_{m}}} \right|} = 4m$. Циклическая группа ${{C}_{m}} = \left\langle U \right\rangle $ порядка ${{C}_{m}} = m$ и группа диэдра ${{D}_{m}} = \left\langle {U,V} \right\rangle $ порядка $\left| {{{D}_{m}}} \right| = 2m$ являются нормальными подгруппами группы $\overline {{{D}_{m}}} $. В случае, когда $m$ – четное, $\overline {{{D}_{m}}} $ имеет еще одну нормальную подгруппу порядка $2$, неединичный элемент которой есть преобразование инверсии относительно начала координат: $(x,y,z) \in {{R}^{3}} \to ( - x, - y, - z)$.

Возможна следующая геометрическая интерпретация преобразований группы $\overline {{{D}_{m}}} $. Пусть $P{{R}_{m}}$ – правильная призма с $m$ боковыми гранями такая, что ее ось совпадает с осью $Oz$, а ребра, параллельные этой оси, делятся плоскостью $xOy$ пополам, причем одно из ребер пересекает положительную полуось $Ox$. Ясно, что при преобразованиях группы $\overline {{{D}_{m}}} $ призма $P{{R}_{m}}$ переходит в себя.

Поскольку мы будем строить квадратурные формулы, симметричные относительно группы $\overline {{{D}_{m}}} $, потребуем, чтобы весовая функция $w(s)$ была инвариантна при этих преобразованиях, и, кроме того, имела вид

где $q(\varphi )$, $\varphi \in [0,2\pi ]$ и $p(z)$, $z \in [ - 1,1]$ – неотрицательные функции, удовлетворяющие условиям инвариантности: Для случая нечетного $m$ на функцию $q(\varphi )$ наложим дополнительное условие Квадратурные формулы для интеграла (2.1), инвариантные относительно группы $\overline {{{D}_{m}}} $, следуя [1], будем искать в виде: Узлы квадратурной формулы (2.3) образуют инвариантное относительно группы $\overline {{{D}_{m}}} $ множество и являются объединением $\overline {{{D}_{m}}} $ – орбит некоторых точек, при этом узлы, лежащие на одной и той же орбите $\overline {{{D}_{m}}} $ имеют одинаковые веса. Количество точек, лежащих на одной и той же орбите $\overline {{{D}_{m}}} $, может быть различным. Для каждой орбиты мы выберем узел, который назовем представителем орбиты; для удобства выберем узел, координаты которого задаются наиболее просто. Зная представителя орбиты, мы можем, применяя преобразования группы $\overline {{{D}_{m}}} $, получить все узлы этой орбиты. В табл. 1 приводятся представители каждой орбиты, при этом ${{r}_{k}} = \sqrt {1 - z_{k}^{2}} $.

В табл. 1 узлы ${{a}^{i}}$ – полюса сферы (точки $\left( {0,0, \pm 1} \right)$); узлы ${{b}^{i}}$, ${{c}^{i}}$, $d_{{0j}}^{i}$ лежат на сфере в плоскости $xOy$, т.е. на экваторе; узлы $b_{k}^{i}$, $c_{k}^{i}$, $d_{{kj}}^{i}$, $k > 0$ лежат на сфере в плоскостях $z = {{z}_{k}}$, $z = - {{z}_{k}}$. Как уже отмечалось выше, при использовании для решения уравнения переноса в $r,\;\vartheta ,\;z$ геометрии квадратура не должна содержать узлов при $\varphi = 0,{\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2},\pi ,{{3\pi } \mathord{\left/ {\vphantom {{3\pi } 2}} \right. \kern-0em} 2}$, на полюсах $\theta = \pm {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2}$ и экваторе $\theta = 0$ сферы. Этому условию удовлетворяют только орбиты $d_{{kj}}^{i},i = 1, \ldots ,4m$.

В дальнейшем будем предполагать следующее упорядочение плоскостей по индексу $k$: $0 < {{z}_{1}} < \cdots < {{z}_{{\overline N }}} < 1$.

Метод получения квадратурной формулы (2.3) заданной степени точности $n$ основан на использовании основной в теории квадратурных формул теоремы С.Л. Соболева [11]. Сформулируем эту теорему применительно к рассмотренному случаю.

Теорема Соболева. Для того, чтобы квадратурная формула (2.3) была точна для следов на $S$ всех многочленов до степени $n$ включительно, необходимо и достаточно, чтобы она была точна для инвариантных относительно группы $\overline {{{D}_{m}}} $ следов многочленов до степени $n$ включительно.

Алгоритм построения всех инвариантных относительно группы $\overline {{{D}_{m}}} $ многочленов основан на следующей теореме, при доказательстве которой были использованы результаты работ [12], [13].

Теорема 1. Алгебра всех инвариантных относительно группы $\overline {{{D}_{m}}} $ многочленов в ${{R}^{3}}$ порождается многочленами

где $r = \sqrt {{{x}^{2}} + {{y}^{2}}} $, $x = r\cos \varphi $, $y = r\sin \varphi $, $C_{n}^{k} = \frac{{n!}}{{k!(n - k)!}}$.

На сфере $S$ ${{\sigma }_{1}} \equiv 1$, поэтому имеет место

Следствие. Любой многочлен, инвариантный относительно группы $\overline {{{D}_{m}}} $, представим на сфере $S$ в виде следа многочлена от ${{\sigma }_{2}} = {{z}^{2}}$ и ${{\sigma }_{3}} = {{\left( {1 - z} \right)}^{{m/2}}}\cos (m\varphi )$.

3. СИСТЕМА УРАВНЕНИЙ ДЛЯ НАХОЖДЕНИЯ ПАРАМЕТРОВ КВАДРАТУРНЫХ ФОРМУЛ, ОБЛАДАЮЩИХ СИММЕТРИЕЙ ГРУППЫ ДИЭДРА

Размерность пространства многочленов на сфере $S$ до степени $n$ включительно ${{W}_{{n.m}}}$, инвариантных относительно группы $\overline {{{D}_{m}}} $, равна

где ${{\chi }_{k}}$ – число решений в целых неотрицательных числах уравнения $2{{\alpha }_{1}} + m{{\alpha }_{2}} = k$. Базис ${{W}_{{n.m}}}$ образует следы на сфере $S$ многочленов ${{e}_{{ji}}} = \sigma _{2}^{j}\sigma _{3}^{i}$, где $i$, $j$ принимают всевозможные значения, такие что $0 \leqslant 2j + mi \leqslant n$. Для удобства дальнейших выкладок разобьем многочлены ${{e}_{{ji}}}$ по группам с фиксированным $i$, группа ${{E}_{0}}$: группы ${{E}_{i}}$: Введем также следующие обозначения: Здесь через $\overline N $ обозначено количество параллелей в полусфере, на которых расположены узлы, за исключением экватора и полюсов. Условия того, что квадратурная формула (2.3) является точной для всех функций пространства ${{W}_{{n.m}}}$, т.е. равенство $I[f] = {{I}_{n}}[f]$ для следов на сфере $S$ всех многочленов групп ${{E}_{i}}$, $0 \leqslant i \leqslant {{i}_{{\max }}}$, эквивалентны системе нелинейных алгебраических уравнений, которая за счет разбиения (3.1), (3.2) и введения обозначений для переменных ${{t}_{{ik}}}$, $k = 0,1, \ldots $, имеет блочно-треугольный вид, подсистема ${{E}_{0}}$: Подсистемы ${{E}_{i}}$, $1 \leqslant i \leqslant {{i}_{{\max }}}$: Обозначим через ${{l}_{i}}$ количество уравнений в подсистеме ${{E}_{i}}$, а через $l$ общее количество уравнений в системе (3.5)–(3.8): Сделаем вспомогательные замечания для решения системы (3.5)–(3.8).

Замечание 1. Поскольку весовая функция $w(s)$ в интеграле (2.1) удовлетворяет условию (2.2), то правые части системы (3.5)–(3.8) равны:

где

Если для интеграла (3.11) существует квадратурная формула определенного порядка точности, то при условии (3.9) сомножитель $J\left[ {{{z}^{{2j}}}{{{(1 - {{z}^{2}})}}^{{\frac{{mi}}{2}}}}} \right]$ в правых частях (3.10) выражается непосредственно по квадратурной формуле, поскольку при четном $m$ функция ${{z}^{{2j}}}{{(1 - {{z}^{2}})}^{{\frac{{mi}}{2}}}}$ является многочленом. При нечетном $m$ функция ${{z}^{{2j}}}{{(1 - {{z}^{2}})}^{{\frac{{mi}}{2}}}}$ не является многочленом, но в этом случае значение интеграла не существенно в силу обращения сомножителя ${{C}_{i}}$ в ноль.

Замечание 2. Для краткости будем называть узлы квадратурной формулы (2.3), расположенные в плоскостях $z = \pm {{z}_{k}}$ (т.е. узлы $b_{k}^{i}$, $c_{k}^{i}$, $d_{{kj}}^{i}$, $k \geqslant 1$), узлами на уровне ${{z}_{k}}$. Переменные ${{t}_{{i0}}}$, ${{t}_{{ik}}}$, входящие в левую часть системы (3.5)–(3.8), были введены по формулам (3.3), (3.4), поэтому для того чтобы получить из них параметры квадратурной формулы, необходимо, учитывая, что ${{\sigma }_{{0j}}} = {{\gamma }_{{0j}}}$ и ${{\sigma }_{{kj}}} = r_{k}^{m}{{\gamma }_{{kj}}}$, $k \geqslant 1$, решить следующие системы уравнений:

a) для параметров узлов на уровне ${{z}_{k}}$, $k = 1, \ldots ,\overline N $:

б) для параметров узлов на экваторе:

Координаты представителей орбит находятся по формулам: Значения ${{i}_{0}}$, ${{i}_{k}}$, $k = 1, \ldots ,\overline N $, определяются из условия разрешимости системы (3.5)–(3.8) и зависят от конкретного вида квадратурной формулы (2.3). В частности, в квадратурную формулу (2.3) могут как входить, так и не входить фиксированные узлы ${{b}^{i}}$, ${{c}^{i}}$, $b_{k}^{i}$, $c_{k}^{i}$. Значение ${{N}_{k}}$, $k = 1, \ldots ,\overline N $, значение ${{N}_{0}}$ выбираются таким образом, чтобы общее число неизвестных в системе (3.12) (в системе (3.13)) было равно числу уравнений.

Замечание 3. Системы (3.12) и (3.13) после несложных преобразований сводятся к моментным системам уравнений от одного переменного. Аналогичные системы возникают при нахождении параметров квадратурных формул Гаусса и Гаусса–Маркова. Рассмотрим метод решения моментной системы. Пусть задана моментная система

где нужно найти ${{B}_{i}}$, ${{x}_{i}}$ по известным величинам ${{c}_{k}}$. Моментная система (3.14) может быть решена следующим образом. Будем считать, что задан отрезок $[a,b]$ и задана весовая функция $p(x)$ такая, что Построим последовательность взаимно ортогональных многочленов $\left\{ {{{F}_{n}}(x)} \right\}_{0}^{p}$ на отрезке $[a,b]$ с весовой функцией $p(x)$: При построении воспользуемся рекуррентными соотношениями где ${{F}_{{ - 1}}} = 0$, ${{F}_{0}} = 1$, ${{\beta }_{1}} = 0$ и в силу ортогональности Поэтому, введя обозначения получаем следующие рекуррентные формулы для нахождения коэффициентов многочленов ${{F}_{n}}(x)$: Далее находим корни многочлена ${{F}_{p}}(x)$, после чего находим неизвестные коэффициенты ${{B}_{i}}$: Вывод соотношения (3.15) для весов моментной системы (3.14) можно найти в [14].

4. ПРЯМАЯ РЕАЛИЗАЦИЯ КВАДРАТУРНОЙ ФОРМУЛЫ С ПОВЫШЕННЫМ КОЛИЧЕСТВОМ ТОЧЕК

Пусть ${{Q}_{n}}[{v}]$ – какая-либо симметричная квадратурная формула для интеграла (3.11) алгебраической степени точности $n$:

Так как она точна для всех многочленов от $z$ до степени $n$ включительно, то в силу (3.10) имеем Поэтому, положив получим одно из решений системы (3.5)–(3.8). После этого, решая $\overline N $ систем (3.12), в каждой из которых ${{i}_{k}} = {{i}_{{\max }}}$, $k = 1, \ldots ,\overline N $, и при наличии узлов на экваторе, систему (3.13)), в которой ${{i}_{0}} = {{i}_{{\max }}}$, находим все параметры квадратурной формулы.

При решении системы (3.5)–(3.8) таким методом число параметров квадратурной формулы (2.3) для узлов, лежащих в каждой из плоскостей $z = \pm {{z}_{k}}$, $k = 1, \ldots ,\overline N $, и $z = 0$ при наличии узлов на экваторе, одинаково и равно $m({{i}_{{\max }}} + 1)$, при этом для полного числа параметров справедлива формула

Коэффициент эффективности этих квадратурных формул сравнительно невелик, для них он Если зафиксировать $n$ и $m$, то минимальным числом узлов будут обладать квадратурные формулы, в которых координатами узлов по оси $z$ являются узлы квадратурной формулы (4.1) Гаусса и Гаусса–Маркова с двумя фиксированными узлами $z = \pm 1$.

5. РЕАЛИЗАЦИЯ КВАДРАТУРНОЙ ФОРМУЛЫ С ПОВЫШЕННОЙ ЭФФЕКТИВНОСТЬЮ

Приведем теперь способ решения системы (3.5)–(3.8), который дает квадратурные формулы типа Гаусса с асимптотикой

Пусть (4.1) – квадратурная формула Гаусса или Гаусса–Маркова с двумя фиксированными узлами $z = \pm 1$. В этом случае решение подсистемы ${{E}_{0}}$ (3.5), (3.6) находим по формулам (4.2). После подстановки значений ${{z}_{k}}$ в остальные уравнения системы (3.7), (3.8) получим ${{i}_{{\max }}}$ линейных подсистем с матицами Вандермонда. Заметим, что число уравнений в подсистемах ${{E}_{i}}$ убывает с ростом номера $i$: ${{l}_{i}} < {{l}_{{i{\text{'}}}}}$, если $i > i{\text{'}}$. Обозначим через $\bar {i}$ номер подсистемы ${{E}_{{\bar {i}}}}$, для которой выполнено следующее: Поскольку мы считаем все ${{t}_{{ik}}}$ неизвестными, $\bar {i}$ – номер последней подсистемы ${{E}_{{\bar {i}}}}$, в которой число неизвестных меньше или равно числу уравнений.

Подсистемы ${{E}_{i}}$, где $1 \leqslant i \leqslant \bar {i}$, являются совместными, их решение находим по формулам (4.3). Рассмотрим теперь подсистему ${{E}_{{\bar {i} + 1}}}$, в которой число неизвестных больше числа уравнений. Выберем какие-либо уровни ${{z}_{k}}$ общим числом $\overline N + {\text{sign}}\left| {{{p}_{0}}} \right| - {{l}_{{\bar {i} + 1}}}$ и определим параметры квадратурной формулы на этих уровнях исходя из имеющихся к этому моменту значений ${{t}_{{ik}}}$, где $i = 0, \ldots ,\bar {i}$. Уровни ${{z}_{k}}$ целесообразно выбрать ближайшими к полюсу, т.е. $k = {{l}_{{\bar {i} + 1}}} - {\text{sign}}\left| {{{p}_{0}}} \right| + 1, \ldots ,\overline N $, чтобы около полюса располагались уровни с наименьшим числом узлов.

Итак, для каждого $k = {{l}_{{\bar {i} + 1}}} - {\text{sign}}\left| {{{p}_{0}}} \right| + 1, \ldots ,\overline N $ решаем систему (3.12), в которой ${{i}_{k}} = \bar {i}$ и находим параметры квадратурной формулы на соответствующих уровнях. После этого найденные параметры подставляем в (3.4) и находим ${{t}_{{ik}}}$ при $i = \bar {i} + 1, \ldots ,{{i}_{{\max }}}$, $k = {{l}_{{\bar {i} + 1}}} - {\text{sign}}\left| {{{p}_{0}}} \right| + 1, \ldots ,\overline N $, которые подставляем в подсистемы ${{E}_{{\bar {i} + 1}}}, \ldots ,{{E}_{{{{i}_{{\max }}}}}}$. После этого в подсистеме ${{E}_{{\bar {i} + 1}}}$ число неизвестных становится равным числу уравнений. Решаем ее и находим неизвестные ${{t}_{{\bar {i} + 1,k}}}$, где $k = 1 - {\text{sign}}\left| {{{p}_{0}}} \right|, \ldots ,{{l}_{{\bar {i} + 1}}} - {\text{sign}}\left| {{{p}_{0}}} \right|$.

Поскольку ${{l}_{{\bar {i} + 2}}} < {{l}_{{\bar {i} + 1}}}$, то мы, действуя аналогично, находим параметры квадратурной формулы на уровнях ${{z}_{k}}$, где $k = {{l}_{{\bar {i} + 2}}} - {\text{sign}}\left| {{{p}_{0}}} \right| + 1, \ldots ,{{l}_{{\bar {i} + 1}}} - {\text{sign}}\left| {{{p}_{0}}} \right|$, в системах (3.12) в этом случае ${{i}_{k}} = \bar {i} + 1$, подставляем  в  (3.4) и  найденные ${{t}_{{ik}}}$, $i = \bar {i} + 2, \ldots ,{{i}_{{\max }}}$, $k = {{l}_{{\bar {i} + 2}}} - {\text{sign}}\left| {{{p}_{0}}} \right| + 1, \ldots ,{{l}_{{\bar {i} + 1}}} - {\text{sign}}\left| {{{p}_{0}}} \right|,$ подставляем в подсистемы ${{E}_{{\bar {i} + 2}}}, \ldots ,{{E}_{{{{i}_{{\max }}}}}}$. После этого в подсистеме ${{E}_{{\bar {i} + 2}}}$ число неизвестных становится равным числу уравнений. Решаем ее и находим неизвестные ${{t}_{{\bar {i} + 2,k}}}$, где $k = 1 - {\text{sign}}\left| {{{p}_{0}}} \right|, \ldots ,{{l}_{{\bar {i} + 2}}} - {\text{sign}}\left| {{{p}_{0}}} \right|$.

Продолжая этот процесс, мы решим все подсистемы ${{E}_{i}}$, $i = 1, \ldots ,{{i}_{{\max }}}$. После решения системы ${{E}_{{{{i}_{{\max }}}}}}$ определяем параметры квадратурной формулы на уровнях ${{z}_{k}}$, где $k = 1, \ldots ,{{l}_{{{{i}_{{\max }}}}}} - {\text{sign}}\left| {{{p}_{0}}} \right|$, в системе (3.12) в этом случае ${{i}_{k}} = {{i}_{{\max }}}$. После этого, если ${{p}_{0}} \ne 0$, решаем систему (3.13) с ${{i}_{0}} = {{i}_{{\max }}}$ и находим параметры квадратурной формулы на экваторе.

Данный метод дает квадратурные формулы типа Гаусса с числом параметров

и асимптотикой (5.1).

В табл. 2 приведены формулы для $\bar {i}$ и формулы, связывающие $n$ и $\overline N $ для четырех возможных случаев.

В разработанной программе задаются следующие параметры: порядок симметрии $m$ и число уровней $q \equiv \overline N $. Параметры $A{\text{'}}$ и ${{p}_{0}}$ полагаются равными нулю. Далее программа вычисляет порядок точности квадратурной формулы $n$ и осуществляет ее построение. При нахождении параметров на уровнях учитываются только орбиты $d_{{kj}}^{i},i = 1, \ldots ,4m$.

Далее порядком $n$ построенной квадратуры Казакова–Лебедева с индексом аксиальной симметрии $m$, которую мы обозначим через $K{{L}_{{n,m}}}$, по аналогии с $L{{C}_{n}}$ ${{S}_{n}}$, $E{{S}_{n}}$ и $LC{{T}_{n}}$квадратурами, мы будем называть число $n = 2q$, где $q$ – число уровней (параллелей) квадратуры для интервала углов $0 < \theta < {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. \kern-0em} 2}$. Точность квадратуры мы будем характеризовать параметром $L$, определяющим максимальный порядок правильно интегрируемых сферических гармоник. Для квадратуры Казакова–Лебедева $K{{L}_{{n,m}}}$ порядок квадратуры $L = 2n - 1$.

В табл. 3 для случая $x,\;y,\;z$ геометрии приведены параметры квадратур Казакова–Лебедева $K{{L}_{{n,m}}}$ порядка $n = 2,4,6, \ldots ,32$, соответствующие минимально возможным значениям индекса симметрии квадратуры $m = 2,3, \ldots $ . Для каждой квадратуры приведены полное число узлов квадратуры $N$, максимальный порядок $L$ правильно интегрируемых сферических гармоник и коэффициент эффективности $\eta $. Начиная с квадратуры порядка $n = 26$ возникает необходимость в постепенном увеличении индекса $m$ для достижения устойчивости алгоритма построения квадратуры. При этом коэффициент эффективности квадратуры начинает постепенно уменьшаться. Прочерками в таблице обозначена невозможность построении квадратуры с данными параметрами $n$ и $m$ (выход косинуса угла ${{\gamma }_{{kj}}} = \cos (m{{\varphi }_{{kj}}})$ за пределы области значения косинуса $\left[ { - 1,1} \right]$). Максимальный порядок квадратуры Казакова–Лебедева рассмотренного вида, который удается построить: $n = 32$.

На фиг. 1–4 для 1-го октанта изображены узлы квадратуры Карлсона $E{{S}_{{16}}}$ (содержащей 288 узлов и имеющую алгебраическую точность $L = 3$), квадратуры Казакова–Лебедева $K{{L}_{{12,2}}}$ (240 узлов, $L = 23$), квадратуры Казакова–Лебедева $K{{L}_{{24,2}}}$ (912 узлов, $L = 47$) и триангулярной квадратуры Лежандра–Чебышёва $LC{{T}_{{16}}}$ (288 узлов) соответственно.

6. СРАВНЕНИЕ ТОЧНОСТИ $K{{L}_{{n,m}}}$, $E{{S}_{n}}$, ${{S}_{n}}$, $L{{C}_{n}}$ И $LC{{T}_{n}}$ КВАДРАТУР ПРИ РАСЧЕТЕ ЧЕТНЫХ УГЛОВЫХ МОМЕНТОВ

В табл. 4 представлена максимальная относительная ошибка $\varepsilon _{{err}}^{k}$ (5.3) при расчете четных угловых моментов (5.4) для пяти типов квадратурных формул порядка $n = 16$.

Здесь где ${{\xi }_{{l,m}}}$, ${{\eta }_{{l,m}}}$ и ${{\mu }_{l}}$ – направляющие косинусы единичного вектора ${{\vec {\Omega }}_{{l,m}}}$ угловых направлений: ${{w}_{{l,m}}}$ – веса квадратурной формулы. Точные значения четных угловых моментов (5.4) согласно [6], [7] равны: Из табл. 4 следует, что квадратуры $L{{C}_{n}}$, $LC{{T}_{n}}$ и $K{{L}_{{n,2}}}$ при расчете четных угловых моментов (5.4) и $n = 16$ имеют примерно одинаковую точность.

7. СРАВНЕНИЕ ТОЧНОСТИ $K{{L}_{{n,m}}}$, $E{{S}_{n}}$, ${{S}_{n}}$, $L{{C}_{n}}$ И $LC{{T}_{n}}$ КВАДРАТУР В ЗАДАЧЕ РАСЧЕТА РАДИАЦИОННОЙ ЗАЩИТЫ ВВЭР-1000

На фиг. 5 и 6 представлены результаты расчета в ${{P}_{3}}$ приближении с точностью сходимости внутренних итераций ${{10}^{{ - 3}}}$ с использованием проблемно-ориентированной библиотеки сечений BGL1000 [15], содержащей 47 групп нейтронов, источника нейтронов деления в активной зоне реактора, нормированного на 3.7447E+19 нейтрон/сек на сектор поворотной симметрии 60°, азимутальной зависимости плотности потока нейтронов с энергией E > 3.0 МэВ на внешней поверхности при $r = {\text{226}}{\text{.75}}$ см корпуса реактора ВВЭР-1000 вблизи максимума аксиального распределения нейтронов $z = {\text{223}}{\text{.8}}$ см (122.2 см от низа АЗ). В табл. 5 приведены значения плотности потока нейтронов с энергией E > 3.0 МэВ в максимуме азимутального распределения в зависимости от используемой квадратуры. В табл. 6 для случая 3D $r,\;\vartheta ,\;z$ геометрии приведены число узлов квадратуры (включая вспомогательные узлы с нулевыми весами) в зависимости от ее типа, порядка и алгебраической точности.

Из табл. 5 и фиг. 5–6 можно сделать вывод, что использование $LC{{T}_{{16}}}$, $L{{C}_{{16}}}$, $K{{L}_{{12,2}}}$, $K{{L}_{{14,2}}}$ и $K{{L}_{{16,2}}}$ квадратур обеспечивает погрешность расчета рассматриваемого функционала в точке максимума азимутального распределения на уровне ≤1%, а ${{S}_{{16}}}$ и $LC{{T}_{{12}}}$ квадратур на уровне ≤1.2%.

8. ЗАКЛЮЧЕНИЕ

В данной работе рассмотрен алгоритм и создана программа для построения $K{{L}_{{n,m}}}$ квадратуры Казакова–Лебедева для сферы, инвариантной относительно группы диэдра, с ограничениями, позволяющих использовать данную квадратуру для решения уравнения переноса в 3D $x,\;y,\;z$ и $r,\;\vartheta ,\;z$ геометриях. Программа оформлена в виде подпрограммы, которая интегрирована в 2D и 3D ${{S}_{N}}$ коды КАСКАД-С и КАТРИН, позволяющие проводить расчеты радиационной защиты РУ в 2D $x,\;z$, $r,\;z$ и $r,\;\vartheta $, и 3D $x,\;y,\;z$ и $r,\;\vartheta ,\;z$ геометриях, а также доступна для независимого использования по ссылке [16]. Выполнено тестирование разработанной программы, которое показало заявленную алгебраическую точность $K{{L}_{{n,m}}}$ квадратуры. Выполнено сопоставление точности и эффективности $K{{L}_{{n,m}}}$ квадратуры с обычно используемыми для решения уравнения переноса в 3D $x,\;y,\;z$ и $r,\;\vartheta ,\;z$ геометриях $E{{S}_{n}}$, ${{S}_{n}}$, $L{{C}_{n}}$ и $LC{{T}_{n}}$ квадратурами, в том числе, при решении практической задачи радиационной защиты.