Журнал вычислительной математики и математической физики, 2019, T. 59, № 8, стр. 1358-1380

1.1. Модель Томаса–Ферми

Согласно модели Томаса–Ферми многоэлектронного атома или иона при нулевой абсолютной температуре (см. [1]–[3, § 70]) безразмерная пространственная плотность заряда подчиняется нелинейному уравнению

где $\Delta $ – оператор Лапласа в ${{\mathbb{R}}^{3}}$. В сферически симметричном случае с учетом естественных требований к поведению $\Phi $ на бесконечности и вблизи точечного ядра (т.е. при $\left| {\mathbf{x}} \right| \to 0$) модель (1.1) приводит к следующей краевой задаче для сингулярного обыкновенного дифференциального уравнения второго порядка: где $r = {\text{const}}\left| {\mathbf{x}} \right|$ – безразмерное расстояние до ядра, имеющего заряд $ + Z$, $R \in (0, + \infty )$ – радиус иона ($R = + \infty $ для атома), коэффициент экранирования $\Psi = \Psi (r)$ определяет внутриатомный электростатический потенциал по формуле $\Phi = \Psi {\text{/}}\left| {\mathbf{x}} \right|$. Полный заряд иона $Z{\kern 1pt} {\text{'}}$, $0 < Z{\kern 1pt} {\text{'}} < Z$, может быть найден в данной модели по формуле Графики решений задачи (1.2), (1.3) при $Z = 1$ и различных значениях $R$ приведены на фиг. 1.

Отметим, что базовые принципы модели Томаса–Ферми находят применение при моделировании неоднородного электронного газа в рамках так называемого метода функционала плотности (см. [4]).

1.2. Уравнение Эмдена–Фаулера и некоторые свойства его решений

Рассмотрим следующую двухточечную краевую задачу для уравнения Эмдена–Фаулера (см. [5, гл. VIII]):

где $R$, $Z$ – некоторые заданные числа. Отметим, что постановка (1.2), (1.3) является частным случаем задачи (1.5)–(1.7) при значениях параметров $\nu = - 1{\text{/}}2$, $\sigma = 1{\text{/}}2$.

Известно (см. [5], [6]), что при $0 < R \leqslant + \infty $, $Z \in (0, + \infty )$ задача (1.5)–(1.7) имеет единственное классическое решение, которое является монотонно убывающей положительной аналитической функцией на интервале $r \in (0,R)$. Известно также (см. [5, гл. VII, п. 7]), что каждое положительное решение $y(r)$ уравнения (1.5) при $\nu > - 2$, определенное на луче , имеет асимптотический вид

где кроме того, уравнение (1.5) инвариантно по отношению к следующей замене переменных:

В частном случае задачи Томаса–Ферми (1.2), (1.3) параметры (1.9) принимают значения $\beta = 3$, ${{f}_{0}} = 12$, и формула масштабирования (1.10) позволяет (см. [2], [7]) выразить решение $\Psi (r)$ этой задачи через решение ${{\Psi }_{1}}(r)$ той же задачи при $Z = 1$ в виде

В этой связи функцию ${{\Psi }_{1}}(r)$ называют стандартным решением задачи Томаса–Ферми для нейтрального атома (см. [8, § 172]). Отметим, что для решения общей задачи (1.5)–(1.7) справедлива аналогичная формула

При отрицательных значениях показателя $\nu \in ( - 1,0)$ левый конец $r = 0$ рассматриваемого в задаче (1.5)–(1.7) промежутка $r \in [0,R]$ является сингулярной точкой уравнения, поскольку коэффициент при нелинейном члене обращается в бесконечность. В отношении задачи (1.2), (1.3) известно (см. [9]), что при любых $R \in (0, + \infty ]$ ее решение в окрестности точки $r = 0$ представляется сходящимся рядом Пюизё по степеням $\sqrt r $:

причем решениям при фиксированном $Z$ и различных $R$ соответствуют различные значения наклона графика в нуле $\Psi {\kern 1pt} {\text{'}}(0) = {{b}_{2}} \in ( - \infty ,B]$, где $B = B(Z)$ – некоторое отрицательное число. Коэффициенты ${{b}_{k}}$, $k \geqslant 3$ определяются (см. [10]) по $Z$ и ${{b}_{2}}$ однозначно из уравнения (1.2). Решению задачи (1.2), (1.3) на полубесконечном промежутке ($R = + \infty $) соответствует в разложении (1.12) критическое значение ${{b}_{2}} = B(Z)$, которое при $Z = 1$ равно Отметим, что величина ${{b}_{2}}$ имеет в модели Томаса–Ферми следующий физический смысл: с точностью до множителя это энергия кулоновского взаимодействия электронного облака с ядром (см. [8, § 174]).

Известно, что решение задачи (1.2), (1.3) при $R = + \infty $ в окрестности бесконечности может быть приближенно найдено по следующей формуле (так называемое представление Коулсона–Марча, см. [11], [12]):

где каждый из коэффициентов ${{a}_{k}}$, $k \geqslant 2$, выражается явно через ${{a}_{0}},\; \ldots ,\;{{a}_{{k - 1}}}$ по формулам рекуррентного типа.

Отысканию численных значений критического наклона (1.13) в разложении (1.12) и константы $A$ в разложении (1.14), а также разработке методов решения задачи (1.2), (1.3) посвящена обширная литература (см. [1], [2], [7], [10]–[45 ]).

1.3. Описание результатов

Для того, чтобы сформулировать результаты настоящей работы, введем в рассмотрение константы ${{\omega }_{0}}$, $\alpha $, зависящие только от параметров $\nu $, $\sigma $ уравнения (1.5), по следующим формулам:

где числа $\beta $, ${{f}_{0}}$ определены равенствами (1.9). Если показатель $\nu $ в уравнении (1.5) является рациональным числом, то будем представлять его в виде несократимой дроби

Основной результат (теорема 1) заключается в следующем параметрическом представлении решения задачи (1.5)–(1.7), (1.17) при $R = + \infty $:

где величины $\beta $, $\alpha $ заданы формулами (1.9), (1.16) соответственно, ${{T}_{1}},{{T}_{2}} = {\text{const}} > 0$, функция $\mathcal{H}(t)$ является положительной и аналитической на всем отрезке $t \in [0,1]$, включая его концы. Отметим, что значение параметра $t = 0$ соответствует пределу при $r \to + \infty $, $y \to 0$, а значение $t = 1$ – точке $r = 0$, $y = Z$. Для коэффициентов ряда Тейлора функции $\mathcal{H}(t)$ в окрестности $t = 0$ получены формулы рекуррентного типа и предложен метод вычисления тейлоровского разложения во внутренних точках отрезка $t \in [0,1]$.

Представление (1.18) развивает известный (см. [38]–[45]) подход к построению решения задачи Томаса–Ферми в параметрическом виде. Отметим также построенный в работе [46] аналитико-численный метод решения сингулярной краевой задачи для нелинейного сингулярно возмущенного уравнения, аналогичного уравнению (1.5). Новизна подхода настоящей работы заключается в использовании обнаруженного ранее (см. [47]) аналитического свойства уравнения Абеля II рода специального вида, к которому приводится исходное уравнение Эмдена–Фаулера, а именно, свойства частичного прохождения модифицированного теста Пенлеве (см. [48]) в узловой особой точке.

Основанный на представлении (1.18) (при $M = 2$, $N = - 1$, $\beta = 3$, $\alpha = (7 - \sqrt {73} ){\text{/}}2$) метод решения задачи (1.2), (1.3) сочетает в себе высокую точность получаемого результата с относительно небольшой вычислительной трудоемкостью, что трудно достижимо в рамках предлагавшихся ранее методов. В основе такого сочетания лежит тот факт, что ряд Тейлора функции $\mathcal{H}(t)$ в точке $t = 0$ демонстрирует экспоненциально быструю сходимость на всем отрезке $t \in [0,1]$ (эмпирическое значение радиуса сходимости ≈1.2; в работе также дана теоретическая оценка снизу этого радиуса), что позволяет эффективно находить решение $\Psi (r)$ и его производную $\Psi {\text{'}}(r)$ в произвольной точке луча $r \in [0, + \infty )$ с практически любой наперед заданной точностью.

Также в работе получено (см. теорему 5) представление типа (1.18) решения двухточечной краевой задачи (1.5)–(1.7) на конечном отрезке $[0,R]$ при $R \in (0, + \infty )$, $Z \in (0, + \infty )$. Такая постановка в рамках модели Томаса–Ферми отвечает случаю положительно заряженного иона.

Следствием полученных параметрических формул (см. теорему 6) является представимость решения задачи (1.5)–(1.7) при $R \in (0, + \infty ]$, $Z > 0$ и рациональном значении $\nu $ в виде сходящегося в окрестности $r = 0$ ряда

по степеням ${{r}^{{1/M}}}$, где $M$ определено формулой (1.17). Заметим, что при нарушении условия (1.17), а также при $\nu \leqslant - 1$ аналитическая структура решения $y(r)$ в окрестности $r = 0$ может иметь более сложный характер по сравнению с представлением (1.19) (см. [49, § 12.4]).

Как и в случае разложения (1.12), обобщением которого является ряд (1.19), значение коэффициента ${{b}_{M}} \in ( - \infty ,B]$ параметризует решения задачи (1.5)–(1.7) при данном $Z$ и различных $R \in (0, + \infty ]$; здесь $B = B(Z,\sigma ,\nu )$ – наклон графика в нуле решения задачи на полубесконечном интервале, т.е. при $R = + \infty $. Из (1.11) следует, что

В данной работе получена формула (см. (2.15)) для критического значения $B(1,\sigma ,\nu )$ производной решения в нулевой точке при $Z = 1$, обобщающая известную формулу Майораны (см. [40]) для величины (1.13). В этой связи отметим построенный в работах [50], [51] аналитико-численный метод нахождения константы Блазиуса, возникающей в теории погранслоя.

Теорема 2 настоящей работы дает следующее представление в окрестности $r = + \infty $ решения задачи (1.5)–(1.7), $R = + \infty $, обобщающее разложение (1.14) и уточняющее асимптотику (1.8):

где коэффициенты ${{a}_{n}}$ не зависят от $Z$, числа $\beta $, $\alpha $ заданы соответственно формулами (1.9), (1.16). Для ${{a}_{1}}(\sigma ,\nu )$ получена формула (см. (2.43)), позволяющая найти константу $A$ в разложении (1.14) с любой заданной точностью (см. (3.5)).

Хорошо известно (см. [52]), что разложение (1.14) сходится лишь при достаточно больших значениях $r$. В теоремах 3 и 4 получены новые представления, аналогичные (1.20), (1.14), для обратной к решению $y(r)$ функции $r(y)$:

причем ряд по $\mathcal{Y}$ в правой части первого равенства (1.22) в случае задачи (1.2), (1.3) демонстрирует экспоненциально быструю сходимость на всем промежутке $\mathcal{Y} \in [0,1]$, т.е. при $y \in [0,Z]$, что делает данную формулу пригодной для вычисления решения на всей области его определения. Приближая сумму этого ряда квадратичной функцией на единичном отрезке (см. п. 3.1), находим

Аналитичность функции $\mathcal{H}(t)$ в представлении (1.18) при $t = 1$ обусловлена некоторым свойством узловой особой точки уравнения Абеля II рода – вспомогательного уравнения первого порядка, к которому приводится исходное уравнение (1.5) в результате процедуры понижения порядка. Это свойство заключается в том, что семейство решений уравнения, проходящих через его узловую точку в начале координат, может быть задано общей формулой, включающей ряд по дробным степеням переменной, причем один из коэффициентов ряда является параметром данного семейства. Как показано в работе [47], такая аналитическая структура решений уравнения Абеля II рода вблизи его узловой особой точки связана с частичным прохождением некоторой модификации теста Пенлеве для данного уравнения. Отметим, что в работах [53]–[55] с помощью указанного свойства уравнения Абеля II рода были получены новые представления квазистационарных решений нелинейных параболических уравнений типа Колмогорова–Петровского–Пискунова.

1.4. Переход к автономному уравнению второго порядка

Положим

где по определению $ln( + \infty ) = + \infty $. Пользуясь обозначениями (1.9), перейдем в задаче (1.5)–(1.7) к новым переменным относительно которых исходное уравнение (1.5) с учетом определения (1.15) принимает следующий автономный вид: Тем самым, каждому решению $y(r)$ задачи (1.5)–(1.7) взаимно однозначным образом сопоставлено решение $\psi (\rho )$ уравнения (1.26), заданное на промежутке $I(R)$ вида (1.24). Обращая замену переменных (1.25), находим

Зафиксируем решение $y(r)$ задачи (1.5)–(1.7) при некоторых $0 < R \leqslant + \infty $, $Z \in (0, + \infty )$ и рассмотрим поведение соответствующего решения $\psi (\rho )$ уравнения (1.26). Из формулы подстановки (1.25) и из условия (1.6) следует асимптотика

а из условия (1.7) и соотношения (1.8) вытекает равенство

В силу автономности уравнения (1.26) его решение $\psi (\rho )$ при сдвиге вдоль оси $\rho $,

переходит в решение этого же уравнения. Из формул (1.25), (1.27) видно, что перенос (1.30) соответствует масштабному преобразованию (1.10) уравнения (1.5). Кроме того, трансляция (1.30) сохраняет величину которая, тем самым, является инвариантом масштабирования (1.10).

Исходя из условий (1.6), (1.7), докажем от противного следующее включение:

В самом деле, из предположения $\mu > 1$ и асимптотик (1.28), (1.29) вытекает наличие у функции $\psi (\rho )$ локального максимума в некоторой точке ${{\rho }_{0}} \in I(R)$, $\psi ({{\rho }_{0}}) > 1$. Тогда в этой точке должны быть выполнены условия которые в совокупности противоречат уравнению (1.26). Таким образом, предположение $\mu > 1$ не верно, и включение (1.32) выполнено.

Следующее утверждение непосредственно вытекает из предложения 5 работы [47].

Предложение 1. Для  каждого $\mu $ из полуинтервала (1.32) существует решение $\psi = {{\psi }_{\mu }}(\rho )$ задачи (1.26)–(1.29) на промежутке $I(R)$ вида (1.24) при некотором $R \in (0, + \infty ]$, имеющее заданное формулой (1.31) значение $\mu $. Функция ${{\psi }_{\mu }}(\rho )$ определена однозначно с точностью до сдвига (1.30), при этом $\mu = 1$ тогда и только тогда, когда $R = + \infty $.

Из проведенных выше рассуждений и из предложения 1 вытекает

Предложение 2. Множество решений $y(r)$ задачи (1.5)–(1.7), $Z \in (0, + \infty )$, $R \in (0, + \infty ]$ находится во взаимно однозначном соответствии с множеством решений задачи (1.26)–(1.29), (1.31), (1.32) при $\mu \in (0,1]$. Это соответствие, заданное формулами (1.25), (1.27), согласовано с преобразованиями (1.10), (1.31) так, что промежуток (1.32) параметризует классы эквивалентности решений $y(r)$ по отношению к масштабированию (1.10).

В дальнейшем нам будет удобно задавать решение $y(r)$ задачи (1.5)–(1.7), не парой $(R,Z)$, где $R \in (0, + \infty ]$, $Z \in (0, + \infty )$, а парой $(\mu ,Z)$, где $\mu \in (0,1]$, $Z \in (0, + \infty )$.

1.5. Понижение порядка уравнения и условие прохождения теста Пенлеве

Определим новую переменную $p = p(\psi )$ в виде

тогда задача (1.26), (1.28), (1.29), (1.31) при $\mu \in (0,1]$ сведется к следующей (сингулярной) задаче Коши относительно $p(\psi )$:

Определим индекс Фукса для узловой точки $\psi = 0$, $p = 0$ уравнения (1.34) с учетом (1.9), (1.15) следующим образом:

Если выполнено условие (1.17) и, тем самым, $K \in \mathbb{Q} \cap (0,1)$, то согласно [47, п. 2.3] уравнение (1.34) частично проходит модифицированный тест Пенлеве. В этом случае через $\tilde {N}$ будем обозначать знаменатель дроби, представляющей $K$:

Введем в рассмотрение переменные $q$, $z$, связанные с $p$, $\psi $ следующим образом:

тогда относительно $q(z)$ задача (1.34) принимает вид

Согласно предложению 1 работы [47] все траектории, представляющие решения задачи (1.34) на плоскости $(\psi ,p)$, приходят в начало координат $\psi = 0$, $p = 0$ с наклоном $\beta $ к горизонтали, а соответствующие решения задачи (1.38), (1.39) обладают свойством $q(0) = \beta $. Применяя теорему 1 той же работы [47], получаем

Предложение 3. Всякое решение $q(z)$ уравнения (1.38) с условием $q(0) = \beta $ является аналитической при $z = 0$ функцией вида

где функция $\mathcal{Q}(z)$ также аналитична при $z = 0$. Решения вида (1.40) образуют семейство, параметризованное величиной $C: = \mathcal{Q}(0)$, т.е. для каждого $C \in \mathbb{R}$ существует единственное решение $q(z)$ вида (1.40). При этом коэффициенты Тейлора функции $\mathcal{Q}(z)$ в точке $z = 0$ рационально выражаются через $\beta $, $M$, $\tilde {N}$, $C$.

Решения $q = {{q}_{\mu }}(z)$ задачи (1.38), (1.39) при $\mu \in (0,1]$ принадлежат указанному семейству, и формула (1.40) определяет для каждого $\mu \in (0,1]$ функцию $\mathcal{Q} = {{\mathcal{Q}}_{\mu }}(z)$. При этом отображение

является монотонным и взаимно однозначным.

Предложение 3 утверждает, таким образом, что при выполнении условия (1.17) переход в уравнении Абеля (1.34) к переменным $z$, $\mathcal{Q}$ в определенном смысле устраняет нерегулярность этого уравнения в узловой точке $\psi = 0$, $p = 0$.

2.1. Решение задачи на полубесконечном промежутке

Рассмотрим задачу (1.5)–(1.7), $R = + \infty $ при некотором $Z > 0$. Согласно предложениям 1, 2 ее решение $y(r)$ выражается по формулам (1.27) через решение $\psi = {{\psi }_{1}}(\rho )$ задачи (1.26), (1.28), (1.29), (1.31) при $\mu = 1$.

2.1.1. Параметрическое представление решения и его производной. Предположим, что для показателя $\nu $ в уравнении (1.5) выполнено условие рациональности (1.17). Тогда теорема 2 работы [47] утверждает, что функция $\psi = {{\psi }_{1}}(\rho )$ имеет следующее параметрическое представление:

где $\tilde {N}$ определено формулой (1.36), функция $H(t)$ является аналитической на отрезке $t \in [0,1]$ и выражается по формулам через решение $\mathcal{Q}(t)$ следующей задачи Коши: причем $\mathcal{Q}(t)$ – убывающая отрицательная аналитическая функция на отрезке $t \in [0,1]$, а входящая в знаменатель правой части уравнения (2.3) функция является неотрицательной и возрастающей на единичном отрезке, $g(0) = 0$, $g(1) = 1$. Здесь переменная $t$ связана с $z$ из (1.37) равенством и $g(t) = {{\beta }^{{ - 1}}}q(t)$, где $q(t)$ имеет вид (1.40).

Подставляя параметризацию (2.1) в формулы (1.27), с учетом равенства

вытекающего из (1.35), (1.36), получаем представление (1.18), где ${{T}_{1}} = {\text{const}} > 0$,

Выберем константу ${{T}_{1}}$ в (2.1) таким образом, чтобы удовлетворить краевому условию (1.6). Вводя обозначение

подставляя $t = 1$ в (1.18) и приравнивая $Z$ полученное выражение для $y$ при $r = 0$, запишем условие (1.6) в следующем виде: откуда с учетом (2.9) найдем где показатель $K$ задан равенством (1.35).

На основании формул (2.2), (2.8), (2.5) найдем производные функций $r(t)$, $y(t)$ в параметрическом представлении (1.18):

Используя полученные выражения (2.12), (2.13), вычислим производную решения $dy{\text{/}}dr$ в зависимости от параметра $t \in (0,1]$:

где константы ${{T}_{1}}$, ${{T}_{2}}$ заданы равенствами (2.11). В частности, при $t = 1$ из (2.14) с учетом (2.10), (2.11) получаем выражение для критического наклона графика $y(r)$ при $r = 0$:

Отметим, что в силу (2.12), (2.13), (1.18) обе производные $dy{\text{/}}dt$ и $dr{\text{/}}dt$ имеют нули порядка $(M - 1)$ при $t = 0$ и не обращаются в нуль в других точках единичного отрезка.

2.1.2. Вычисление тейлоровских коэффициентов и оценка радиуса сходимости. Найдем коэффициенты разложения введенных выше (2.2)–(2.4), (2.8) функций $\mathcal{Q}(t)$, $\mathcal{H}(t)$ по степеням переменной $t$. Используя обозначение (2.5), перепишем уравнение (2.3) в виде

и подставим в получившуюся систему (2.5), (2.16) формальные разложения

Приравнивая коэффициенты при ${{t}^{n}}$, $n \geqslant 2$, в обеих частях равенства (2.16), получаем

Перепишем последнее равенство в виде где правая часть зависит лишь от величин ${{\mathcal{Q}}_{j}}$, ${{g}_{j}}$, $j = 0,\; \ldots ,\;(n - 1)$, и не зависит от ${{\mathcal{Q}}_{n}}$. Найдем множитель при ${{\mathcal{Q}}_{n}}$ в левой части (2.19):

Решая (2.19) как линейное уравнение относительно ${{\mathcal{Q}}_{n}}$ последовательно при $n = 2,3,\; \ldots $, с учетом (2.21) получаем

Формулы (2.18), (2.20), (2.22) позволяют вычислить каждый из коэффициентов разложения Тейлора (2.17) функции $\mathcal{Q}(t)$ за конечное число арифметических операций.

Найдем коэффициенты ряда Тейлора

функции $h(t)$ на основании (2.2), (2.5), (2.7): Подставляя полученные выше разложения (2.17), (2.18) в выражение (2.24), находим где ${{G}_{n}}$ – коэффициент при ${{t}^{n}}$ в числителе дроби (2.24), откуда получаем

Для тейлоровских коэффициентов ${{H}_{n}}$ функции

имеем

Ряд Тейлора функции (2.8) получаем, раскрывая скобки и приводя подобные члены в правой части формулы

где коэффициент при ${{t}^{n}}$ представляет собой многочлен от ${{H}_{1}},\; \ldots ,\;{{H}_{n}}$.

Оценим радиус сходимости рядов (2.17), (2.18), (2.23), предполагая дополнительно, что значение параметра $\sigma $ уравнения (1.5) лежит в диапазоне

Сначала покажем, что решение $p = {{p}_{1}}(\psi )$ задачи (1.34), $\mu = 1$, является аналитической функцией в круге $\mathcal{O}: = \left\{ {\left| {\psi - 1} \right| < 1} \right\}$. Согласно предложению 4 работы [47] достаточно показать, что свободный член $F(\psi ) = {{f}_{0}}\psi (1 - {{\psi }^{\sigma }})$ этого уравнения обладает следующими свойствами: во-первых, $F(\psi )$ является голоморфной внутри круга $\mathcal{O}$ функцией, не обращающейся в этом круге в нуль нигде, кроме его центра $\psi = 1$, и, во-вторых, коэффициенты разложения функции $F(\psi ){\text{/}}\psi $ в точке $\psi = 1$ по степеням переменной $\tau : = (1 - \psi )$ образуют убывающую последовательность положительных чисел. Первое из этих условий, очевидно, выполнено; проверим второе условие:

что с учетом условия (2.30) дает требуемый результат.

При замене переменной $\psi $ на t при помощи подстановок (1.37), (2.6) круг $\mathcal{O}$ сходимости функции ${{p}_{1}}(\psi )$ конформно отображается на выпуклую “каплевидную” область $\mathcal{D}$ аналитичности функций $\mathcal{Q}(t)$, $q(t)$, $h(t)$, $H(t)$, $\mathcal{H}(t)$. Область $\mathcal{D}$ имеет угловую точку с раствором угла $(\pi \sigma {\text{/}}\tilde {N})$ в точке $t = 1$. Расстояние от начала координат $t = 0$ до границы $\partial{ \mathcal{D}}$ равно

Таким образом, получен следующий результат.

Предложение 4. При выполнении условия (2.30) функции $\mathcal{Q}(t)$, $q(t)$, $h(t)$, $H(t)$, $\mathcal{H}(t)$ голоморфны в области $\mathcal{D}$ и величина (2.31) ограничивает снизу радиус сходимости рядов (2.17), (2.18), (2.23).

Получив значение функции $\mathcal{Q}(t)$ в некоторой внутренней точке единичного интервала $t \in (0,1)$, можно найти в этой точке тейлоровские разложения функций $\mathcal{Q}(t)$, $g(t)$ методом неопределенных коэффициентов из уравнений (2.5), (2.16), и затем вычислить соответствующие разложения функций $h(t)$, $H(t)$, $\mathcal{H}(t)$ по формулам (2.24), (2.27)–(2.29).

Суммируем результаты, полученные в пп. 2.1.1, 2.1.2, в виде следующей теоремы.

Теорема 1. Рассмотрим задачу (1.5)–(1.7), $R = + \infty $ при $Z \in (0, + \infty )$ и при условии (1.17). Тогда для решения $y(r)$ этой задачи и его первой производной $dy{\text{/}}dr$ справедливо параметрическое представление (1.18), (2.11), (2.14), т.е.

где константы $\beta $, ${{f}_{0}}$, ${{\omega }_{0}}$, $\alpha $, $K$, ${{\mathcal{H}}_{e}}$ зависят только от $\sigma $, $\nu $ и задаются формулами (1.9), (1.15), (1.16), (1.35), (2.10), функция $\mathcal{H}(t)$ является положительной аналитической на отрезке $t \in [0,1]$ и выражается через решение $\mathcal{Q}(t)$ задачи (2.3), (2.4) по формулам (2.24), (2.27), (2.29). Разложения (2.17), (2.29) функций $\mathcal{Q}(t)$, $\mathcal{H}(t)$ в ряды Тейлора в точке $t = 0$ можно получить по формулам (2.18), (2.20), (2.22), (2.25)–(2.29). Если $\sigma \in (0,1]$, то радиус сходимости этих рядов ограничен снизу величиной (2.31).

2.1.3. Исключение параметра t. Рассмотрим решение $y(r)$ задачи (1.5)–(1.7), (1.17), $R = + \infty $ при $Z \in (0, + \infty )$ в виде (2.32), (2.33) и исключим параметр $t$ из этого представления. Возводя в степень $\sigma {\text{/}}\tilde {N}$ обе части равенства $\psi = (\mathop {{{f}_{0}}}\nolimits^{ - 1/\sigma } {{r}^{\beta }}y)$, т.е. первой формулы (1.25), запишем результат с учетом (2.1), (2.7) в следующем виде:

Выразим параметр $t$ через $r$ из равенства (2.32). Разделив обе части (2.32) на ${{T}_{1}}$ и возведя их в степень $\alpha $, введем обозначения соответственно $\tilde {r}$ и $\xi $ для преобразованных левой и правой частей:

где функция $\xi = \xi (t)$ является аналитической по $t$ на единичном полуинтервале.

Проверим, что $d\xi {\text{/}}dt \ne 0$ при $t \in [0,1)$. В самом деле, в силу (2.13) с учетом (2.5) имеем

таким образом, задаваемое формулой (2.36) соответствие $t \mapsto \tilde {r} = \xi (t)$ монотонно и взаимно однозначно отображает полуинтервал $t \in [0,1)$ на луч $\tilde {r} \in [0, + \infty )$. Тогда существует обратная к $\xi (t)$ возрастающая положительная аналитическая функция $\theta = \theta (\tilde {r})$, отображающая полуось $\tilde {r} \in [0, + \infty )$ на промежуток $t \in [0,1)$, и параметр $t$ в представлении (2.32) с учетом связи (2.36) между $\tilde {r}$ и $r$ выражается через $r$ в виде

Вычислим разложение функции $\xi (t)$ по степеням переменной $t$. В окрестности $t = 0$ справедливы формулы

перемножая которые, находим где коэффициенты ${{\xi }_{n}}$, ${{\theta }_{n}}$ полиномиально зависят от известных величин $\alpha $, $M$, ${{H}_{1}}$, $ \ldots $, ${{H}_{{n - 1}}}$.

Подставляя выражение (2.37) для $t$ в формулу (2.35), получаем следующие представления функции $y = y(r)$:

причем в окрестности $r = + \infty $ представление (2.41) с учетом (2.40) принимает вид (1.20):

Таким образом, получена следующая

Теорема 2. В условиях теоремы 1 параметр $t \in (0,1]$ выражается через $r$ в виде (2.37). Для решения $y(r)$ задачи (1.5)–(1.7), $R = + \infty $, справедливо представление (2.41), принимающее в окрестности $r = + \infty $ вид (2.42), (2.43).

Выразим теперь параметр $t$ через $y$ из равенства (2.33). Возводя в степень $ - \alpha {\text{/}}\beta $ обе части этого равенства и вводя по аналогии с (2.36) обозначения $\tilde {y}$, $\Xi $, находим

В силу (2.44), (2.12), (2.5), (1.18), при $t \in [0,1]$ имеем следовательно, производная $d\Xi {\text{/}}dt$ положительна при $t \in (0,1]$ и имеет нуль порядка $(M - 1)$ при $t = 0$. Таким образом, заданное равенством (2.44) отображение $t \mapsto \tilde {y} = \Xi (t)$ определяет взаимно однозначное (в силу монотонности) соответствие между отрезками $t \in [0,1]$ и $\tilde {y} \in [0,\mathcal{H}_{e}^{\alpha }]$, и существует обратная к $\Xi (t)$ возрастающая аналитическая на полуинтервале $\tilde {y} \in [0,\mathcal{H}_{e}^{\alpha })$ функция $\Theta = \Theta (\tilde {y})$, имеющая (как функция комплексной переменной $\tilde {y}$) $M$ – листное ветвление при $\tilde {y} = \mathcal{H}_{e}^{\alpha }$. С помощью функции $\Theta (\tilde {y})$ параметр $t$ выражается через $y$ в виде

Для функции $\Xi (t)$ в окрестности $t = 0$ справедливо разложение (2.38), откуда находим разложение Тейлора в точке $\tilde {y} = 0$ обратной функции:

Подставляя $t$ в виде (2.45) в формулу (2.35) и разрешая ее относительно $r$, получаем

откуда с учетом разложения (2.46) вытекает представление (1.21) в окрестности $r = 0$, т.е.

Таким образом, справедлива

Теорема 3. В условиях теоремы 1 параметр $t \in (0,1]$ выражается через $y$ в виде (2.45), причем функция $t = t(y)$ имеет в комплексной области $M$ – листное ветвление при $y = Z$. Для обратной к решению функции $r = r(y)$ верно представление (2.44), (2.47), принимающее вид (2.48) в окрестности $y = 0$.

Круг сходимости ряда в правой части (2.48) заведомо не покрывает точку ветвления $y = Z$ аналитической (по $y$) функции (2.45), соответствующую точке $M$ – листного ветвления $\tilde {y} = \mathcal{H}_{e}^{\alpha }$ функции $\Theta (\tilde {y})$. Устраним эту особую точку, введя, во-первых, переменную $\mathcal{Y} \in [0,1]$ по формулам

и, во-вторых, аналитическую на единичном отрезке функцию взаимно однозначно отображающую отрезок $\mathcal{Y} \in [0,1]$ на промежуток изменения параметра $t \in [0,1]$. Тогда формулу обращения (2.45) с учетом (2.49), (2.50) запишем следующим образом: Отметим, что в силу (2.46), (2.50) разложение $\epsilon (\mathcal{Y})$ в окрестности $\mathcal{Y} = 0$ имеет вид

Подставляя (2.50) в (2.47), получаем зависимость $r$ от $y$ в виде

Возводя обе части равенства (2.54) в степень $M$ и подставляя разложение (2.52), после раскрытия скобок получаем совпадающее с (1.22) представление в окрестности $y = 0$ для обратной к решению функции

Таким образом, получена

Теорема 4. В условиях теоремы 1 параметр $t \in (0,1]$ выражается через $y$ в виде (2.51), (2.49), где функция $\epsilon (\mathcal{Y})$ является аналитической функцией на замкнутом отрезке $\mathcal{Y} \in [0,1]$. Для обратной к решению $y(r)$ функции $r = r(y)$ верно представление (2.53), принимающее вид (2.54) в окрестности $y = 0$.

2.2. Решение задачи на отрезке

Получим решение задачи (1.5)–(1.7), (1.17) при $Z,R \in (0, + \infty )$ в предположении, что известно соответствующее (см. предложение 2) значение параметра $\mu \in (0,1)$.

Теорема 3 работы [47] утверждает, что решение $\psi = {{\psi }_{\mu }}(\rho )$ задачи (1.26)–(1.29), (1.31) при заданном $\mu $ допускает следующее параметрическое представление:

где функция $S(u)$ является аналитической на полуинтервале $u \in ({{u}_{1}},1]$ и выражается в виде через решение $\mathcal{P}(u)$ следующей задачи Коши: Также в работе [47] показано, что $\mathcal{P}(u) < 0$ при $u \in ({{u}_{1}},1]$, а знаменатель правой части уравнения (2.57) не обращается в нуль при $u \ne 0$.

Отметим, что если справедливо включение

то уравнение (2.57) удовлетворяет условию теоремы Коши (см. [56, § 3]) при $u = - 1{\text{/}}\tilde {N} = {{u}_{1}}$, следовательно, функция $\mathcal{P}(u)$ в этом случае является аналитической на всем отрезке $u \in [{{u}_{1}},1]$, включая его левый конец.

Подставляя параметризацию (2.55) в формулы (1.27) и полагая

с учетом (2.7) получаем

Значение постоянной ${{U}_{1}}$ найдем из сопоставления (2.62) и краевого условия (1.6), т.е. $y = Z$ при $u = 1$:

где $K$ задано формулой (1.35).

Вычислим производную $dy{\text{/}}dr$ как функцию параметра $u$. Из формул (2.56), (2.60)–(2.62) с учетом (2.7) находим

откуда получаем следующее выражение для производной:

Таким образом, доказана

Теорема 5. Решение $y(r)$ задачи (1.5)–(1.7), (1.17) при $Z,R \in (0, + \infty )$ допускает следующее параметрическое представление:

где числа $K$, $M$, $\tilde {N}$ заданы формулой (1.35), ${{u}_{1}} = ( - 1{\text{/}}\tilde {N})$, функция $\mathcal{S}(u)$ является аналитической на полуинтервале $u \in ({{u}_{1}},1]$, ${{\mathcal{S}}_{e}} = \mathcal{S}(1)$, и справедливо представление (2.56), (2.60) функции $\mathcal{S}(u)$ через решение $\mathcal{P}(u)$ задачи (2.57), (2.58). Значение параметра $\mu \in (0,1)$, входящего в уравнение (2.57), связано с решением $y(r)$ формулами (1.31), (1.25).

Если выполнено условие (2.59), то функции $\mathcal{P}(u)$, $\mathcal{S}(u)$ являются аналитическими на всем отрезке $u \in [{{u}_{1}},1]$, включая его левый конец $u = {{u}_{1}}$.

Из формулы (2.68) при $u = 1$ получаем значение наклона графика решения $y(r)$ в концевой точке $r = 0$ в следующем виде, аналогичном (2.15):

Из формул (2.64), (2.65) и (2.61), (2.62) следует, что каждая из производных $dr{\text{/}}du$ и $dy{\text{/}}du$ представляется в виде произведения множителя ${{(1 - u)}^{{M - 1}}}{{(\tilde {N}u + 1)}^{{1/\sigma - 1}}}$, задающего кратность нулей в концах отрезка параметризации $u \in [{{u}_{1}},1]$, и некоторой функции, не обращающейся в нуль нигде на этом отрезке.

Тейлоровские разложения функций $\mathcal{P}(u)$, $s(u)$, $\mathcal{S}(u)$ в точке $u = 0$ и других точках промежутка $u \in ({{u}_{1}},1]$ могут быть вычислены, исходя из уравнения (2.57) по схеме, изложенной в п. 2.1.2.

2.3. Поведение решения вблизи начала координат

Покажем, что для решения $y(r)$ уравнения (1.5), (1.17), $Z \in (0, + \infty )$, с условием (1.6) справедливо разложение вида (1.19) в окрестности $r = 0$. Для этого воспользуемся параметрическим представлением этого решения с параметром $z$, определяемым подстановкой (1.37).

Из второй формулы (1.25), подстановок (1.33), (1.37) и равенства (2.7) находим

где функция $G(z)$, имеющая с учетом (1.40) вид является аналитической при $z = 0$ и определена с точностью до аддитивной константы. Из равенств (2.70), (1.37), (2.7) и первой формулы (1.25) получаем искомую параметризацию в окрестности $z = 0$ ($r = 0$, $y = Z$), охватывающую представления (2.32), (2.33) и (2.66), (2.67). Отметим, что параметрическое представление (2.72) в частном случае задачи (1.2), (1.3) было получено в работе [40].

Используя обозначение $C: = \mathcal{Q}(0)$ (см. предложение 3), вычислим начальные члены тейлоровских разложений величин $G$, $\mathcal{G}$, ${{r}^{{1/M}}}$, $y$ по степеням $z$, исходя из формул (2.71), (2.72):

где символ Кронекера $\delta _{M}^{1}$ равен единице при $M = 1$ и нулю в противном случае. Согласно предложению 3 коэффициенты рядов внутри скобок в правых частях равенств (2.73)–(2.75) рационально выражаются через $\beta $, $M$, $\tilde {N}$, C. Сопоставляя равенство (2.75) и условие (1.6), находим

Выражая $z$ через $r$ из равенства (2.74), получаем в окрестности $r = 0$ разложения

где коэффициенты рядов в правых частях также являются рациональными функциями от $\beta $, $M$, $\tilde {N}$, C. Подставляя равенства (2.77), (2.78) в формулу (2.75), с учетом (2.76) получаем в окрестности $r = 0$ представление (1.19): причем в силу предложения 3 каждому $C \in \mathbb{R}$ соответствует решение $y(r)$ вида (2.79).

Решение сингулярной задачи Коши для уравнения (1.5) с условием (1.6) и заданной производной при $r = 0$ существует и единственно по теореме Каратеодори (см. [57, гл. 2, § 1]). Получаем следующее уточнение этого результата.

Теорема 6. Для всякого $b \in \mathbb{R}$ задача Коши для уравнения (1.5), (1.17) с условиями (1.6), $dy{\text{/}}dr(0) = b$ однозначным образом разрешима на некотором отрезке $[0,\varepsilon ]$, $\varepsilon > 0$, при этом ее решение $y(r)$ представимо в окрестности $r = 0$ сходящимся рядом (2.79), где ${{b}_{M}} = b$.

Коэффициенты ${{b}_{n}}$ разложения $($2.79) при заданном ${{b}_{M}} = b$ можно найти методом неопределенных коэффициентов, подставив его в уравнение (1.5). Кроме того, в силу (2.76) и указанных выше свойств коэффициентов разложений (2.77), (2.78), при каждом $n \geqslant M$ величина ${{b}_{n}}\mathop {{{f}_{0}}}\nolimits^{n/\tilde {N}} {\text{/}}{{Z}^{{1 + (\sigma n/\tilde {N})}}}$ рационально выражается через $\beta $, $M$, $\tilde {N}$ и $C = {{b}_{M}}\mathop {{{f}_{0}}}\nolimits^K {\text{/}}(\beta {{Z}^{{1 + (1/\beta )}}})$.

3.1. Решение задачи Томаса–Ферми на луче

Представление (2.32)–(2.34) решения $\Psi (r)$ задачи (1.2), (1.3) при $R = + \infty $ имеет вид

где функции $\mathcal{Q}(t)$, $\mathcal{H}(t)$ являются аналитическими на отрезке $t \in [0,1]$, при этом $\mathcal{Q}(t)$ есть решение задачи (2.3), (2.4), т.е. и $\mathcal{H}(t)$ выражается через $\mathcal{Q}(t)$ по формулам (2.2), (2.8). Отметим, что уравнение (3.2) было получено Э. Майораной (см. [40]) непосредственно из уравнения (1.2).

Разложения функций $\mathcal{Q}(t)$ и $\mathcal{H}(t)$ в ряды Тейлора в окрестности $t = 0$ могут быть вычислены по формулам (2.17), (2.18), (2.20), (2.22) и (2.25)–(2.29) соответственно, т.е.

причем радиус сходимости этих рядов ограничен снизу величиной (2.31): $\sqrt[6]{2} - 1 \approx 0.12$. Эмпирический радиус круга сходимости при этом составляет ≈1.2, т.е. этот круг полностью покрывает область изменения параметра $t$ – единичный отрезок $t \in [0,1]$.

Представление (2.54) обратной к решению функции принимает вид

причем ряды в правой части (3.4) демонстрируют еще более быструю сходимость, чем тейлоровские разложения функций $\mathcal{Q}(t)$, $\mathcal{H}(t)$, с эмпирическим радиусом сходимости ≈1.8 (см. табл. 1 и фиг. 2). При помощи значений коэффициентов ${{\epsilon }_{n}}$, представленных в табл. 1, значение функции $\epsilon (\tilde {\Psi })$ может быть найдено с точностью до 12 знаков после запятой. Графики функций $\mathcal{Q}(t)$, $\mathcal{H}(t)$, $\epsilon (\tilde {\Psi })$ представлены на фиг. 3 и 4.

Представление (2.42) принимает вид (1.14) с константой $A$, заданной формулой (2.43):

где в силу рекуррентных соотношений (3.3) имеем

Сформулируем полученный результат.

Теорема 7. Решение $\Psi (r)$ задачи (1.2), (1.3) при $R = + \infty $ допускает параметрическое представление (3.1). Коэффициенты Тейлора входящих в него функций $\mathcal{H}(t)$, $\mathcal{Q}(t)$ могут быть вычислены по формулам (3.3), (2.8). Обратная к решению функция $r = r(\Psi )$ допускает представление (3.4), где ряды в правой части сходятся в окрестности $\tilde {\Psi } = \Psi = 0$. Для константы $A$ в разложении (1.14) справедлива формула (3.5).

Подчеркнем отличие формул (3.5), (3.6), основанных на суммировании экспоненциально сходящегося числового ряда, каждый член которого вычисляется явно за конечное число операций из соотношений рекуррентного типа (3.3), от аналогичной формулы (см. [42, § 2.12]), опирающейся на параметризацию вида (2.72) и требующей численного разрешения неопределенности вида $(0 \cdot \infty )$, где бесконечно большой множитель является результатом численного интегрирования неограниченно растущей функции.

Отметим, что выражение для константы (1.13), даваемое формулой (2.15) при $Z = 1$,

было получено Э. Майораной (см. [40], [42, § 2.12]).

Вычислить решение $\Psi (r)$ и его производную $d\Psi {\text{/}}dr$ в заданной точке ${{r}_{0}} \in (0, + \infty )$ с помощью представления (3.1) можно, найдя соответствующее значение параметра $t = {{t}_{0}}$ и подставив его в формулу $\Psi = \Psi ({{t}_{0}})$. Для обращения функции $r = r(t)$ методом Ньютона нужно учесть, что в силу формулы (2.13) эта функция имеет кратный нуль порядка $M = 2$ при $t = 1$, и ее производная не обращается в нуль более нигде на отрезке $t \in [0,1]$. Таким образом, эффективность метода Ньютона обеспечена для задачи обращения не самой функции $r(t)$, а функции $\sqrt {r(t)} $ (или функции ${{r}^{{1/M}}}(t)$ в общем случае представления (2.32)), при этом начальное приближение для ${{t}_{0}}$ можно выбрать на отрезке $[0,1]$ произвольно.

Принимая в (3.4) приближение $\epsilon (\tilde {\Psi }) \approx \tilde {\Psi }$, получаем формулу (1.23). Аппроксимируя функции $\epsilon (\tilde {\Psi })$, $\mathcal{Q}(t)$, $\mathcal{H}(t)$ тем или иным способом, можно получить из представлений (3.1), (3.4) другие эффективные формулы приближенного решения задачи в явном аналитическом виде.

3.2. Решение задачи Томаса–Ферми на отрезке

Представление (2.66)–(2.68) решения $\Psi (r)$ задачи (1.2), (1.3) при $Z,R \in (0, + \infty )$ имеет вид

где функция $\mathcal{P}(u)$ является решением задачи (2.57), (2.58), т.е. аналитическая при $u \in ({{u}_{1}},1]$ функция $\mathcal{S}(u)$ выражается через $\mathcal{P}(u)$ в виде (2.56), т.е.

Численное решение $\Psi (r)$ задачи (1.2), (1.3) при $R,Z \in (0,\infty )$ с помощью представления (3.7) проводится аналогично случаю $R = + \infty $, рассмотренному в п. 3.1. Круг сходимости ряда Тейлора функции $S(u)$ с центром $u = 0$ может не покрывать всего отрезка изменения параметра $u \in [ - 1{\text{/}}3,1]$ – в этом случае можно построить тейлоровские разложения функций $\mathcal{P}(u)$, $S(u)$ в нескольких других точках этого отрезка с помощью уравнения (3.8).

Значение параметра $\mu $, соответствующее заданным $R,Z \in (0, + \infty )$, может быть приблизительно установлено c помощью интерполяции согласно табл. 2 (вычисления проведены на основе формул (3.7)–(3.9), (2.69)).