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

Вестник СПбГУ. Сер. 4, 2004, вып. 1

ДВУХЧАСТОТНАЯ ДВУХПОЗИЦИОННАЯ ФУНКЦИЯ КОГЕРЕНТНОСТИ ПОЛЯ СФЕРИЧЕСКОЙ ВОЛНЫ В ДИФФУЗИОННОМ МАРКОВСКОМ ПРИБЛИЖЕНИИ

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

К числу не решенных до конца проблем в диффузионном марковском приближении относится задача построения двухчастотной двухпозиционной функции когерентности. Общее точное аналитическое решение соответствующего уравнения имеет ряд трудностей и не найдено до сих пор. Численно эта проблема была исследована для степенного и гауссова спектра ионосферных флуктуаций [1]. Удалось построить точное аналитическое решение в случае распространения плоской волны .в случайной среде при квадратичной аппроксимации структурной функции [2]. В работе [3] этот результат был обобщен на случай распространения сферической волны. При этом структурная функция оставалась по-прежнему квадратичной. Авторы работ [4-7] предложили использовать метод разделения переменных и представили решение проблемы в виде разложения в ряд по поперечным собственным функциям при произвольной структурной функции. Однако подобный подход имеет некоторые трудности, связанные, например, со сходимостью рядов при некоторых начальных условиях, а также он, естественно, не работает, если переменные не делятся.

В работе [8] предложен асимптотический метод для построения функции когерентности в случае падения плоЯФЙ волны и произвольной структурной функции флуктуаций среды. Решение уравнения строилось по аналогии с классическим методом геометрической оптики (ГО), используя .представление о комплексных траекториях. В настоящей работе данный метод обобщается на случай распространения сферической волны в случайной среде с однородным фоном и флуктуациями, описываемыми произвольной структурной функцией. При применении квадратичной модели структурной функции предложенный метод приводит к результатам, полученным в [3].

2. Постановка задачи. Рассмотрим однородную среду (холодную изотропную плазму); занимающую все пространство. Пусть N.о - концентрация электронов. Пусть в левом полупространстве среда без флуктуаций, а справа с флуктуациями электронной концентрации на однородном фоне: А^г) = Л^+ДЛ^г) , где АИ - флуктуирующая добавка. Введем декартову систему координат. Границей раздела пусть будет плоскость 2 = 0. Поместим в левое полупространство, в точку с координатами (0,0, —г^, источник сферической волны. Будем считать, что флуктуирующие неоднородности обладают пространственной статистической однородностью. Рассмотрим случайное поле

ею

© А. А. Битюков, Н. Н. Зернов, 2004

УДК 537.86:519.2

А. А. Битюков, Н. Н. Зернов

*( \ ААГ(г) Ло

в котором г = {ж, у, г}.

Будем считать, что неоднородности крупномасштабны: Ц > 1, здесь - минимальный масштаб случайных неоднородностей в произвольном направлении (в дальнейшем просто /), а к - волновое число в фоновой среде. Используем для поля £(г) диффузионное марковское приближение, т.е. примем, что флуктуации электронной концентрации обладают дельта-корреляцией в направлении распространения [9].

Тогда, представляя распространяющееся поле в рде и = V{р,г)е*кх, где р = {х,у}, функция взаимной когерентности для комплексной амплитуды

г(Рг ~ Ръ г,кикг) = (у{р1,киг)у*{р2,к2,г)) будет подчиняться уравнению [3]

1 о к4 1 1

в котором кр = ипл- циклическая плазменная частота для фоновой концентрации Nо: и>„ч = гп и е - масса и заряд электрона соответственно; £о - диэлектричес-

кая проницаемость вакуума; к8 и кд - суммарная и разностная частоты: кя = , ка = к\ — к-2] Vс1 и У8 - дифференциальные операторы по поперечным разностным и суммарным координатам соответственно: р = рх — р2 ий = Р1+Р?; = {X, У}.

3. Преобразование уравнения. В случае падения плоской волны в направлении оси г функция когерентности Г не зависит от суммарных координат: У5Г = 0, и уравнение (1) упрощается. В рассматриваемом случае - падение сферической волны - производные по суммарным координатам остаются. Следуя работе [3], преобразуем уравнение (1). Будем искать функцию Г в виде Г =ТоГх, где Го - точное решение уравнения (1) в среде без флуктуаций (тогда А(0) = А(р) = 0):

1 (гк$ [Хх + Уу) + гка [1/2(Х2 + Г2) + 1/8(х2 + у2)] \ Го = т—■—г? ехР -:--- •

Подставляя Г = Г0Г1 в (1) и. учитывая, что Го есть его решение в среде без флуктуаций, получим уравнение для Г15 куда входят производные и Так как для Г! граничное условие .Гх^ -= 1 не зависит от Х и У, то и функция Г1 не зависит от суммарных координат: = = 0. В дальнейшем будем обозначать оператор градиента по разностным координатам без значка <1\ Vа = Ух- Так же как и в [3], представим Гх в виде Гх = Г2Г3, где

Г3 = ехр

■1/8А(0 )гк;

П2'

к h k

Подставляя Ti = Г2Г3 в уравнение (1) и полагая, что ks — та k2$ и ~ j^r , получим окончательно уравнение для функции Г2:

ж+++SÍA(0) - =(2)

которое для завершения постановки задачи дополняется граничным условием

4. Квазиклассическое приближение для уравнения (2). В цитировавшейся выше работе [3] было получено решение уравнения типа уравнения (2) лишь для квадратичной модели корреляционной функции флуктуаций. Применим для решения этого уравнения регулярную асимптотическую процедуру, разработанную в [8], которая позволяет строить решение для любых реалистических моделей корреляционной функции флуктуаций.

Для построения решения уравнения (2) используем метод формально аналогичный классическому методу ГО. Обозначим структурную функцию 2[Л(0) — А(р)], входящую в уравнение (2), как D(p):

D(p) = 2[A(0)-A(p)}.

Характерным масштабом изменения для нее будет величина I. Мы будем искать ' решение уравнения (2) в виде асимптотического ряда по малому параметру р:

Р = -j^j, ksl » 1.

Перейдем к безразмерным координатам:

- П = pksr\ zit = pkszt; г = {p,z}; n = {pi,*i}. В новых координатах уравнение (2) примет вид

7ЯГ+ - л)+ = (3)

где Vix ~ дифференциальный оператор по поперечным координатам рх. Как и в [10], будем искать решение уравнения (3)

Г-2 = и{ п)ехр[»ф(п)].

Функции U(t\) и Ф(г1) зависят от ks и kd- Функцию $(ri) удобно представить следующим образом:

*ы = Ms), .

Р

а функцию U (ri) разложим в ряд по степеням параметра р:

Гг = ехр

■ *i(r 1Г °°

т=0

Подставим ряд (4) в уравнение (3) и, приравнивая коэффициенты при одинаковых степенях р, получим уравнения, которые будут аналогичны уравнениям эйконала и

переноса соответственно в методе ГО. В исходной размерной системе координат и с учетом Ф = эти уравнения будут выглядеть так:

^ - а(УхФ)2 + • р) ~ ^ЬИ(р) = 0, (5)

дг г + гг 2

^- - 2о(Ух^оУ±Ф) - аЩУ 1Ф + —(Ъхи0 ■ р) = 0, (6)

ог г +

^-2а(УхС/тУхФ)-а[/тУ1Ф + -^— {Vхит ■ р) = ак$Ю\ит-и . (7) с/2; 2 +

к к* где а = 2*2"! Ь = т 6 N.

. Уравнение (5) является уравнением типа Гамильтона-Якоби. Для его решения формально применим метод характеристик [11]. Каноническая система в данном случае будет такой:

% = (9)

• (10)

■(РхР), (П)

др^

дт (2 + ^)2 Яф 1

— = р2 - 2а(рх)2 + —(рхР), (12)

здесь рж = дФ/дх, ру = дУ/ду, р3 = д$/дг, р/ = {рх,ру}-

Необходимо дополнить каноническую систему (8)-(12) начальными условиями. Обозначим начальное значение вектора р(т) как р0:

Р(0) = Р0 = {®о,Уо}.

Начальные условия для функции Ф (г) и вектора р± найдем из условия Гг | __0 = 1:

Фо = 0, Рох = 0.

Параметр ро, - начальное значение функций р5(т) - можно определить, подставляя в уравнение (5) величины векторов рл и р на граничной плоскости г = 0:

. рог = (г/2)Ш(р0).

Преобразовывая систему (8)-(12), можно получить

2 т

о

"íf = —iabV ±D(p), (15)

где уравнение (15) дополняется начальными условиями

Р( 0) = Ро. (16)

—р». а?)

drJT=о zt

Можно считать, что уравнение (15) является аналогом уравнения лучей в классическом методе ГО и определяет, вообще говоря, комплексную траекторию в трехмерном комплексном пространстве. Из равенства (13) следует, что параметр вдоль луча-траектории можно взять чисто вещественным. Величины р(0) = р0 являются комплексными координатами точки выхода луча-траектории с начальной поверхности Е, определяемой условием 2 = 0 (при г = 0). Точки выхода траекторий р0 должны подбираться таким образом, чтобы при фиксированном лучевом параметре ти траектория приходила бы в точку наблюдения рн с вещественными координатами: рн <Е R2[ll]. При этом при 0 < т < тн векторная функция р(т) принимает комплексные значения: р = р'+ гр".

Уравнение (14) определяет функцию Ф (т) вдоль луча, вышедшего из точки р0, где интеграл вычисляется вдоль луча-траектории р(т). Таким образом, решив уравнение лучей (15) и вычислив квадратуру в (14), найдем Ф как функцию лучевых координат: {р0, т}. Для того чтобы вернуться к исходным координатам г = {р, zj, необходимо разрешить систему уравнений р = R(p0Относительно р0. Причем, как уже отмечалось, при фиксированном значении т, равном z, необходимо выбирать такие комплексные значения р0, чтобы р было бы вещественным: р € R2. ■ Л ••■

Перейдем к вопросу об амплитуде U. Для нахождения членов ряда асимптотического разложения амплитуды необходимо разрешить систему рекуррентных уравнений (6)^(7). В дальнейшем для амплитуды U будем использовать только нулевое приближение, т.е. уравнение (6), которое с помощью уравнения лучей (15) приводится к обыкновенному дифференциальному уравнению. Решая его с учетом начального условия Z7o(0) = 0, найдем амплитуду Uq вдоль луча, т.е. представим Uq как функцию лучевых координат {р0,т}:

"£>(0)11/2

[/0 = l±ü

zt.

д(х,у)

[ D(r)

(18)

ад = ■

Считая, что в соотношениях (14), (17) и (18) zt —> оо, получим уравнение лучей с начальными условиями, функцию Ф и амплитуду Uq, соответствующие случаю падения плоской волны [8].

4.1 ■ Случай квадратичной структурной функции. Вначале посмотрим работу предлагаемого метода в случае квадратичной аппроксимации структурной функции. Считая неоднородности статистически изотропными, разложим функцию А(р) в ряд и ограничимся квадратичным членом этого ряда:

А(р) = А0 + А2р2.

Тогда структурная функция

D(p) = -2 А2р2. (19)

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

В случае структурной функции, выраженной формулой (19), переменные х и у в уравнении лучей (15) разделяются, и оно принимает простой вид

j2

— AiabAzp. (20)

атг

Решая уравнение (20) с начальными условиями (16) и (17), получим

р(т) = -i-(ztSchör + sMr)p0, (21)

ZtO

где 5 = 2{abA2)1/2ein/4.

В силу четности (21) как функции 5 ветвь корня в выражении для 5 может быть выбрана произвольным образом. Соотношение (21) описывает комплексную траекторию р(т) в четырехмерном пространстве р' + гр". Подставляя зависимость (21) в уравнение (14) и вычислив интеграл, можно найти функцию Ф вдоль траектории, т.е. представить Ф как функцию лучевых координат: {р0, т}. Для того чтобы вернуться к исходным координатам г = {р, z], необходимо разрешить соотношение (21) относительно р0, т.е. получить зависимость р0 — р0(р); р при этом следует брать вещественным - р G R2:

Ф(р, г) = -L {-L. _ + l р». (22)

' 4a\z + zt SztchSz + sMz J . '

Подставляя уравнение траекторий (21) в (18), можно найти определитель D(t) и амплитуду Uq(t). Возвращаясь к исходным координатам, для амплитуды нулевого приближения С/о будем иметь

Ц<»«> = Jcutl'L, _ <23>

Следует отметить, что в данном случае - случае квадратичной стуктурной функции - решения уравнений для амплитуд высших приближений в системе (6), (7) равны нулю. Тогда искомую функцию Г-2 можно записать следующим образом:

Г2 (г, kd,k,) = i/o (г) ехр[гФ(г)], (24)

где Uo(r) дается формулой (23), а Ф(г) есть выражение (22). Решение (24) совпадает с точным решением, приведенным в [3].

Таким образом, предлагаемый метод в случае квадратичной структурной функции D(p) = —2А-гр2 уже в первом приближении дает точный результат, а поправки, определяемые последующими приближениями, равны нулю.

Если дополнительно zt —> оо, то получим решение для плоской волны в среде с квадратичной структурной функцией, построенное в [2].

4-2. Общий случай статистически изотропных неоднородностей. КогДа флуктуирующие неоднородности статистически изотропны: А(р) = А(р), удобно перейти в уравнении (5) и системе (6), (7) к цилиндрическим координатам {р, <р, z}, совмещая ось z с направлением распространения волны. "Уравнение (5) для функции Ф в этих координатах будет иметь вид

89 (9Ф\2 1 дФ .1 , , ,

д7~аШ +7Т7/^-г2Ь1){р) = 0- (25)

Применив к уравнению (25) метод характеристик, получим систему уравнений, похожую на систему (13)—(15):

2 т *<т> = laThi ('(Т))2 - Чг " + */ДМ*,

й2г2 ... йр Начальные условия к уравнению (28) будут следующие:

' г = т, ' ' " (27)

(28)

Р(0) = Ро, (29)

т=О

Поскольку в цилиндрической системе координат радиус р не м'ожет быть отрицательной величиной, то в данном случае следует выбирать такие траектории, конечные точки которых не только вещественны, но и положительны. Решая задачу Коши - уравнение (28) с начальными условиями (29), (30) - находим

•> ол

(31)

где в = —2аЪ.

Ветвь корня в выражении (31) определим условием в точке ро:

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

Для расчета вдоль луча амплитуды С/о перейдем в уравнении (6) к цилиндрическим координатам и, используя вышеприведенные преобразования (см. п. 4), находим

1 1/2 ро 1 1

Р ОД

(32)

где £>(т) = ■

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

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

4.3. Случай экспоненциальной корреляционной функции. Одной из возможных реалистических корреляционных функций флуктуаций среды, для которых удается вычислить интегралы в формулах (26), (31), будет экспоненциальная функция вида

А{р) = <т21е~р/1. (33)

В этом случае при р оо структурная функция £>(р) стремится к положительной константе 2а21.

Решение соответствующей задачи Коши (28)-(30), задающей жомплексную траекторию р(т), в данном случае будет следующим:

р = ро + 211п (сЬ^ + , _ (34)

где <7 = Iй; 7 = у/21ва21е~р°/1 + д2. Переходя в4соотношения (34) к пределу гг -> оо, получим выражение для комплексных траекторий р(т), соответствующих задаче падения плоской волны, которое методом тождественных преобразований можно привести к формуле (34) работы [8].

Подставляя зависимость (34) в уравнение (26) и вычислив квадратуру, имеем для функции Ф

*<-»> = Ш " ^ +О+*+Ш-«)

(35)

-у т а -ут у т о "у7"

здесь 1и{т) = сЬ— + -вЬ —, и(т) = эЬ— + -сЬ —.

Выражение (48) задает функцию Ф вдоль луча, вышедшего из точки ро, т.е. представляет Ф как функцию лучевых координат {ро, г}. Амплитуда С/о будет определяться уравнением (32).

Как известно, функция 1п IV в случае комплексного переменного ш многозначная. Поэтому соотношение (34) следует понимать как аналитическое продолжение на комплексной плоскости -ш исходного элемента-ростка функции 1п ъи, определенного в некоторой окрестности Ую0 точки то условием 1пгио = 1п1 = 0, вдоль кривой 7, задаваемой комплеквнозначЦой функцией ш(т) вещественного переменного т, принимающей начальное значение ги(0) — гио = 1.

К сожалению, из соотношения (34) выразить ро невозможно. Поэтому в данном случае аналитически вывести выражение для функций Ф, С/о и соответственно для функции Г2 в исходных координатах {р, <р, г} не удается.

Переходя к пределу кл —»• 0 (а —»■ 0), получим случай одночастотной функции когерентности. Траектории здесь будут задаваться уравнением

р = ро(1 + ^

что соответствует решению уравнения (28) при а = 0. Функция одночастотной пространственной когерентности тогда будет следующей:

к$ л, [, , г + гь_-рц (л р/1-

' г + гг 4&2 р/1

)}]

(36)

где р = \р\ = л/х2 + у2] К = {X, У} - суммарные координаты (см. п. 2); к - волновое число в фоновой среде, соответствующее круговой частоте ш. Следует отметить, что такое решение отвечает нулевому члену асимптотического ряда (4), а поправки, определяемые остальными членами этого ряда, так же как и в случае квадратичной структурной функции, будут равны нулю.

В монографии [12] приведен метод построения точного решения для функции одночастотной пространственной когерентности, основанный на преобразовании Фурье исходного уравнения по суммарным координатам. Применив данный метод в случае корреляционной функции неоднородностей (33), получим точно такое же решение, т. е. функцию (36).

4-4- Рациональная корреляционная функция. Теперь рассмотрим другую функцию корреляции неоднородностей среды, также пригодную для описания умеренных и слабых флуктуаций поля. Это функция вида

2 >

л/ \ — а 1 А{Р> ~ 1 + а2(р/02'

Она имеет дополнительный масштаб а/1 в плоскости, поперечной направлению распространения волны. Структурная функция в данном случае

D(p)

2а la (р/1) 1 + а2(р/1)2 '

.При р —> оо структурная функция D(p), как и в предыдущем случае, будет стремиться к положительной константе 2а21.

Структурная функция D(p) на комплексной плоскости р будет аналитической в Точке р = 0 и бесконечно удаленной точке. Введем безразмерный параметр /3 = а{р/1). Разложим структурную функцию D в ряд Тейлора в окрестности точки /3 = 0 и ограничимся первым членом этого ряда:

D(p) « 2a2 la2 (р/1)2.

Структурная функция такого вида уже исследовалась (см. п. 4.1).

Теперь разложим функцию D в ряд по степеням 1/В и опять оставим только первый член этого разложения:

^'"^ф-гчЫ'- (37>

Применим предлагаемый асимптотический метод в случае структурной функции (37). Вычислив интеграл в формуле (31), для комплексной траектории т(р) получим

т =

/

{f/zt)2 + 2i6a2l \\\р0

■ (з8>

где / = а(ро/1). Ветвь корня в соотношении (38) определим в точке р = ро условием

/

zt'

Переходя к исходным координатам {р,г}, выразим ро из уравнения (38), задающего траектории, и, учитывая соотношение (27), получим алгебраическое уравнение восьмой степени относительно ро'-

(39)

Nt* + Mt3 + Pt2 + Ft + Q = 0,

в котором t = (po/l)2, P = 2i6a2l[l4- 2т/zt(1 + т/zt)], N = (al/zt)2(l + r/zt)2, F s= -2%Qo2p2/l,

M

-(ap/zt)2,Q = -

40 V

Способ построения решения уравнения (39) известен [13]. В случае падения плоской волны, при гь оо, уравнение (39) упрощается и становится биквадратным. Решив его, имеем выражение для ро, совпадающее с формулой (42) из работы [8].

Вернемся к случаю падения сферической волны. Функция Ф в исходных координатах {р, г} будет иметь вид

*> = Г-Т- - Ф- ('■-У + {1 + ^тЛ

Ааг + гг 4а \гь J \ а2(р0/1)2)

In z ( - +-^гЛ + I -In >,

4a a t .Л«« аРо / . . <*Ро.7 .)

(40)

где ?7 = 2aVableiv/4.

Амплитуда Uo будет задаваться соотношением

Uo{p,z) =

Z + Zt Zt

1/2

.р др/дрь.

(41) 31

В формулы (40), (41) входит величина ро как функция исходных координат {р, г), которая определяется из уравнения (39).

Рассмотрим случай одночастотной пространственной когерентности, и для этого перейдем к пределу ка -» 0 (о 0). Функция одночастотной пространственной когерентности, полученная предлагаемым асимптотическим методом в случае структурной функции (37) будет такой:

Г (р, R. z, к) = у——гх ехр (z + zty

z + zt

(42)

что совпадает с результатом, полученным строгим способом [12].

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

Работа выполнена при финансовой поддержке Министерства образования РФ (грант № Е02-3.5-128).

Summary

■ Bitjukov A.A., Zernov N.N. Two-frequency two-position coherence function of the field of a spherical wave in the diffusive Markov's approximation.

. In this work asymptotic method to solve Markov's parabolic equation for second order spaced position and frequency coherence function is reported. We consider the case of spherical wave propagation in a medium with homogeneous background and fluctuating inhomogeneities statistical properties of which are known. For obtaining asymptotic solution of the problem the complex geometrical-optics method is applied. The proposed asymptotic technique of solving Markov's parabolic equation works in the case of statistically anisotropic fluctuations of medium inhomogeneities as well as in the case of statistically , isotropic ones. Using this asymptotic technique we can obtain coherence function of the field for realistic structure functions of fluctuations. These enable to describe the cases of both weak and strong fluctuations. In the case of quadratic approximation of the structure function this method yields a previously known exact result.

Литература

1. Liu C.H., Yeh К. C. // Radio Science. 1975. Vol. 10, N 12. P. 1055-1061. 2. Sreenivasiah /., Ishimaru A., Hong S. Т. // Radio Science. 1976. Vol. 11, N 10. P. 775-778. 3. Knepp D.L. // Radio Science. 1983. Vol. 18, N 4. P. 535-549. 4. Oz J., Heyman E. // Radio Science. 1996. Vol. 31, N 6. P. 1907-1917. 5. Oz J., Heyman E. // Waves in Random Media. 1997. Vol. 7, N 1. P. 79-93. 6. Oz J., Heyman E. // Ibid. P. 95-106. 7. Oz J. // Ibid. P. 107-117. 8.'Bitjukov A. A,, Gherm E., Zernov N. N.//Radio science. 2002. Vol. 37, N 4. Art. 1066. 9. Liu С. H., Yeh К. C. U Proc. IEEE. 1982. Vol. 70, N 4. P. 324-360. 10. Кравцов Ю. А., Орлов Ю. И. Геометрическая оптика неоднородных сред. М., 1980. 11. Кравцов Ю.А. // Аналитические методы в теории дифракции и распространения волн/ Под ред. С. В. Бутаковой. М., 1970. С. 257-363. 12. Рытое С. М., Кравцов Ю. А., Татарский В. И. Введение в статистическую радиофизику: В 2 ч. Ч. II: Случайные поля. М., 1978- 13. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М., 1978.

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