Вычисление ньютоновского потенциала гравитирующей конфигурации с поверхностью, близкой к сфероиду, с помощью символьных и численных методов Текст научной статьи по специальности «Физика»

ф

ф

Математическое моделирование

УДК 519.6,517.9

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

Е. В. Беспалько, С. А. Михеев, В. П. Цветков, А. Н. Цирулев, И. В. Пузынин

В работе получено точное аналитическое представление ньютоновского гравитационного потенциала неоднородной возмущённой эллипсоидальной конфигурации во внутренней точке Фт. С помощью составленной программы в системе символьной математики MAPLE и численных вычислений регуляризованным аналогом метода Ньютона получено выражение для потенциала Ф;п в виде полинома от координат. В квадруполь-ном приближении получено выражение для потенциала во внешней точке Ф0ut.

Ключевые слова: ньютоновский гравитационный потенциал, возмущённая эллипсоидальная поверхность, регуляризованный аналог метода Ньютона, система символьной математики MAPLE.

Известно, что задача о потенциалах однородного эллипсоида возникла первоначально как предмет теории тяготения. Её решение позднее нашло применение в различных физических приложениях. Например, формулы для ньютоновского гравитационного потенциала Ф используются в гидродинамике в задачах о потенциальном течении несжимаемой идеальной жидкости вокруг эллипсоида, о силе сопротивления, действующей на медленно движущийся в вязкой жидкости эллипсоид, и т. д. Причём в задачах о конфигурациях гравитирующих систем необходимо аналитическое представление внутреннего потенциала, а в задачах, связанных с описанием движения одной системы относительно другой, — внешнего потенциала, так как потенциал входит в уравнения, определяющие конфигурацию, и в уравнения, описывающие динамику гравитирующих масс.

В теории ньютоновского гравитационного потенциала возникает три типа задач по трём типам симметрии. Наиболее простой случай, когда конфигурация не вращается; тогда она будет иметь сферически симметричную форму. Этот случай наиболее хорошо изучен [1]. Но часто необходимо учитывать вращение, тогда конфигурация будет иметь сферически несимметричную форму [2]. В этом случае мы приходим к необходимости использовать более сложные поверхности, в частности эллипсоидальные, граница которых представляет собой эллипсоид [3,4]. Изучением эллипсоидальных фигур равновесия занимались Маклорен, Якоби, Ляпунов и многие другие. Ими были получены точные аналитические представления потенциала Ф в случае, когда плотность представлена в виде степенных функций координат. Внутренний потенциал является тогда многочленом по степеням координат, коэффициенты которого представляют собой эллиптические интегралы.

Статья поступила в редакцию 15 июня 2005 г.

Кафедра общей математики и математической физики Тверской государственный университет Россия, 170100, г. Тверь, ул. Желябова, 33

Введение

е-

ф

ф

ф

Однако в реальных случаях для неоднородных конфигураций возникает задача, в которой поверхность конфигурации представляет из себя более сложную структуру. В работе [5] предложен метод аппроксимации поверхности псевдоповерхностью, а именно возмущённой эллипсоидальной поверхностью, параметры которой определяются из условия минимума квадрата плотности на этой поверхности. Для этих целей эффективно может быть использован метод разложения в ряд Бурмана-Лагранжа по малому параметру.

Задача вычисления ньютоновского потенциала возмущённых эллипсоидальных конфигураций актуальна в связи с изучением нелинейных эффектов в астрофизике, геофизике и так далее, она представляет собой трудную математическую задачу [6].

Ввиду чрезвычайной сложности как аналитических, так и численных расчётов возникает необходимость использования компьютерных методов. С этой целью нами был выбран пакет символьной и численной математики MAPLE, в частности, опробованный в работе [7].

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

Задачу вычисления ньютоновского гравитационного потенциала Ф можно разбить на две задачи: отыскание его во внешней точке Фоиь и задачу отыскания потенциала во внутренней точке Ф;п [1].

Ньютоновский гравитационный потенциал Ф при этом удовлетворяет уравнениям:

ЛФопь = 4 nGp, ЛФ-п = 0.

Решением этих уравнений является

г p(r') dV-

ф = -g/p^v- , (1)

J |r - r

D

где D — область, занимаемая конфигурацией.

1. Выбор параметров возмущённой эллипсоидальной поверхности и метод рядов Бурмана-Лагранжа

Граница конфигурации Е находится из условия равенства плотности на границе нулю:

p(x,y,z) = 0, (x,y,z) G Е. (2)

Гравитационный потенциал Ф явно зависит от формы границы (2). В явном виде аналитически его удаётся вычислить только для простейших форм поверхности (шар, эллипсоид), поэтому в общем случае точную границу конфигурации Е мы заменим псевдограницей 5D, форма которой зависит от неизвестных пока параметров Zjk.

Потенциал Ф будет удовлетворять на границе условиям: Фои^Е = ФтЕ и (^out)z = (^Ф-п)Е.

Следуя работам [8], выберем поверхность Е в виде возмущённой эллипсоидальной поверхности

SD : xl + x2 + x3 Ziiiikx\x{xk3 = 1, (3)

i,j,k

X У 7

где х\ = ) Х2 = ^, жз = 7, а1' аз — полуоси эллипсоида вращения, которые наряду с Zijk параметризуют ёП, Ь — максимальная степень многочлена по координатам хх, Х2, хз. Параметр сплюснутости е = •

Условия близости (2) и (3) можно сформулировать введением функционала Л:

Л = 1

2 J Р2 dn. (4)

SD

Очевидно, параметр nsD = Л2 будет представлять меру погрешности в наших уравнениях при замене точной поверхности конфигурации Е на SD. Под псевдоповерхностью ÓD будем понимать поверхность, на которой среднее значение квадрата плотности не превосходит достаточно малой величины по сравнению с единицей. Эта величина определяет точность решения поставленной задачи.

Условие минимума Л приводит к системе алгебраических уравнений относительно ai, a3 и Zijk:

дЛ

Wjk = = o, dZijk

д

Wi = ai д^Л^к = 0) = 0, (5)

д

W2 = аз дО^Л^к = 0) = 0.

Представим плотность конфигурации р в виде полинома степени Р

P

Р(Р)= Е Ра,Ъ,сХаАх3. (6)

a,b,c

Если выбрать Р достаточно большим, то с высокой степенью точности выражение (6) будет аппроксимировать плотность реальной конфигурации.

После перехода к сферическим координатам R, в, ф: Хк = Rak, ¿i = sin в cos ф, а,2 = sin в sin ф, <5з = cos в выражение для возмущённой эллипсоидальной поверхности (3) примет вид:

L Ri+j+k

R =1 — Е Zijk ^—rajk = 1 + W (R, ak, Zijk). (7)

rrf R +1

ijk

Дальнейшие вычисления будут основаны на использовании варианта теоремы Лагранжа [6].

Теорема. [6]. Пусть f (z) и W(z) аналитические функции z на контуре C, окружающем точку а, и внутри него и для всех z на C выполняется неравенство: k\W(z)|c ^ |z — a\c. Тогда уравнение z = а + kW(z) имеет один корень z = £ внутри C и f (£) разлагается в степенной абсолютно сходящийся ряд

Ж Ks ds-1

f (£) = f (a) + Е ^das-r [f'(a)W(a)s] . (8)

s=1

Ряд (8) получил в литературе название ряда Бурмана-Лагранжа.

В нашем случае £ = R, f (£) = £h+2, a =1, W(a) = W(a, Zijk, ak). Тогда получим

Ж 1 ds-i г т

Rh+2 = ah+2 + (h + 2)E^T^zy ah+iW(a)s . (9a)

s=i

Более детально (9а) будет выглядеть следующим образом:

вЬ (_1)в . нв-1 aí+j+k+h+1

Е^+Ь = ан+2 + (Г + 2) £ £ ^Zгok(фЛ¿-г а , (9Ь)

з=1 гЛк ^ '

где Zijk(в) —многочлен от Zijk степени, не превосходящей в.

Рассмотрим более подробно а +а+1+т + . После замены а = у + 1 (у близко

к нулю) получаем ¿ут-г -—(Г+ж)!-

Для определения явного вида Фв понадобится соотношение [9]:

= -1Ф

2s ФS '

y=0 2

Фs = S (a,b)

(1 + y)b l(a)

(1 + y )a+1

= П(Ь + 1 - 2r),

y=0 r=1

где a = s — 1, b = i + j + k + h + 1. Из него и (9b) следует:

s=œ sL (_i)sh©(s)

Rh ^E -^ï-Zijk (s)à1àJ2à3fc Фs(i + j + k + h — 1), (9c)

s=0 ij 3

где ©(s) = 0 при s = 0 и ©(s) = 1 при s > 0. Фs = 1 при s = 0 и s = 1. Квадрат плотности будет выглядеть следующим образом:

P2(P) = Е Е PabcPa'b'CХ^'x2+b'Х^'. (10)

abc a'b'c'

Так как мы перешли к сферическим координатам, то, учитывая (9c), выражение для квадрата плотности примет вид:

<» Р Р вЬ -.Лвг^И

Р2(Р) = £ £ £ £ ( .Г1 РаЬсРа'ЪС X

в=0 аЪс а'Ъ'е'

X ZгJk(в)aa+a'+гaЬ2+ь'+jа1+с'+кФв(Г + г + 3 + к - 1), (11)

где Г1 = а + а' + Ь + Ь' + с + с'.

Учитывая (4) и (11), Л можно представить в явном виде, удерживая в разложении Л члены Zijk степени не выше вт, где вт — максимально учитываемая степень по в :

P sm sL

а, л (—1)sh0(s) ^(sm) = -2ss!-PabcPa'b'c' Zij3 X

abc,a'b'c' s=0 ij3

x Фs(hl + i + j + k — 1)QA1;A2,A3 , (12)

QU^a = J «A1 «A2«A3 d^, (13)

где Ai = a + a' + i, A2 = b + b' + j, Ai = c + c' + k. Интегралы (13) могут быть упрощены:

- (a + a' + i — 1)!!(b + b' + j — 1)!!(c + c' + k — 1)!!

Qa+a'+i,b+b' +j,c+c'+3 = 2n (hi + i + j + k — 1)!! . (14)

Трёхмерные массивы неизвестных в системе (5) Zijk обозначим через {ут} (т = 1,2, ..N1).

При Р = 6, Ь = 2 имеем N1 = 5. Связь между переменными будет выглядеть следующим образом:

У1 = у, У2 = е, уз = Zoo2, у 4 = Z020, У5 = Z200,

где I — характерный масштаб распределения масс гравитирующей системы. В этом случае систему уравнений (5) удобно записать в векторном виде:

/ (У,е) = 0, (У1,У2,...У5). (15)

Для численного решения системы (15) использован регуляризованный аналог метода Ньютона [10] с параметром регуляризации а = 10-3:

У(п+1)(е) = У(п)(е) - гп[а/2(У(п)(е),е) + ?(У(п)(е), е)/' (У(п)(е), е)]-1 х

х 7(У(п)(е),е)/(У(п)(е),е), (16)

где п — номер итерации, тп (00 ^ тп ^ 1) —оптимальный шаг на п-й итерации, причём из соображений сходимости 00 ^ 0,1, / (У(п)(е),е) — матрица Якоби, /'(У(п)(е),е) — транспонированная матрица Якоби. Величина /2(У(п)(е),е) = 6п(т) представляет собой квадрат невязки. Оптимальный шаг

М0)+М1),

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

Распределение плотности считаем заданным, т.е. выбираем конкретные раЪс:

(а+Л

РаЪс = Д у Д Ра+Ъ,еУа+Ъ+СУС2 , (17)

\2/^ъ]'

где У1 = ^, У2 = е = О3. Мы взяли аксиально симметричный случай а,2 = аь Коэффициенты ра+ъ,с для Р = 6 приведены в табл. 1 и взяты из [7] для конкретной конфигурации, соответствующей модели быстро вращающейся нейтронной звезды с уравнением состояния Бете-Джонсона.

Таблица 1

Коэффициенты pa+b с для P = 6

а + b 0 0 0 0 2 2 2 4 4 6

c 0 2 4 6 0 2 4 0 2 0

Pa+b,c 1 -0, 729 0, 634 -0, 912 -0, 690 1, 294 -2, 765 0, 683 -2,845 -0, 996

Выбранная ньютоновская схема была реализована в рамках пакета MAPLE. Она оказалась достаточно эффективной, так как число итераций для достижения невязки 10-30 не превысило 10.

В результате получены численные значения Zijk и yi, y2. Они приведены в табл. 2. Для ра+ъ,с , приведённых в табл. 1, 5D будет порядка 10-6.

—ф

4

Таблица 2

Численные значения Zijk и y i, y 2

yi У2 Z002 Z020 Z200

1 0, 59237 — 1, 35538 ■ 10-6 6,46815 ■10-7 6,46815■10-7

2. Внутренний гравитационный потенциал

Фп(Р, L)

В работе [9] с использованием рядов Бурмана-Лагранжа доказана теорема, согласно которой внутренний гравитационный потенциал Ф-ш может быть представлен абсолютно и равномерно сходящимся степенным рядом по коэффициентам Zijk. При этом возникают некоторые ограничения на Zjk и коэффициенты при них, которые являются полиномами степени P + s(L — 2) + 2 координат Xk. Здесь s — номер члена ряда Бурмана-Лагранжа.

После перехода к обобщённым сферическим координатам R, 0, впервые предложенным в работе [8], со смещением центра координат в точку наблюдения Xk имеем:

xk = Xk + «kR,

1 (аЛ ■ а 1 (аЛ ■ а ■ 1 f aA а

«i = — — sin 0 cos a2 = — — sin 0 sin a3 = — — cos 0,

a \ai/ a \ai/ a \a3/

/ a2 a2 a2 ^ 1 (18) I an . 2 л 2 an . 2 л • 2 an 2 л

a = I sin 0 cos ^ +—2 sin 0 sin ^ +—° cos 0

ai ai a3

a0 = (a1,a2,a3)3, 0 ^ R ^ R.

В этом случае выражение для возмущённой эллипсоидальной поверхности (3) перепишется в следующем виде:

L

R2 + 2RT — Q = ^ Zijk(xi + aiR)¿(x2 + a2R)j(хз + a3R)k, (19)

где

33

T = 53 akXk, Q = 1 — xk.

k=1 k=1

Решив квадратное уравнение относительно R в левой части (19), получим уравнение для нахождения R:

D D (xi + aiR)i(x2 + a2 R)j (x3 + a3R)k

R = Ro — -R + T + U-=

ij'k (20)

= Ro + ^(R, ak, xk, Zi^k);

Ro = U — T, U = (T2 + Q) 2.

Элемент объёма интегрирования в новых координатах будет равен

dV' = a3 Í?2dÍR , d^ = sin 0 d0 d^.

a3

Дальнейшие вычисления будут основаны, как и в случае с Л, на использовании теоремы Лагранжа.

—е е

Согласно (8) в нашем случае £ = Я, /(£) = £н+2, а = Я0, Ф(а) = Ф(Я0,ак,Хкк). Тогда получим

вЬ

Ян+2 = Ян0+2 + (И + 2) ЕЕ

(-1)в

Zijk

в=1 ijk

г 3 к

п т I

Х

i — n j — m k — l n m l

в1

Х

2

Х

з

а1 а2 а

3 аяг1

г>Ь,+т+п+1+1 я0_

(Я0 + Т + и )в

(21)

х

Для сокращения записи здесь и далее используется введённый в [8] оператор

[ Щ п2 ... пт ]

суммирования [кг к2 ... кт\, действие которого заключается в суммировании по каждому нижнему индексу от нуля до значения стоящего над ним индекса с одновременным умножением на соответствующие биномиальные коэффициенты.

Для определения явного вида

Я

Н + т + п + 1 + 1

(Яо+Т+и )3

произведём замену Я0

и — Т + Уи, ёЯ0 = иёу и воспользуемся ранее приведённой леммой:

в1

ёЯ,

в1

г>Ь,+т+п+1+1 Я0_

(Я0 + Т + и )в

И + п + т + I + 1 Л

Яо=и-Т+уи Р

И + п + т + I + 1 Л

в-1

(—1)

2е ёув-1

(1 + у)Ь+П+ТО+1-Р+1

(—1)Р

иЬ,+п+т+1-2(в—1) — ртР х

х Фв(И + т + п + I — / + 1).

Учитывая полученный результат, Я примет вид:

3 — 1

а

0

в

^ вЬ (_1)в+р ЯН+2 = я£+2 + (И + 2) Е Е ЦЬ— Zijk(s)

в=1 ijk

в\

3 к И + п + т + I + 1 т I /

х x\-nx22-mxk-lаnаmаl3Uь+п+т+1-2(в-1)-ртрфв. (22)

х

п

Выражение для распределения плотности (6) в новой системе координат примет вид:

р

р(Р) = раЪс

аЪс

а Ь с

Л / д

5а-ёхЪ-/хС-9 Яа+/+9 а! а/

а1 а2 аз.

(23)

Ф-

Ф

-Ф

Ф

Выражение для ньютоновского гравитационного потенциала (1) в новых координатах с учётом (23) примет вид:

P

Ф(Р,£) = -Gag]T

abc

-Gag £

abc

abc

d f g

abc

d f g

xa-dxb-fxc-g 1 2

j R d+f+g+1djR| ^ff dß

3

a-d^b-f^c-g Rd+f+g+2

Ж /у» /у»

12

3

0

2 dR

ada2 af

a

2

dß-

P P sL

- Gag(h +

s=1 abc ijk

a b c i d f g n

j k h + m + n + l + 1 m l ß

_i)s+M

(-1)

a-d+i-n b-f +j-m c-g+k-l

2ss!

1

Ж.

2

3

d+n f+m g+l

x иh+n+m+l-2(s-1)-^T^ф a1 a2

j a2

a

3

dß.

(24)

h = d + / + g. Таким образом, (24) можно представить в виде:

Ф(Р,Ь) = -Gag(Fo + £ Fs)

s=1

(25)

где Fo и Fs соответственно первый и второй члены в (24). Fo — член аналитиче-

s=1

ского представления внутреннего потенциала, соответствующий невозмущённому эллипсоиду.

Для явного представления воспользуемся формулой бинома Ньютона для Е^+Ь:

Rh+2 _ (U - T)h+2 =

h + 2 ^ + 1 h + 2 - 2v Л v т С

ß

Р т С X

х х?(т-i)+h+2-2v-^x2(i-x)+A-pxfx+Pah+2-2v-Aa^af. (26)

Учитывая (26), Fo примет вид:

P

Fo = V Pa,b,c(-ir+T X

a,b,c

a b c h + 2 ^^ + 1 h + 2 - 2v Л v т С

d f g

ß

Р т С X

X Ж^+2-2^-A+a-d+2(r-£) ^ ^b-f+A-p+2(£-x) ^ xc-g+p+2x ^ QA A A (27)

A1 = d + h + 2 - 2v - Л; Ag = / + Л - p; A3 = g + p.

Здесь и далее Qf^f _ / afl af2af3 dar.

Воспользуемся, как и ранее, формулой бинома Ньютона. Получим

Uh+n+m+l-2(s-1)-^ _ у^ _

X

X

X

X

Л

v

X

X

Л

v

N-р 2

- (s - 1) V т £ N + 2 - 2s - 2v Л V т £ X Л Р

(-1)т х

X x

N-2(s-1)-2v-X+2(r-£) Х-р+2(£-х) р+2х N-2(s-1)-2v Х-р р

x

x

ал

а0

а

з-

(28)

Здесь и далее N = Л + I + т + п.

Учитывая (28), находим аналитическое представление для Рз:

(_1)м+т+з ^^ р = —^—

2ss!

abc ijk

a b c i j k N+1 d f g n m l и

(N- (s - 1) N + 2 - 2s - 2v Л v т £ v Л P т ? X

х ^sZijk (s) ■ x'a

a-d+j-n+N+2-2v-2 s-Х+2т -2£ b-f +j-m+X-p+2£-2x

x2 x

c-g+k-l+p+2x з

■ X

ЯЛгЛ^Лз

(29)

A1 = d + n + N - 2s + 2 - 2v - Л; A2 = f + m + Л - p; A3 = g + l + p;

s— 1

Ф1 = 1; Фэ>1 ^n(N +1 - л + 1 - 2r).

r=1

Из (27) и (29) следует представление Ф, содержащее все степени координат до Р включительно:

Р

Ф(Р,Ь) = -2nGpoa2J2 Ф abc xax2x3.

abc

Интегралы QA1A2A3 могут быть упрощены:

Q2

a,2b,2c

= 2п

I

2(a+b),2c = e

(2a - 1)!!(26 - 1)!!

2a+b(a + 6)!

-2c , (1 - x2)ax2c dx

" J (1 + (f - 1)x2)a+2c+1 ; -1 2

2(a+b),2c :

(30)

(31)

e =

озЛ

a1 /

В результате нами получено аналитическое представление Ф для реализации в системе символьных вычислений MAPLE.

Вычисление Фabc для значений Р = 4, 6, 8 и L = 2, 4, 6 нереально без использования метода символьных вычислений на компьютере. Как следствие нами была составлена программа в рамках того же пакета MAPLE для аналитического представления внутреннего потенциала в виде разложения по степеням координат.

Задача сильно упрощается, если задано распределение плотности, т. е. известны коэффициенты pabc. Мы взяли pabc, приведённые в табл. 3.

Подставив коэффициенты Zjk, мы получили потенциал как функцию координат, обозначив r2 = x2 + x2:

ФШ(Р = 6, L = 2) = -Gpoa2 (2, 7216 - 1, 6191r2 +

+ 0, 3540r4 - 0,1776r6 - 1, 068437x3 + 0, 20631x4 - 0, 079436x6 -

- 0, 3852x3r4 - 0,2954x3r2 + 0, 5279x^r2). (32)

Соответствующий график приводится на рис. 2, x3 = X3.

Поскольку программа оказалась большой и сложной, возникает необходимость её тестирования. Нами были выполнены два теста. Проводился предельный переход при e ^ 1 к сферически симметричному случаю, который легко рассчитывается по независимой программе. Кроме того, АФт(Р = 6, L = 2 = 4)= 4nGp(P = 6) и, действуя на Ф-ш оператором Лапласа, мы получаем распределение плотности, умноженное на 4nGp, коэффициенты в которой совпадают с заданными с

точностью 10-24.

1

2

х

х

Таблица 3

Коэффициенты раъс

а 0 0 0 0 0 0 0

Ь 0 0 0 0 2 2 2

с 0 2 4 6 0 2 4

РаЪс 1 -0, 729 0,634 -0, 912 -0, 690 1, 294 -2, 765

а 0 0 0 2 2 2 2

Ь 4 4 6 0 0 0 2

с 0 2 0 0 2 4 0

РаЪс 0, 683 -2,845 -0, 996 -0, 690 1, 294 -2, 765 1,367

а 2 2 4 4 4 6

Ь 2 4 0 0 2 0

с 2 0 0 2 0 0

РаЪс -5, 690 -2,989 0,683 -2, 845 -2, 989 -0, 9964

е—

е

—е

е

3. Потенциал во внешней точке (М,Р,Ь)

Наибольший интерес представляет потенциал во внешней точке Р, Ь),

когда г > «4, где г — расстояние до притягиваемой точки. Наиболее эффективный метод нахождения внешнего потенциала — метод разложения по мультипольным моментам. Внешний потенциал имеет вид:

Фои. = -С / , (33)

■> |г — г'|

где ёг' = ёж' ёу' ёг'.

Разложим функцию _1 ^ в ряд Тейлора:

|r - r^ - x')2 + (y - y')2 + (z - z')2

ж ( 1)а+в+7 / Яа+в+7 1 \

= Y" ^_ —__1 х'ау'вz'Y (34)

^ а'в'7! dyedzY ^x2 + y2 + z2 J У ' ()

Заменим в (34) ( —('axSf^ZY ^2+^2+^2) на опеРатоР Мав7> а ЧеРе3 1

' x2+y2+z

обозначим f 1 ==. Тогда (34) перепишется в виде:

л/x2+y2+z2

1 ^ ^ 1

—. = £ М«в7^Х'ау'вZ'Y. (35)

|r ' 1 ав1

Используя (35), получим выражение для внешнего потенциала:

^ i „

<£out = -G Е r J P(r')x'ay'ez'Y dx' dy' dz'. (36)

a0Y

Заменив в (36) интеграл на Dafí^, получаем:

^out = -G £ 1 DaeY. (37)

ав^

Учитывая, что xi = x, X2 = y, X3 = , Dae7 будет выглядеть следующим образом:

DaeY = aa+e+2aY+1 J p(r'Kxex3 dx1 dx2 dx3. (38)

Как и в случае внутреннего потенциала представим плотность конфигурации р в виде полинома (6) степени P.

Перейдём к сферическим координатам: xk = r«k, ai = sin в cos «2 = sin в sin «з = cos в, где 0 < r < Rmax, и, как в случае с вычислением Л, воспользуемся теоремой Лагранжа. После этого выражение для возмущённой эллипсоидальной поверхности (3) в данном случае будет иметь вид:

ж sL (_i)s

Rtax = 1 + h££ ^Г"Zijk+ i + j + k - 1)aia22«k.

s=1 ijk

Положив Qabc = / a-, «2, «3 d^ получаем

1

1

N P

Фоп^,Р,Ь) = (-poG) Е Е PabcMaen X

а@ Y abc

1 го ijk ( — i)s

X ' ho + 3 <Qa+a,e+b,Y+c + E E — ss| Zijk(s)(^s(h + 1)Qa+a+i,e+b+j,T+c+fc

(39)

s- 1

r=1

где h0 = a + в + Y + a + a + c, h = h0 + i + j + k, а Фs(h + 1) = П (h + 2 — 2r); в

случае s =1 Ф., = 1.

Нетрудно доказать, что:

Q abc =

(a — 1)11(6 — 1)!!(c — 1)!! (a + b + c + 1)!! '

(40)

Учитывая (40), выражение для внешнего потенциала можно представить в виде разложения по мультипольным моментам следующим образом:

NP

ФоиЬ^,Р,Ь) = — 4npoG Е ^РаЬеМавл X

ав y abc

1 (a + a — 1)!!(в + b — 1)!!(y + c — 1)!!

ho + 3

(ho + 1)!!

+

го ijk (_i)s

^^e Zijk(s)фs(h + 1)x

(41)

^ 2ss!

s=1 sL

(a + a + i — 1)!!(в + b + j — 1)!!(y + c + k — 1)!! X (h + 1)!!

где N — порядок мультиполя.

На основе (41) была составлена программа в системе символьных вычислений MAPLE. Как сказано выше, распределение плотности считаем заданным, т.е. известны коэффициенты pa+b,c, зная которые и используя (10), находим pabc.

Коэффициенты pa+bc приведены в табл. 4 на основании [7], в которой нами рассмотрены три уравнения состояния ядерной материи: Бете-Джонсона (BJ), Оппенгеймера-Волкова (OV) и Рейда (R).

Коэффициенты ра+Ъс

Таблица 4

OV OV BJ BJ R R

a + b c pa+b,c pa+b,c pa+b,c pa+b,c pa+b,c pa+b,c

0 0 1 1 1 1 1 1

0 2 —0,46117 —0,46375 —0,72181 —0, 72895 —0,68902 —0,69547

0 4 0, 51814 0, 51818 0, 63580 0,63397 0, 62306 0, 62165

0 6 — 1, 05713 — 1,05863 —0, 91425 —0, 91169 —0,93429 —0,93256

2 0 —0,45871 —0,44935 —0, 71544 —0, 69125 —0,68320 —0,66110

2 2 1, 03862 1,04169 1,27984 1,29401 1,25331 1,26543

2 4 —3,17536 —3,17357 —2, 75605 —2, 76435 —2,81462 —2, 82108

4 0 0,520818 0, 53272 0, 64485 0,68239 0, 63099 0, 66416

4 2 —3,18038 —3,19872 —2, 77129 —2, 84334 —2,82819 —2, 89133

6 0 —1,06216 — 1,08540 —0,92951 —0,99437 —0,94788 —1,00615

—0

4

Используя данные табл. 4 и выражение для внешнего потенциала (41), в рамках того же пакета MAPLE получено выражение для внешнего потенциала сложной эллипсоидальной конфигурации в квадрупольном приближении для N = 2 в сферических координатах:

^out(N = 2, Р = 6, L = 2) = -^ (l + - 3cos2 ^ + ..., (42)

где m — масса гравитирующей конфигурации.

В работе [7] показано, что а2 = 2nGp2y(e), где Ро —плотность в центре конфигурации, Ро —давление в центре конфигурации, y(e) — рассчитанные функции параметра сплюснутости e и учитывающие быстроту вращения конфигурации. При изменении угловой скорости вращения функции Ро и ро также будут зависеть от параметра сплюснутости e

Ро = PooY (e), Ро = Роох(е), (43)

С учётом (43) имеем:

^ = D(e) А ' ^Ш) = ^' (44)

где

n (e) n(e) Y(e) K(р ) роо(Роо)

Dl(e) = n(e)хеш, K(роо) = ^оо^рю.

После проведённых преобразований выражение для внешнего потенциала возмущённой эллипсоидальной конфигурации (42) примет вид:

^out(N = 2, P = 6, L = 2) = - ^ (l + Di(e)K (роо)(1 - ^ + .... (45)

График функции D (e) представлен на рис. 2, а значения K(роо) приводятся в табл. 5.

Рис. 2. Зависимость коэффициентов от параметра е: кривая 1 соответствует уравнению

состояния ОУ, 2 — БЛ, 3 — И.

Как видно из рис. 2, при е =1 обращается в ноль, а с уменьшением параметра е увеличивается, т. е. коэффициент разложения в квадрупольном приближении монотонно зависит от параметра сплюснутости конфигурации. Таким образом, чем больше сплюснутость конфигурации, тем существенней вклад квадрупольного члена.

Вестник РУДН, Серия Математика. Информатика. Физика. № 1. 2008. с. 28-42 41

Таблица 5

Значения K(poo)

Р00, г/см3 K (OV), см2 K (BJ), см2 K(R), см2

4 ■ 1014 6 ■ 1014 8 ■ 1014 1.418 ■ 1011 1.284 ■ 1011 1.160 ■ 1011 1.911 ■ 1011 2.379 ■ 1011 2.636 ■ 1011 1.200 ■ 1011 1.454 ■ 1011 1.670 ■ 1011

Заключение

В результате нами получено аналитическое представление ньютоновского гравитационного потенциала возмущённых эллипсоидальных конфигураций на внутреннею точку при фиксированных значениях Р и L в виде абсолютно сходящихся рядов по параметру возмущения (25). В рамках пакета MAPLE была составлена программа для представления потенциала Фт(Р, L) в виде функции координат (32) для произвольных значений Р, L и реализована для случая Р = 6, L = 2. Потенциал возмущённой эллипсоидальной конфигурации на внешнею Фоиь(^, Р, L) точку мы представили в виде разложения по мультипольным моментам и получили его конкретное аналитическое представление (41) в квадру-польном приближении. С помощью составленной программы в системе символьной математики MAPLE и численных вычислений регуляризованным аналогом метода Ньютона были получены выражения для потенциала Ф1П(Р = 6, L = 2) и = 2, Р = 6, L = 2) в виде функций координат (45). Результаты численных расчётов представлены в виде таблиц и графиков.

Литература

1. Сретенский Л. Н. Теория ньютоновского потенциала. — М.: Гостехиздат, 1946. — С. 316.

2. Тассуль Ж. Л. Теория вращающихся звезд. — М.: Мир, 1982.

3. Муратов Р. З. Потенциалы эллипсоида. — М.: Атомиздат, 1976.

4. Антонов В. А., Тимошкина В. И., Холшевников К. В. Введение в теорию ньютоновского потенциала. — М.: Наука, 1988.

5. Masjukov V. V., Tsvetkov V. P. Nonlinear Effect in Theory of Equilibrium Gravitating, Rapidly Rotating, Magnitized Barotropic Configurations and the Gravitational Radiation from Pulsars // Astron. and Astrophys. Transactions. — Vol. 4. — 1993. — Pp. 41-42.

6. Лаврентьев М. А., Шабат Б. В. Метод теории функции комплексного переменного. — М.: Наука, 1987. — 688 с.

7. Математическая модель гравитирующей быстровращающейся сыерхплотной конфигурации с реалистическими уравнениями состояния / Е. В. Беспалько, С. А. Михеев, И. В. Пузынин, В. П. Цветков // Препринт Р11-2005-35. — Дубна: ОИЯИ, 2005.

8. Цирулев А. Н., Цветков В. П. Вращающиеся постньютоновские конфигурации однородной намагниченной жидкости, близкие к эллипсоидам, I,II. // Астроном. — Т. 59. — 1982. — С. 476-482, 666-675.

9. Цветков В. П., Масюков В. В. Метод рядов Бурмана-Лагранжа в задаче об аналитическом представлении ньютоновского потенциала возмущенных эллипсоидальных конфигураций // ДАН СССР. — Т. 313, № 5. — 1990. — С. 1099-1102.

10. Ермаков В. В., Калиткин Н. Н. Оптимальный шаг и регуляризация метода

Ньютона // Жур. вычисл. матем. и мат. физ. — Т. 21, № 2. — 1981. — С. 419-

497.

—0 4

ф

ф

UDC 519.6,517.9

Calculation of Newton Potential of the Gravitating Configurations with the Surface Close to the Spheroid by Symbolic and Numerical Methods

E.V. Bespalko, C. A. Mikheev, V. P. Tsvetkov, A.N. Tsyrulev,

I. V. Puzynin

An exact analytical representation of the heterogeneous perturbed ellipsoidal configuration of the Newton gravitational potential to the inner point <Pin is obtained. An expression for the potential <Pin is obtained in the form of a polynomial of the coordinates by the program prepared in the system of symbolic mathematics MAPLE, and by numerical calculations based on the regularized analogue of the Newton method. A formula for the potential to the outer point <Poat is found by the quadrupole approximation.

Department of General Mathematics and Mathematical Physics Tver State University 33, Zhelyabova st, Tver, 170100, Russia

e— e

—e e