Применение многомерных сферических координат для численного интегрирования и некоторых других задач Текст научной статьи по специальности «Математика»

УДК 519.95

ПРИМЕНЕНИЕ МНОГОМЕРНЫХ СФЕРИЧЕСКИХ координат ДЛЯ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ И НЕКОТОРЫХ ДРУГИХ зддач

С. Д. Субочев,

канд. техн. наук

Санкт-Петербургский государственный университет аэрокосмического приборостроения

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

The numerical method of integration of functions in multidimensional spherical coordinates, which is convenient for the practical applications, has been worked out. Application of method of numerical integration for calculation of probabilities and moments on multi-dimensional distributions is described. The possible use of multidimensional spherical coordinates for extreme tasks on surfaces and some other tasks are shown.

Многомерные сферические координаты

Общеизвестно соотношение между трехмерными декартовыми прямоугольными координатами х0, х-|, х2 и сферическими координатами г, ос0, :

х0 = г/0, х1 = rlb х2 = rl2

/0 = cos cos а0, А| = cos а-| sin а0, /2 = sin а0

(1)

где г-длина радиус-вектора, /0, /-, и /2- его направляющие косинусы; а0и оц - углы долготы и широты, соответственно.

Обобщая по индукции проектирование в соответствующие подпространства и на координатные оси (этот процесс здесь не приводится в силу сложности его описания и физического представления), приведем окончательные связи между углами и декартовыми координатами для л-мерных пространств:

х,= rlh / = 0, 1, 2,л-1, (2)

где /, - направляющие косинусы.

Направляющие косинусы равны соответственно

/0 = cos an_2cos an_3...cos c^cos c^cos a0 /-) = cos an_2cos an_3...cos a2cos a-|Sin a0 12 = cos an_2cos an_3.. .cos a2 sin 0|

/3 = cos an_2cos an_3...sin a2 [■, (3)

In-2 = cos an_2sinan_3

/n_1 = sin a

n-2

где некоторые углы ау - обобщенные физические аналоги углов долготы и широты.

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

е = [/0,/1,/2> ...,/л_1]т;Н = [х0,х1>х2, ...,хп_1]г=ге, (4) несложно показать, что нормы векторов равны

IIе! = '\//02 + /12 + --- + /л-12 = 1;

(5)

И = >/х02 + х12 + ... + хл„12 =г,

т. е. концы векторов е и Р лежат на сферах радиусов / и г, соответственно. Следовательно, параметр г > О является текущей длиной радиуса-вектора Я, а параметры г, «0, ос-1,..., осп_2в целом являются л-мерны-ми сферическими координатами. Любая линейная координата х,- может принимать два значения (положительное либо отрицательное). Несложно убедиться, что при изменении углов в пределах

О < а0 < +2л, -к/2 < < +к/2, ..., -к/2 < ап_2 < +я/2, (6)

координаты х0, хь х2, ..., хп_1, рассчитанные по формулам (3), принимают 2п всех возможных сочетаний положительных и отрицательных значений, при этом единичный радиус-вектор е полностью заметает конус всех возможных направлений в л-мерном пространстве, попадая на одно и то же направление лишь единожды.

Бесконечно малый элемент объема в п-мерных сферических координатах

Этот элемент равен произведению определителя матрицы О - якобиану преобразования координат (3) на соответствующие дифференциалы:

dV = Ое\[О№гіа0даі..^ап_2. (7)

Столбцами матрицы в являются векторы - производные радиус-вектора В по аргументам г, (ц, щ,

а,,. 2.

G =

т ди_ дв_

дг ’ да0' да-\

9R

да

п- 2

(8)

Можно показать, что все вектор-столбцы в матрице й попарно ортогональны, при этом детерминант матрицы равен произведению норм этих векторов:

det[G]

ЭП 3R 3R

II Эг! • Эа0 Эа-, Эап_2

(9)

Расписывая эти нормы по индукции и подставляя выражение (9) в (7), получаем

сП/п = со5п_2ап_2... х х ... соэ1со5°а0/'л“1 drda0da.1 ... с/ап_2. (10)

Пусть область интегрирования функции Г (г, а0, а1» •■•> ап-2) в сферических координатах ограничена гладкой замкнутой гиперповерхностью Б, и при этом радиус-вектор, исходящий в любом направлении, пересекает эту поверхность изнутри наружу один раз, т. е. угол между нормалью к внешней стороне поверхности и радиус-вектором всегда меньше я/2. В этом случае с учетом полных пределов изменения углов (6) представим интеграл от этой функции в виде

л/2 /= J COS

-п/2

к/2

п~2ап_2 f cos'7_3an_3 -л/2

2л

я/2 2л Re

X... J cos1 ос 1 Jcos°tto J rn~^x

-л/2 о о

X F{r, (Xq, <*!,..., an_3, an_2)drda0daA ... dan_3dan_2l (11)

где Re = Rend (a0, a,..an_3, an_2) - длина радиус-век-

тора, исходящего от начала координат в направлении, заданном углами а0, а-,, ..., ап_3, ап_2, до конечной точки пересечения с поверхностью.

Чтобы, например, вычислить объем л-мерного шара радиуса г0, надо принять в формуле (11) длину радиус-вектора и удельную плотность постоянными, соответственно равными

Fte = r0 й F(r> <*о> «1, Оп-З. ап-2) = 1 •

После необходимых преобразований получим объем л-мерного шара в виде

2(2я)м/2г0п[1 • 3 • 5 •... • (л - 2) • л]-1, если л-нечетное,

(2я)п/2г0п[2 • 4 • 6 •... • (л - 2) • п]~\ если л-четное.

1/С =

(12)

Результаты выражение (12) совпадают с известными формулами для объема л-мерного шара [1,2], полученными «формальным» и косвенным способами, соответственно, что подтверждает достоверность вышеприведенных выводов и формул, полученных автором путем использования многомерных сферических координат.

Необходимые свойства многократных сферических интегралов

Возьмем в выражении (11) (у + 1) раз интегралы по аргументам г, а0, ос,, ..., ау_1 и обозначим интеграл, зависящий от параметров ау, ау+1, ..., ап_2, как

п/2

// = //(ау, аг-Ъ •••■ ап-2) = | СОБНам ... X

-л/2

л/2

2л

R е

х... [ cos1a1 [cos°a0 f rn 1х

-л/2 0 0

(13)

х Р(г, а0, оц,..., сх/_1/сху, сху+1,..., осп_2) с/гс(а0с/а1... с/ау_1,

где Яе = Иепс, (а0, аь ..., а^/ау, ау+1.....ап_2).

Отметим, что кратность интеграла (13) по аргументам - углам ос0, аь ..., ау_-| равна;. Запишем связь между интегралами, имеющими кратности по углам (У+1)иу:

/у+1(ау+11ау+2,..., ап_2) =

л/2

= | со5;а///(ау/а/+1, ау+2,..., ап_2)с/ау. (14)

-п/2

Обозначим подынтегральную функцию по^глу ау

как

^^у(ау/ау+1,..., ап_2) = соб(ху/у(ау/ау+1,..., ап_2),(15) а численные значения производных на границах интервалов интегрирования как

f (2Аг-1) _ -------

у Эсх

2А:

к= 1,2,3...

ау = ± я/2

і (2т). 7

д2т1

(16)

да

2т , т = 1, 2, 3, ...

| ау- = ± я/2

Отметим необходимые далее свойства единичного вектора направления е и составляющих его направляющих косинусов в окрестности точек ау = ± я/2 ± Дау:

е| «0, а-|, а2, ..., ау = ±(я/2-Дау); е I а0-я, -а-|, -а2, ..., ау = ±(я/2+Дау). (17)

Нетрудно убедиться, что направляющие косинусы (3) будут равны для отсчетов углов, записанных слева и справа в выражении (17), и что такое равенство соблюдается по углам с любыми номерами (у =

= 1,2, л-2), взятыми в качестве полюсов осу = ±ге/2,

если проинвертировать значения углов с меньшими номерами а, = -ос, (0 < / < у), а угол с нулевым номером изменить как ос0 = о.0-к. Следовательно, если единичный вектор направления пересекает полюс по углу ау = ±к/2, перемещаясь на величину ±Асху, то для этого вектора единственным образом находится равный вектор направления, при этом диапазоны интегрирования по углам с меньшими номерами приводятся к прежним диапазонам

-л/2 <-а, < тс/2, -л<а0-я< к. (18)

Кроме того, частичный элемент объема в подпространстве размерности 0+1)

dVM = cos'1 ocy_i... х х... cos1 gctCOS0ос0гп-1с/гсУос0с/ос1..,dan_2 (19) при такой замене углов не меняется.

Из перечисленных свойств следуют свойства четности интегралов /; (ос,) относительно полюсов осРу= ±я/2:

/у (a.pj+До,) = /у (аРу - Дау). (20)

Из свойств четности интегралов /у (о,-) несложно доказать полезные (и необходимые далее) свойства равенства нулю всех производных нечетного порядка от подынтегральных функций по углам с четными номерами

f/2*-1> = 0, ^ = ±71/2,7=2,4, 6...; /с=1, 2, 3... (21) и ряда нечетных производных по углам с нечетными номерами

Эо.3 -°' s Эа5 - dal

dfj_ ^7 : = 0,...

Эау Эау

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

Численный метод интегрирования функций внутри замкнутой гладкой области в многомерных сферических координатах

0. Зафиксируем углы а0, аъ ..., ап_2как параметры и допустим, что найдена численно или аналитически длина радиус-вектора от начала координат до точки пересечения с поверхностью, ограничивающей область интегрирования Яепс) = Яепс)(а0, щ, ..., ап_2).

Обозначим интегрируемую функцию по радиусу г, как по аргументу. В ней углы а0, аъ ап_2 являются параметрами:

/го(/') = /го(''М)> «1. •••> осп_2) = /л"1Р(г/а0, си.ал_2). (23)

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

~)2/<-1р л2/с-1гг

Ро1*'"(0)= “ЭГГ и ’№„<,)= г• (24)

|г=0 I г = А,,„0

Представим внутренний интеграл по радиусу, применяя формулу Эйлера-Маклорена [3], в виде приближенных конечных сумм из N ненулевых отсчетов подынтегральной функции и т уточняющих поправок, учитывающих разности нечетных производных на концах интервала интегрирования:

/о — /о(оо. «1.ап_2) =

я_е Л/-1 1

= | Г0 (г) dr~ ]Г г0 (/Д г)Дг+ - Я0(Яепс/)Дг +

0 '=1

т о

+ (2Л)1 Д'21‘[^о|21‘'1)(0) - Ро12*"'1 (Я»м)]. (25)

где Дг = ЯепсУ/Л/ - шаг суммирования отсчетов функции.

При этом на первом шаге (и далее) положим, что используется формула Эйлера-Маклорена.

1. Последующий интеграл по углу а0 (при фиксированных углах (а,......а„_3, а„_2) представим в виде

конечной приближенной суммы из Л/0 отсчетов подынтегральной функции:

Мон, <Х2....а„_2) =

2л

= | соз0а0/0(а0/а1, а2,..., ап_2)с/а0 * о

N0-1

“ £ *о (/Аос0/а-|, а2,.... ап_2)Да0, (26)

/=0

где Да0 = 2тс/Л/0.

Уточняющих поправок не требуется, так как в начале и конце интервала интегрирования при а0 = 0 и а0 = 2л; радиус-вектор имеет равные направляющие косинусы (3), попадая на одно и то же пространственное направление. Следовательно, подынтегральная функция является периодической по углу а0 (период равен 2тт) и поэтому разности всех производных и соответствующих им уточняющих поправок равны нулю.

2. Следующий интеграл по углу ат представим в виде конечной суммы ИЗ (Л/т - 1) отсчетов и, например, пяти уточняющих поправок:

л/2 МН п

/2(а2, а3,ап_2) = | » £ ^1 (- 2 + /Дсц) +

-71/2 '=1

^ В

+ 2(адДа12Л[^(2/<-1>(-7г/2)-^(2/с-1)(+7с/2)]( (27) где Дат = к/Л/,;

Ма-,) и /:1<2'(~1) (±я/2) - подынтегральная функция и ее производные, определенные в соответствии с (23) и (24).

Здесь и далее ограничимся практическим и очевидным случаем, когда заведомо нулевые отсчеты отбрасываются: //(«, = ± к/2) = 0, так как cos' (±л/2) = 0.

Производные нечетных порядков от подынтегральной функции (а-,) - cos1 a-i/^a-i), применяя формулы многократного дифференцирования произведений, удобно рассчитывать через производные четных порядков (на один порядок ниже) функции l-iia-i). Приведем эту связь между численными значениями производных до 9-го и 8-го порядков при а, = ± к/2:

f1ti)=-/1l /р=-3/1(2) + /1> f1<5>=-5/1(4) + 10/1(2,-/1l f}7) = -7//6) + 35/|4) - 211\2) + ,

= —9/1<S) + 84/,(6) -126/1{4) + 36/1(2) - lv

(28)

На основе формул (28) можно рассчитать для приближенной суммы (27) пять уточняющих поправок на ошибки до 10-го порядка малости, включительно.

3. Из свойства (21) следует, что все уточняющие поправки в применяемой формуле Эйлера-Маклорена равны нулю и следующий интеграл по углу а2 заменяется конечной суммой из (N2-1) отсчетов:

/з = /з(«з.«п-з. <V2) =

л/2

= J cos2a2/2(a2/a3, ..., a„_3, an_2)da2 *

-я/2

N2-1 л _

« ^ cos2 (- 2 +/Да2)/2(- 2 +/Да2)Да2, (29)

/=1

где Дос2 = я/Л/2.

4. Интеграл по углу а3 вычислим как конечную сумму, взяв пять уточняющих поправок на ошибки до 12-го порядка малости, включительно:

/4(а4, а5, ..., ап_2) =

к/2 л/3-1 я

= j f3 (a3)da3 - X f3 (- 2 +/Ааз) +

-я/2 '=1

+ Х^Даз2*^2*-1^-^)-^-1^^)], (30)

где Да3 = тс/ Л/3.

При этом, как и на втором шаге, можно воспользоваться аналогичными связями между производными при а3 = ± л/2:

f3(1) =0, f3(3)=-6/3, f3(5)--60/3(2) + 60/3, f3(7) = -210/3(4) +1260/3(2) - 546/3,

/з(9) = -504/3(6) + 7560/3(4) -19656/3(2) + 4920/3, ^(31) f3(11) = -990/3(8) + 27720/3(6) - 80180/3(4) + +27060/3(2)-44286/3.

Необходимые производные вида /т*2"7’ по углу а-, на предшествующем втором шаге, 13{2т) по углу а3 на данном четвертом шаге и т. д. (по углам а5, а7, ... на последующих четных шагах) удобно вычислять численно по четырем узлам, используя четность функций - интегралов /1, /3, /5, ... относительно полюсов:

/у(2) - [8064/у (1) -1008/у (2) +128/у(3)- 9/у(4)]/(25205а2),

/ ;(4) = [-19521: (1) + 676/,- (2) - 96/. (3) + 7/. (4)]/(1205а?),

\ (32)

/у(6> ~ [116/у (1) - 52/у (2) +12/у (3) - /у (4)] /(25а®),

/у(8) - [-112/у(1) + 56/у (2) -16/;(3) + 2/у(4)]/5ау,

где /у(/) = /у (±л/2 + /5ау) - первый отсчет функции, четной относительно ау = ±я/2; у = 1, 3, 5, ...;/ = 1, 2, 3, 4; 5ау-выбранный шаг численного дифференцирования.

5. Следующий интеграл по углу а4 представим суммой отсчетов без уточняющих поправок, аналогично, как и на третьем шаге по углу а4, и т. д.

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

Итак, в предлагаемом алгоритме интегралы по радиусу - г и углам с нечетными индексами - аь а3, а5, ... рассчитываются как суммы с уточняющими поправками, а интегралы по углам с четными индексами а0, а2, а4, ... - без этих поправок. При этом циклы вычислений по углам ау организуются как внутренние циклы по отношению к циклам по углам ау+1, самый внутренний цикл вычислений - по радиусу. Необходимый переход от многомерных прямоугольных координат к сферическим осуществляется по формулам (2) и (3).

Предлагаемый метод численного интегрирования в многомерных сферических координатах удобен тем, что во многих практических случаях гладкая замкнутая область интегрирования в угловых пределах является параллелепипедом и при разбиении не имеет элементарных граничных участков, частично не попадающих в нее. В сравнении с методом Монте-Карло [4], который дает лишь вероятностные оценки искомых интегралов, в предлагаемом методе не требуется громоздкого статистического моделирования распределений случайных величин. Обнаруженные свойства четности подынтегральных функций относительно полюсов по углам с четными номерами позволяют рассчитывать интегралы в виде конечных сумм [рассмотрение погрешностей из-за несходимости рядов вида (25) выходит за рамки статьи]. Результаты изложены по убывающей общности, также возможны изменения и развития алгоритма в интересах пользователя. Так, при численном интегрировании по нечетным углам можно использовать не только формулу Эйлера-Маклорена, но и другие известные способы уточнения.

Примеры численного интегрирования в многомерных сферических координатах

Приведем два примера решения задач, которые иллюстрируют использование метода и подтверждают его достоверность. Каждая задача решается дважды на основе известного метода и предлагаемого. Известный метод является точным и эталонным, предлагаемый - приближенным и тестируемым. Результаты решения предлагаемым методом сходятся к точным результатам.

Пример 1. Вычислить полярный момент второго порядка для четырехмерного эллипсоида с единичной равномерной плотностью. Уравнение эллипсоида:

(х0 -0.9)2 , (Хї-1,1)2 , (х2-1,5)2 , (х3 -1,7)2 ,

„9 + . -О + . _0 ” » ■

10"

12"

16"

(33)

Центр этого эллипсоида смещен в точку с координатами

Кс= [хсо> хсь хс2> хсз]Г= [0,9; 1,1; 1,5; 1,7]г. (34) Полярный момент второго порядка рассчитывается относительно полюса - начала координат и равен интегралу от квадрата радиус-вектора по объему, заключенному внутри поверхности эллипсоида:

МР= \p(R)R2dV

(35)

где р(Я) - удельная плотность, как функция радиус-вектора (далее для упрощения примера примем р (Я) = const = 1).

Можно показать, что полярный момент второго порядка смещенного эллипсоида равен

МрА - MpQA + VA ||яс||2

(36)

где МРСА - полярный момент центрированного эллипсоида; 1/л-его объем; ||ЯС||2 =(хс0)2 + (хс1)2 + (хс2)2+ + (хс3)2- квадрат нормы вектора Р?с

Данный эллипсоид (33) получается из шара с радиусом г0 = 2 при коэффициентах растяжения по соответствующим осям /с0 = 3, /с1 = 5, к2 - 6 и к3 = 8 раз. Тогда, используя свойства подобия и формулу (12) при п = 4, получим

Мрл'ст = ку(МК/А)

{к02 + к,2 + к22 + к32) + к^фс12, (37)

где ку = к0к^к2к3 -- коэффициент увеличения объема;

\/5 = 2л2г04/4 и МРЗ = 2тс2г06/6 (38)

- объем и полярный момент соответствующего шара. Подставляя в формулы (36) и (37) необходимые величины, рассчитаем численное значение полярного момента эллипсоида, которое и примем за точное (корректное):

МРА{С0Й) = 5 485 541,917... . (39)

Полярный момент в сферических координатах несложно привести к интегралу:

л/2 л/2

МрА= J COS2 «2 J COSO^X

-п/2 -п/2

2п

J [Rend(aо- av a2)f da0da^da2, (40)

о

где Яепс, (ос0, а1( а2) - конечная длина радиус-векто-ра, исходящего из начала координат в направлении, заданном углами а0, аъ а2, до точки его пересечения с поверхностью эллипсоида. Чтобы найти это значение, необходимо уравнения связи (3) при п = 4 подставить в уравнение эллипсоида (13) и, решая получившееся квадратное уравнение относительно г, выбрать положительное значение корня (это тривиальное решение здесь не приводится). 0-й шаг алгоритма - вычисление [Rend (ос0, ос-,, ос2)]6 в этом примере производится аналитически без шагов по Дг. Производя оставшиеся шаги по углам ос0, сц, а2 согласно описанному алгоритму, вычислим приближенное (аппроксимированное) значение полярного момента эллипсоида МРА(ЛРРЯ).

Пример 2. Вычислить объем шестимерного эллипсоида

*о

*1

х2 (х3 -1,2)"

+ —+ —---

3,9 3,9 3,9 5,1"

(х4 -1,5)2 (х5 -2,1)"

6,9"

5,7"

= 1.

(41)

Этот эллипсоид получается из шара радиусом г0 = 3 при его растяжении по соответствующим осям в к0= /ст = к2 = 1,3 и к3 - 1,7, кл = 2,3, к5 = 1,9 раз. Тогда объем данного шестимерного эллипсоида больше объема соответствующего шестимерного шара \/5 в к0к^...к5 раз. Подставляя необходимые величины п = 6 и г0 = 3 в формулу (12), рассчитываем численное значение объема эллипсоида, которое будем принимать за точное:

^(соя) = к0к,к2к3к4к5У5= 61 487,425.... (42)

Во втором примере уравнение эллипсоида (41) специально задано так, что длина радиус-вектора ЯепсУ не зависит от углов а0 и «!, а зависит только от углов а2, сс3 и а4 (что несложно показать). Кроме того, в обоих примерах центры эллипсоидов специально смещены относительно начала координат, чтобы уничтожить осевую симметрию, которая препятствовала бы проверке свойств четности интегралов относительно полюсов. Записывая соответствующую формулу для объема и взяв первые три внутренних интеграла по г, а0 и аналитически, (т. е. выполняя 0-й, 1-й и 2-й шаги алгоритма), получаем

2 к

3

л/2 п/2

I cos4 ос4 | cos3 а3 х

-л/2 -тс/2

тс/2

X | соэ2а2 [Яепс1(ос2, а3, а4)]6с1а2с1а3с1а4. (43)

-к/2

(Аналитическое интегрирование введено для того, чтобы снизить количество необходимых вычислений на три порядка и соответственно иметь возможность по быстродействию рассчитать этот пример на ПК). Далее, производя 3-й. 4-й и 5-й шаги алгоритма - численное интегрирование по углам а2, а3 и а4] вычисляем приближенное (аппроксимированное) значение объема эллипсоида ]/А{АРРЯК Результаты решения примеров сведены в таблице.

Искомые величины - полярный момент четырехмерного и объем шестимерного эллипсоидов - вычислялись с точностью до ошибок от 2-го до 12-го и от 4-го до 14-го порядков малости, соответственно. Фактические относительные ошибки вычислений при порядках малости 2к соответственно равны

МРА^-МРА{СОП\

$Яе/-р(М)- м (СОЯ)

Мрл

^(ДРРЯ)_^(СОЯ)

5яе/-р('/)= 1/ (СОЯ)

м

Минимальные относительные ошибки достигают десятичных порядков 10-14— 10-15.

Величины, которыми аппроксимировались фактические относительные ошибки вычислений 2/с-го порядков малости (и которые можно принять за относительные теоретические ошибки), составляют

5«е/-г= В2кЬа?к[^2к-'\-к/2) - у2к~'\+к/2)]/(2к)\, (44)

где у = 3 иу = 5 при вычислении полярного момента и объема, соответственно. Также в таблице приводится относительное расхождение между истинными относительными ошибками и соответствующими аппроксимирующими величинами:

(8яе/-г~ $яе/-г)( §яе/-я) 1100 %.

Эти расхождения не превышают -6,47 % для первого и -5,18 % для второго примеров, т. е. фактические относительные ошибки с порядком мало-

сти 2к можно аппроксимировать первым отброшенным членом вида (44) (если, разумеется, предшествующие интегралы с меньшей кратностью вычислены более точно).

Шаг численного дифференцирования в обоих примерах выбран равным 80ц = 8«3 = 1 /256. Количества циклов по углам для первого и второго примеров соответственно равны

Ы0 = 128, Л/т = Л/2 = 64; Л/2 = Л/3 = Л/4 = 64. Относительные ошибки с повышением порядка малости уточняющих поправок уменьшаются на два порядка.

О применении численного интегрирования в многомерных сферических координатах для некоторых задач статистической радиотехники

Разработанный автором метод численного интегрирования был практически использован и апробирован в вышеназванных применениях. Приведем две задачи в общем виде, для решения которых предлагаемый метод оказался удобным и целесообразным.

Задача!. Найти вероятность Р попадания многомерного вектора ошибки И = [х0, хь ..., хп_1]тв определении параметров сигнала в некоторую область I/, заключенную внутри поверхности Б = 5(х0, х1( ..., хп_-|) = С, при функции плотности распределения вероятностей ошибок /"= Цх0, хь ..., хп_-,), которая в общем случае не интегрируется аналитически:

Р = Р(йс\/) = Р[3(Я)<С] =

= _|7(х0,х 1, ...,хл_1)с(х0^х1...с/хп_1. (45)

V

В статистической радиотехнике, теории автоматического управления,навигации и других областях часто рассчитывают вероятность того, что модуль многомерного вектора ошибок не превысит заданный радиус. При этом поверхность Б, ограничивающая область интегрирования, превращается в ^-мерную сферу и очевидно удобство перехода к многомерным сферическим координатам, заключающееся в том, что пределы интегрирования по радиусу и по

Порядок малости ошибки 2к Ошибки вычисления МРА, п = 4 Ошибки вычисления \/д, п = 6

§я е,-т(М) 5бЯе; (М), % §яе/-г('/) «я 55яе,(Ю,%

2 -6,12 10-4 -6,11-10-4 -0,09 - - -

4 -5,36-10-7 -5,33-10"7 -0,41 -7,48-10"7 -7,43-10'7 -0,70

6 -2,19 10-9 -2,17-10"9 -0,87 -5,25-10-9 -5,19 10-9 -1,13

8 — 1,90 10-11 -1.87-10-11 -1,44 —5.91-10-11 —5,81 -10-11 -1,70

10 -2,74-10”13 -2,56-Ю”13 -6,47 -1,01 Ю-13 -9,80-10-13 -2,46

12 — 1,84-10-14 - - -2,47-10'14 -2,35-10-14 -5,18

14 - - - -1,30-1 о-15 - _

углам становятся независимыми. Если ограничивающая поверхность не является сферой, то следует от уравнения поверхности

БСхо, х^,...,х„_1) = С перейти к уравнению ее радиус-вектора как функции углов:

= ^е/)с/(^0’ 1 2)’

Соответствующий переход делается и в подынтегральной функции:

ЯСХО.Х!, ...,х„) -»Р[Я(ао, аъ ..., €*„_■,)].

Задача 2. Найти моменты для распределения на области определения I/:

М(р0,рь ...!рл/_1) =

= |х0Рох1Р1...хл/_1Рл/"1/:(Хо.х1, ...)хл/_1)с(х0сУх1...с(хл/_1> (46)

I/

гдер0, Рт...Рл/-1 = о, 1, 2, 3, ... .

Эта задача встречается в статистике и ее технических приложениях. При не интегрируемой аналитически функции распределения искомые моменты можно определить предлагаемым численным методом . В общем случае область интегрирования Vможет иметь и бесконечную протяженность, при условии сходимости соответствующих интегралов (45), (46).

Возможные применения многомерных сферических координат для поиска экстремума на поверхностях и для некоторых других задач

Многомерные сферические координаты можно рассматривать как частный случай многомерных криволинейных координат. Переход к каким-либо криволинейным координатам от прямоугольных может определяться различными соображениями, в частности, удобством поиска экстремума функций и функционалов на поверхностях. Эти задачи оптимизации при ограничениях, встречающиеся, например, в теории управления и навигации [5], решаются методом неопределенных множителей Лагранжа [6, 7]. Однако, когда минимизируемая (или максимизируемая) функция характерно зависит от радиального направления, можно перейти к сферическим координатам и решать задачу прямым методом. При этом количество аргументов поиска экстремума становится равным числу степеней свободы в экстремальной задаче, т. е. устраняется неопределенность.

Пример. Найти минимальное и максимальное собственные числа матрицы квадратической формы и соответствующие им собственные векторы:

Г(Р) = 0,5 <Я,МИ>, (47)

где В = [х0, хь х2, х3]г - четырехмерный вектор прямоугольных координат;

4,1908 2,7032 -0,8016 1,2733

2,7032 5,7507 2,2245 0,9414

М «

-0,8016 2,2245 8,1563 -0,1829

1,2733 0,9414 -0,1829 9,2522

- матрица квадратической формы, которая сформирована на основе канонического жорданового разложения так, что ее собственные числа в точности равны

Х,= {1,530 6,240 9,210 10,37}, / = 0, 1, 2, 3.

При этом собственные нормированные векторы (| |е,| | = 1), соответствующие минимальному и максимальному собственным числам, равны

е0 = [0,731231 -0,614461 0,293668 -0,038707]г;

е3= [0, 327727 0, 495632 0, 318361 0,738641]г.

Подставляя в равенство(47)уравнения связи(2) и (3) при п = 4. получаем уравнение квадратической формы в виде

F(R, ocq, au a2) = 0,5R2FS,

гдеFs = Fs(a0, a1( a2) = <e, Me)-значение квадратической формы на единичной четырехмерной сфере;

е = [CS2CS-|CS(), CS2CS-,Sn0> CS2Sn1t sn2]T - единичный вектор направления, где cs,= cos ос,-, sn,= sin a,-, / = 0, 1, 2.

Квадратическая форма приводится к каноническому виду

F=0,5(X0u% + A., u2 + Х2 и2 + Х3 и2)- (48)

Дифференцируем дважды (47) и (48) и сопоставляем

Э 2F ^F_

^2 = Fs(a0. «1, a2) и э 2 = К (49)

Из формулы (49) можно заключить, что минимальное значение второй производной от квадратической формы по радиусу Fs будет получено тогда, когда единичный вектор направления е совпадет с направлением оси 0Ц-, которой соответствует минимальное собственное число Xj = X.mjn; максимум Fs также будет достигаться при аналогичном совпадении:

min{Fs(a0, сц, a2)} = Xmin и max{Fs(a0, a,, a2)} = Хтах.

Таким образом задача сводится к поиску экстремумов квадратической формы на единичной сфере. Минимум (максимум) несложно найти градиентным методом, соответственно как

V-т +1 — + кreg^ F.m ’

где a = [a0, 0^, a2]r - вектор углов;

VF«

ЭFSf dFS' dF5 da0 ‘ da1 ’ da2

- вектор первых произ-

водных (аналитические выражения для производ-

ных здесь не приводятся); кгед- регулировочный коэффициент, т - номер итерации.

В этой задаче минимальное и максимальное собственные числа находятся, например, за 16 и 32 итерации, соответственно, с относительными погрешностями

5а0) = 0,618 10“78(А3) = 0,119 10~7,

а нормированные собственные векторы находятся с угловой ошибкой в направлении

8е0 — 0,142 10“3 рад, бе, = 0,100-10“3 рад от нулевого приближения, в качестве которого было выбрано направление по оси 0Х0 (а0 = а1 = а2 = 0). Минимальное собственное число для корреляционной матрицы шума может использоваться, например, для расчета обобщенной дисперсии [8]:

где Е- энергия сигнала.

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

Следует упомянуть еще об одной интересной возможности приложения формул для направляющих косинусов в виде (3). В теории управления или навигации может встретиться задача аппроксимации траектории движения системы в многомерном пространстве состояний (или параметров), или аналогичная задача управления системой. При этом вектор состояния или параметров системы И! = [х0, Хт, х2,Хп_-|]гявляется функцией времени, где в общем случае могут быть неизвестны (или безразличны) как функциональная зависимость, так и отсчеты времени.

Введем некоторые аппроксимирующие функции от времени а0(0, 0^(0, а„_2(Г) таким образом, что

обобщенные скорости системы х/ = —- будут вы-

сУ?

ражены как направляющие косинусы

х0' = cos а0(0 cos a^{t)... cos an_2(t), х/ = cos a0{t) cos oti(f)... sin an_2(f),

Хн'^'П^).

При таком задании появляется удобное и полезное свойство

т

5(Т)= \M\dt =Т, (50)

о

где Б(Г) - путь, проходимый системой в пространстве состояний за условное время Г;

I М I = [(х0,)2 + (^1)2+--.+ (хп_1')2]1/2-модульскоро-сти. Известно, что одной траектории в пространстве может соответствовать множество функций времени Р(Г). Свойство (50) позволяет автоматически ввести нормировку и избавиться от неопределенности.

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

Автор выражает благодарность коллеге по работе - доценту кафедры №42 СПбГУАП С. Н. Воробьеву, который инициировал автора к написанию статьи и дал необходимые ценные советы.

1. Смирнов В. И. Курс высшей математики. - Т.2. - М.; Л.: ГИТТЛ, 1953. -627 с.

2. Возенкрафт Дж.5 Джекобе И. Теоретические основы техники связи. - М.: Мир, 1969. - 640 с.

3. Демидович Б. П., Марон И. А. Основы вычислительной математики. - М.: ГИФМЛ, 1968. - 659 с.

4. Кетков Ю. Л., Кетков А. Ю., Шульц М. М. , MATLAB 6. х: Программирование численных методов. - СПб.: «БХВ-Пе-тербург», 2004. - 662 с.

5. Алексеев В. И., КориковА. М., Полонников Р. И., Тарасенко В. П., Экстремальная радионавигация. - М.: Наука, 1978. - 279 с.

6. Булдырев B.C., Павлов Б. С. Линейная алгебра и функции многих переменных. - Л.: Изд-во ЛГУ, 1985. - 496 с.

7. Васильев Ф. П. Численные методы решения экстремальных задач. - М.: Наука, 1980. - 518 с.

8. Нестерук В. Ф. О влиянии формы сигнала на его обнаружение при нормальных коррелированных помехах // Радиотехника и электроника. -1963. - № 8. -С. 1319-1325.

9. Peter F., Swaszek, John В. Thomas Multidimensional spherical coordinats. Quantization // IEEE. Transfction of information theory. - Vol. lt-29. - N 4. - July 1983. - C. 570-576.