У ДК: 576.8.21 ОСОБЕННОСТИ ВЫЧИСЛЕНИЯ ИНТЕГРАЛОВ ФОКА В.И. Вьюнник, А.А. Звягинцев Харьковский национальный университет имени В.Н. Каразина, Харьков, пл. Свободы 4, 61077 Поступила в редакцию 10.12.2008 В работе рассматривается вычисление интегралов Фока, имеющих прикладное значение в задачах дифракции на выпуклых телах. Описан метод вычисления интегралов и проведена оценка их сходимости. Определен диапазон значений комплексного параметра, характеризующего поверхностный импеданс, для которого применим данный метод . Проведена оценка погрешности вычислений. КЛЮЧЕВЫЕ СЛОВА: интеграл Фока. численное интегрирование, погрешность. Интегралы Фока широко используются в решении задач дифракции на выпуклых телах. Впервые они были введены Фоком для решения задачи распространения радиоволн над земной поверхностью [1]. В дальнейших своих работах (см. в сборнике [2]) Фок обобщил их применение на случай произвольного выпуклого тела. Замечательным выводом Фока является обнаруженное им свойство локальности поля на поверхности в полутеневой зоне. А именно, значение поверностной плотности тока, выражаемой через интеграл для точки, лежащей вблизи границы свет-тень, определяется значениями подинтегральной функции в окрестности данной точки. Отсюда следует, что любые тела с плавно изменяющейся кривизной, будут иметь в полутеневой зоне при совпадении их радиусов кривизны одинаковые значения токов на их поверхности. Таким образом, решив задачу для дифракции плоской электромагнитной волны на параболоиде, Фоком были получены универсальные формулы для полей на выпуклых телах, связывающие значение полного поля на поверхности тела с полем падающей волны, умноженным на некоторую универсальную функцию, названную Фоком множителем ослабления. В зависимости от поляризации падающей волны, множитель ослабления может быть выражен через введенные Фоком G и F функции, которые принято называть g и f интегралами Фока. Несмотря на то значение, которое имеют эти интегралы для практических приложений, вычислению собственно самих интегралов уделено недостаточно внимания. Помимо таблиц значений, приведенных в [2], на русском языке существует, насколько нам известно, только одна работа [3], подробно описывающая расчет интегралов Фока. Однако в ней рассматривается вычисление интеграла при q = 0, и не приведены формулы для более общего случая, когда q = 0. В связи с этим не исследован вопрос о границе применимости самого метода вычисления интеграла. Наиболее полно вопрос вычисления интегралов Фока описан в работе [4]. В ней исследуется вопрос о сходимости интегралов и обосновывается возможность применения Filon-trapezoidal метода численного вычисления этих интегралов. В работе также исследуется вопрос точности вычислений, причем для оценки верхнего предела интегрирования используется аналитический подход. При всей его простоте и изяществе, в то же время он не позволяет оценить влияние параметра q на ошибку вычислений. В настоящей работе мы попытались охватить вышеуказанные проблемы. ОПРЕДЕЛЕНИЯ И СВОЙСТВА Будем называть интегралом Фока интеграл вида 1 F (z, y, q ) = √ π Γ eizt Φ(t, y, q ) dt. (1) Здесь Φ(t, y, q ) = v (t − y ) − либо Φ(t, y, q ) = ′ ′ v (t) − qv (t) w1 (t − y ) ′ w1 (t) − qw1 (t) ′ ′ (2) i 2 w2 (t − y ) − w2 (t) − qw2 (t) w1 (t − y ) , ′ w1 (t) − qw1 (t) (3) где v , w1 , w2 , w1 и w2 -функции Эйри и их производные в нотации Фока (называемые функциями ЭйриФока), причем v = (w1 − w2 )/2i, а контур интегрирования Γ на комлексной плоскости проходит вдоль 2π луча ei 3 от ∞ до 0 и от 0 до ∞ вдоль действительной оси ОХ. Этот интеграл описывает поле вблизи поверхности выпуклого тела (на расстоянии y от нее), и может непосредственно быть получен в задачах дифракции на круговом цилиндре и сфере, путем асимптотического выражения функций Бесселя и Ханкеля через функции Эйри [5]. При y = 0 интеграл (1) может быть представлен в виде 1 F (z, q ) = √ π Γ eizt dt. w1 (t) − qw1 (t) ′ (4) Для случая идеально проводящего тела, Фоком также были получены формулы для касательных составляющих магнитного поля на его поверхности. Эти составляющие могут быть выражены через интеграл 1 g (z ) = √ π Γ eizt dt, ′ w1 (t) (5) в случае H -поляризации и через интеграл 1 f (z ) = √ π Γ eizt dt w1 (t) (6) в случае E -поляризации. Интегралы (5) и (6) носят название g – и f – интегралов Фока. Несложно заметить, что интеграл (4) превращается в g -интеграл Фока, если положить q = 0. В то же время, если cчитать, что q = ∞ то интеграл (3) сведется к интегралу 1 F (z, q ) = √ π Γ − eizt dt, q w1 (t) (7) который на поверхности равен нулю. Видно, что выражение (7) в −q раз меньше f –интеграла Фока, что связано с тем, что эти интегралы представляют поля соответственно E и H касательных компонент электромагнитного поля. Касательные же компоненты, в случае E –поляризации, связаны соотношением Леонтовича-Щукина Etg = −δ = − Htg µ ε (8) ka 3 Здесь m = , k –волновое число, a–радиус кривизны поверхности в данной точке . Магнитная 2 проницаемость принята равной единице. Входящая в (9) диэлектрическая проницаемость представляет в общем случае комплексную величину ε=ε +i ′ ′ где δ –приведенный поверхностный импеданс, а µ и ε–соответственно магнитная и диэлектрическая проницаемости. Отсюда видна физическая суть параметра q как величины, связанной с поверхностным импедансом. Этот параметр равен  1  imδ = im √ в случае H -поляризации; ε (9) q=  im 1 = im√ε в случае E -поляризации. δ 1 4πσ . ω (10) В выражении (10) величина ε есть диэлектрическая проницаемость среды, обусловленная эффектом поляризации (токами смещения), а второе слагаемое-величина, обусловленная токами проводимости (σ представляет электропроводность, а ω –частоту). Можно заметить, что с учетом возможных значений диэлектрической проницаемости, электропроводности и частоты, величина δ -комплексная величина, с положительной действительной и отрицательной мнимой частью. Более точное исследование диапазона возможных значений опирается на тот факт, что при отсутствии источников поля внутри облучаемого тела, поток мощности должен быть направлен внутрь тела, а не наоборот. Несложно показать (см. например [6]), что в этом случае Re(δ ) 0 т.е. целиком расположена в правой полуплоскости комплексной плоскости. Следует заметить, что значениям δ расположенным вне области −π/4 < arg (δ ) 0 соответствуют σ и ′ ε которые в этом случае могут принимать отрицательные значения. Как показано в [6], такие импедансы можно реализовать используя многослойные (двухслойные) покрытия. С учетом вышеизложенного и анализируя (9) несложно заметить, что параметр q может принимать значения, соответствующие 0 arg(q ) π . Функции Эйри-Фока, входящие в состав подинтегральных выражений, представляют собой интегралы вида 1 w1,2 (t) = √ π Γ1,2 2π etz− 3 z dz 1 3 (11) где Γ1 представляет собой контур от ei 3 ∞ до 0 и от 0 до ∞ вдоль действительной оси OX а Γ2 – его зеркальное отражение от действительной оси. Они, вместе с "классическими"функциями Эйри (функциями Эйри в нотации Миллера) являются решениями дифференциального уравнения w (t) = tw(t) и могут быть представлены в виде w1 (t) = u(t) + iv (t); w2 (t) = u(t) − iv (t), где u и v выражаются через функции Эйри-Миллера как (13) ′′ (12) u(t) = √ πBi(t); v (t) = √ πAi(t), (14) где Ai и Bi функции Эйри первого и второго рода (которые в общем случае комплексных значений t имеют комплексный характер). Удобными формулами, связывающими эти функции являются: 2π π √ (15) w1,2 (t) = 2 πe±i 6 Ai te±i 3 для функций и π 2π 2π √ ′ ′ w1,2 (t) = 2 πe±i 6 e±i 3 Ai te±i 3 (16) для их производных. Для целей нашего исследования полезными будут также соотношения, связывающие функции Эйри друг с другом (так называемые соединительные формулы): w1 (tei 2π 3 ) = ei 3 u(t) − iv (t) = ei 3 w2 (t); w1 (tei π π 4π 3 ) = 2ei 6 v (t). π (17) Подробное описание функций Эйри и их свойств может быть найдено в [7] ВЫЧИСЛЕНИЕ ИНТЕГРАЛОВ ФОКА Для вычисления интегралов Фока могут использоваться различные методы. Возможно как непосредственное численное интегрирование, так и представление интеграла в виде суммы вычетов. В данной работе используется второй подход. Он применим, если z не является большим отрицательным числом (реально при z > −3). Рассмотрим интеграл (4). Чтобы перейти к численному интегрированию, нам необходимо преобразовать путь интегрирования, чтобы он целиком был расположен на действительной оси. Для этого представим интеграл в виде суммы двух интегралов: F (z, q ) = F1 (z, q ) + F2 (z, q ), где 0 (18) 1 F1 (z, q ) = √ π и ei2π/3 ∞ ∞ eizt dt; w1 (t) − qw1 (t) ′ (19) 1 F2 (z, q ) = √ π 0 eizt dt. w1 (t) − qw1 (t) ′ (20) Интеграл (20) есть интеграл с вещественными пределами. Первый же интеграл мы преобразуем к ним. Осуществляя замену переменной t = t1 e i 2π 3 (21) и учитывая первое соотношение соединительных формул (17), получим 0 1 F1 (z, q ) = − √ π eizt1 e ′ i 2π 3 2π 3 ∞ w2 (t1 ) − qei w2 (t1 ) 1 dt1 = √ π ∞ ′ eizt1 e i 2π 3 2π 3 0 w2 (t1 ) − qei dt1 . w2 (t1 ) (22) Аналогичным образом может быть приведен к интегралам с вещественными пределами интегрирования и интеграл (1). Причем, как показано в [3] выражение (3) следует использовать для пути интегрирования от ei 3 ∞ до 0, а выражение (2) для пути от 0 до ∞. Мы приведем здесь конечный результат преобразования интеграла (1) к вещественным пределам 2π 1 F1 (z, y, q ) = √ π ∞ e 0 2π i izte 3 v (t1 − 2π ye−i 3 ) − v (t1 ) − q ei ′ 2π 3 v (t1 ) w2 (t1 ) w2 (t1 ) − ′ 2π q ei 3 w2 (t1 − ye−i 2π 3 ) dt. (23) Оценим насколько пригодны (20), (22) для численного интегрирования. Можно заметить, что числитель подинтегрального выражения (20) представляет собой периодическую (тригонометрическую) величину и, следовательно, сходимость интеграла определяется знаменателем. Знаменатель же представляет в общем случае разность двух быстрорастущих величин. Чтобы оценить скорость роста этих величин, воспользуемся асимптотическими представлениями функции Эйри-Фока и ее производной. Имеем w1,2 (t) ∼ ′ e3t 3 2 2 ′ 1 ; w1,2 (t) = e 3 t t 4 . 3 2 2 1 (24) t4 Величина w1 (t) растет быстрее чем w1 (t). Несмотря на это, несложно заметить, что при определенных q разность будет равна 0, что приводит к бесконечно большому значению подинтегрального выражения. Что же касается интеграла (22), то замена переменной приводит к тому, что показатель экспоненты содержит в себе действительную часть. При отрицательных значениях z это может приводить к росту подинтегрального выражения. Асимптотическая оценка подинтегрального выражения (22), выполненная без учета влияния параметра q (например, приравнивая q нулю), приводит к следующему выражению integrand = 1 e 2 izt e− √ 3 3 2 2 2 zt+ 3 t t− 4 . 1 (25) Зависимость (25) от t для различных отрицательных z представлена на рисунке 1. Видно, что сначала преобладает рост, за счет действительной части показателя экспоненты, но затем начинает доминировать асимптотика Эйри функции, что приводит к быстрому уменьшению подинтегрального выражения. Отметим также сильную зависимость первоначального роста от величины отрицательного значения z . Рис. 1. Зависимость асимптотического представления подинтегрального выражения (25) от t для различных значений величины z . Вернемся к оценке зависимости величины знаменателя от параметра q . Для сходимости интеграла ситуация когда знаменатель обращается в нуль должна быть исключена. Для нахождения соответствующих условий, необходимо найти нули уравнения w1 (t) − qw1 (t) = 0, ′ (26) т.е. найти его корни ts (q ) где s- номер корня ( s = 1, 2 . . . ). Прежде всего заметим, что при q = 0 уравнение (26) сведется к уравнению w1 (ts0 ) = 0, а при q = ∞ к уравнению w1 (ts∞ ) = 0, ′ (27) (28) где ts∞ и ts0 –нули функции Эйри-Фока и ее производной соответственно. С учетом выражений (15), (16) можно показать, что нули функции Эйри-Фока и ее производной связаны с нулями функции ЭйриМиллера Ai и ее производной следующим соотношением tsw1,2 = −tsAi e±i 3 ′ π (29) ′ т.е все нули w1 и w1 расположены на луче ei 3 , а нули w2 и w2 расположены на луче e−i 3 . С учетом (12), дифференцируя (26) несложно получить дифференциальное уравнение для нахождение корней ts как функции от q : dt 1 . = dq t − q2 Подобное же дифуравнение dt ei 3 = 4π , dq t − q 2 ei 3 можно получить и для уравнения w2 (t) − qei Решая эти уравнения с начальными условиями: t(q = 0) = ts0 ; или t(q = ∞) = ts∞ , (34) (33) ′ π π (30) 2π (31) 2π 3 w2 (t) = 0. (32) для различных q , получаем функцию ts (q ), которая является решением этих уравнений. Решения уравнений (30) и (31) представлены на рисунках 2 и 3. На них изображена комплексная плоскость значений ts (q ), причем непрерывными линиями изображены ts (q ) при arg(q ) = const, а точечными – ts (q ) при abs(q ) = const. Представлены решения для первых двух корней уравнений, соответствующих t10 , t1∞ и t20 , t2∞ . Из рисунка 2 видно, что случаю arg(q ) = 0 соответствует кривая (значений корней), начинающаяся с корня t10 и уходящая на бесконечность по вещественной оси OX . Поскольку интегрирование производится по t вдоль этой самой оси, в случаях, когда arg(q ) = 0 (т.е. когда q –чисто действительная величина) указанные выше формулы вычисления интеграла Фока неприменимы. Что же касается знаменателя интеграла (22), то его корни, как видно из рисунка 3, нигде не пересекают путь интегрирования (полуось OX ). Рис. 2. Корни уравнения (26) как функция параметра q . Рис. 3. Корни уравнения (32) как функция параметра q . ОЦЕНКА ПОГРЕШНОСТИ ВЫЧИСЛЕНИЙ Выражения (20), (22) и (23) представляют собой интегралы от 0 до ∞. Численное интегрирование в полубесконечных пределах хотя и возможно, но не всегда оправдано. Представляется важным оценить погрешность вычислений, связанную с усечением верхнего предела интегрирования. Одним из возможных подходов является метод аналитической оценки, изложенный в работе [4]. Будем рассматривать интегрирование интеграла (22) в интервале от 0 до A. Величина остатка представляет собой интеграл 1 ǫ(a) = √ π ∞ ′ eizte i 2π 3 2π 3 A w2 (t) − qei dt, w2 (t) (35) величина которого определяет ошибку, вызванную усечением. Чтобы оценить ее, заменим функцию Эйри и ее производную их асимптотиками, и принимая q = 0 получим выражение, аналогичное (25) (мы откинули экспоненту с мнимым показателем, так как она не оказывает влияния на рост интегранда) 1 ǫ(a) = √ π ∞ e A − √ 3 3 2 2 2 zt+ 3 t t− 4 . 1 (36) Поскольку величина ошибки не должна зависеть от z , а как было показано ранее (см. рисунок 1) при отрицательных значениях z первоначально происходит даже рост подинтегрального выражения, поступим следующим образом. В асимптотическом представлении подинтегрального выражения удалим составляющую, зависящую от z , и введем новую величину B , которая и будет являтся нижним пределом интегрирования 1 ǫ(a) = √ π ∞ e− 3 2 2 3t t− 4 . 1 (37) B При положительных z , интеграл (37) сходится несколько медленнее чем (36) (в предположении, что B = A, и значит будет несколько завышать ошибку, а при отрицательных z он сходится быстрее, и будет уменьшать (занижать) ее. Чтобы компенсировать это уменьшение, будем считать, что при z > 0, B будет равняться A, а при z < 0, B должно равняться расстоянию между A и точкой пересечения асимптотической кривой оси OX . Отсюда для B имеем:  A − 27 z 2 ; при z < 0. (38) B= 16 A; при z > 0. Интеграл (37) может быть легко приведен к выражению ǫ(A) = 3 1√ 1√ 6 1 − erf 6B4 3 3 , (39) которое позволяет легко определять ошибку усечения для любого, наперед заданного A. Другой возможный подход состоит в том, чтобы применить численную оценку ошибки усечения. Как известно, интеграл с бесконечными пределами интегрирования, путем соответствующей замены переменной может быть приведен к интегралу с конечными пределами, который затем вычисляется численными y методами. Например, замена t = преобразует пределы интегрирования от 0 − ∞ к 0 − 1. 1−y Сравнение значений интеграла полученного таким методом, со значениями интеграла, полученными посредством усечения верхнего предела позволяет оценить ошибку усечения. Этот метод также позволяет оценить влияние значений параметра q на ошибку усечения. (а) (б) Рис. 4. Абсолютное значение ошибки усечения верхнего предела интеграла (20): (а) верхний предел = 5, (б) верхний предел = 6. На рисунке 4 приведены результаты зависимости от q = x + iy такой оценки для интеграла (20), для случая усечения верхнего предела к t = 5 и t = 6 соответственно (z = 1). Легко видеть, что случаям малых y (y → 0) соответствуют максимальные ошибки вычислений, которые в окрестности корней не позволяют оценить не только ошибку усечения но и сами интегралы. (а) (б) Рис. 5. Абсолютное значение ошибки усечения верхнего предела интеграла (22): (а) верхний предел = 5, (б) верхний предел = 6. На рисунке 5 представлены результаты зависимости оценки усечения от q = x + iy для интеграла (22) (верхний предел t = 5 и t = 6 соответственно). Можно заметить, что ошибка усечения возрастает при y → 0, особенно в области x ≈ −1 , хотя и не достигает таких значений, как в случае первого интеграла. ИНСТРУМЕНТЫ, МЕТОДИКА РАСЧЕТА И РЕЗУЛЬТАТЫ Расчет интегралов Фока может быть реализован в различных программных средах. Наиболее удобным представляется использование для этих целей программных продуктов типа Matlab или Octave. Реализованные нами программы в Matlab использовали для численного интегрирования встроенную процедуру quad. Она использует метод адаптивной квадратуры Симпсона. Численное интегрирование проводилось в пределах от 0 до 20. Сравнение с табличными данными, опубликованными в [2] показало полную идентичность. На рисунке 6 приведены результаты вычисления по формулам (20), (22) модуля и фазы интеграла Фока π для различных значений комплексного параметра q (при arg(q ) = ). 3 Как видно из рисунка, параметр q оказывает существенное влияние на численные значения как модуля так и фазы интегралов Фока. Для исследования погрешности вычислений и обеспечения наперед заданной точности использовалась система компьютерной алгебры Maple, обладающей широкими возможностями не только символьных но и численных вычислений. Численное интегрирование проводилось с использованием адаптивного double-exponential метода (Dexp). Для поиска нулей уравнений (26) и (32) комплексная переменная q была представлена в виде произведения действительного и мнимого числа q = αq и уравнения решались относительно действительной переменной α (подход, описанный в [8]). Таким образом, вместо уравнений (30) и (31) решались уравнения (а) (б) Рис. 6. Зависимость модуля и фазы интеграла Фока от z для различных значений |q | (arg(q ) = π/3): (а) модуль интеграла, (б) фаза интеграла. q dt ; = dα t − α2 q 2 и dt qei 3 = 4π ; dα t − α2 q 2 ei 3 2π (40) (41) соответственно. Уравнения решались численным методом, с использованием Maple функции dsolve. Данная функция позволяет использовать различные методы решения ODE . Нами использовался Livermore Stiff ODE метод. Заметим, что данные дифференциальные уравнения крайне чувствительны к точности вычислений, что связано с чрезвычайно быстрым ростом как самой функции Эйри-Фока, так и ее производной. В силу этого, вычисление нулей функции Эйри и ее производной, которые являются начальными граничными условиями для уравнений (40) и (41) желательно выполнять с максимальной точностью. Для их вычисления нами была использована программа на языке FORTRAN, реализующая алгоритм вычисления функций Эйри и производных, а также их нулей с точностью до 34 знаков (quadruple precision), описание которого приведено в [9]. ЗАКЛЮЧЕНИЕ Реализован алгоритм вычисления интегралов Фока при ненулевых значениях комплексного параметра q . Исследованы границы применимости метода численного интегрирования с учетом различных значений указанного параметра и проведена оценка точности метода. 1 СПИСОК ЛИТЕРАТУРЫ [1] Фок В.А. Дифракция радиоволн вокруг земной поверхности. — Москва: Изд. АН СССР, 1946. [2] Фок В.А. Проблемы диффракции и распространения электромагнитных волн. — Москва: Советское радио, 1970. [3] Иванов В.И. Табулирование функций В.А. Фока // Записки научных семинаров ЛОМИ:Математические вопросы теории распространения волн. — 1979. — Т. 89, № вып. 10. — С. 97–111. — Ленинград, Наука. [4] Pearson L. A scheme for automatic computation of fock-type integrals // IEEE Transactions on Antennas and Propagation. — 1987. — Vol. AP-35, no. 10. — Pp. 1111–1118. [5] Syed H., Volakis J. An asymptotic analysis of the plane wave scattering by a smooth convex impedance cylinder: Tech. Rep. 025921-9-T: University of Michigan Radiation Laboratory, 1990. — January. — 31 pp. [6] Гюннинен Э.М., Макаров Г.И. Поле точечного диполя над импедансной поверхностью // Проблемы дифракции и распространения волн. — 1966. — no. вып. 5. — Pp. 97–120. — Л., Изд-во Ленингр. ун-та. [7] Frank W. J. Olver. Chapter 9. Airy and Related Functions. — NIST Digital Library of Mathematical Functions. — http://dlmf.nist.gov/. [8] Гюннинен Э.М., др. Распространение электромагнитных импульсов и их гармонических составляющих над земной поверхностью // Проблемы дифракции и распространения волн. — 1964. — no. вып. 3. — Pp. 5–191. — Л., Изд-во Ленингр. ун-та. [9] B. R. Fabijonas. Algorithm 838: Airy Functions. — ACM Transactions on Mathematical Software (TOMS). — http://portal.acm.org/citation.cfm?id=1039819.