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

УДК 621.371

А.А. Сенченко, Ю.П. Саломатов

Сравнение квадратурных методов решения интегрального уравнения Хаффорда

Получены формулы численного решения интегрального уравнения Хаффорда при использовании различных квадратурных формул Ньютона-Котеса. Приводятся результаты численных расчетов функции ослабления над различными однородными трассами для сферической модели Земли. Проведено сравнение результатов численных расчетов с результатами расчета функции ослабления по ряду нормальных волн и приближенным решением интегрального уравнения Хаффорда для малых расстояний.

Ключевые слова: численное решение интегральных уравнений, интегральное уравнение Хаффорда, интегральное уравнение Фейнберга, квадратурные формулы Ньютона-Котеса, функция ослабления, распространение радиоволн.

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

[2] и в других областях [3].

Решение интегрального уравнения Хаффорда является не единственным методом учета влияния свойств подстилающей поверхности на распространяющуюся вдоль поверхности Земли волну. Для однородных трасс можно использовать формулу для плоской Земли [4], формулу Фока [5], решение интегрального уравнения Хаффорда в виде ряда по полуцелым степеням [6]. Для кусочнооднородных трасс возможно использование формул Калинина и Фейнберга [7]. Также во всех перечисленных случаях можно использовать численное решение интегрального уравнения Фейнберга для трасс с электрофизическими неоднородностями [7] и трасс с электрофизическими и геометрическими неоднородностями [8]. Выбор того или иного метода расчета зависит от выбранной модели трассы распространения радиоволн и её протяженности. Наиболее универсальным методам при этом является численное решение интегральных уравнений Фейнберга или Хаффорда, чему посвящено много работ (см., например, [9, 10]). Оба уравнения относятся к интегральным уравнениям

Вольтерра второго рода и в своем ядре содержат множитель 1/д/5(х-5), где 5 - переменная интегрирования, х - точка, в которой вычисляют значение функции ослабления. Поэтому к ним применимы одни и те же методы численного решения, но ядро интегрального уравнения Хаффорда не требует вычисления функции ослабления над эталонной трассой, что делает его более удобным для разработки алгоритма численного решения, по сравнению с интегральным уравнением Фейнберга. По этой причине в рамках данной работы мы ограничились решением только интегрального уравнения Хаффорда.

Интегральное уравнение Хаффорда для модели гладкой однородной сферической Земли можно записать как:

, ,| . 5 . X-5 . X

I у- х У 2акI вт-----Нет---------вт—

Ж(х) = 1 + у Ц- (Ж(я)е [ 2а 2а 2а

'2п 0

, . X - 5

о + вт---------

2а

Ж, (1)

где х - расстояние по прямой между точками излучения и наблюдения; а - радиус Земли.

Использование данной модели позволит нам сравнить результаты численного решения интегрального уравнения Хаффорда с результатами расчетов по формулам для однородной сферической Земли.

Для уменьшения объемов вычислений при решении интегрального уравнения Хаффорда воспользуемся следующей аппроксимацией показателя экспоненты в (1) [11]:

. ( 5 | . (х - 5| . ( х | 5х(х - 5)

81п! —1 + 8ш|—— I-вт I —I*-----------------------------------------------—. (2)

У2а) У 2а ) У2а) 16а3

В этом случае ядро интегрального уравнения будет иметь вид

5х(х-5)

У- 2

е 8а

(3)

* ')■ у

-^s(Ууs)

Разница между значениями, получаемыми с помощью точного выражения ядра и с использованием аппроксимации (2), как показывают расчеты, не превышает 0,052 % на расстояниях до 1000 км для любых трасс. Поэтому вместо точного выражения ядра интегрального уравнения мы будем использовать формулу (3).

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

( ,--1 Л

Ж ■

/(1-ХЛцЯц), (4)

/, +^Х ЛцКуЖу

У У1

где Д, = Дх,), х, = х0 + (, - 1)к,, = 1, 2, ..., N к - шаг интегрирования; N - число точек интегрирования; х = х^ Ку = К(х,, 5у), Ж, = Ж(х,), Лу - набор коэффициентов для выбранной квадратурной формулы интегрирования. Первое значение функции определяется соотношением Ж0 = /0.

Ядро интегрального уравнения Хаффорда (как и ядро интегрального уравнения Фейнберга [7]) имеет особенности на концах промежутка интегрирования, что не позволяет напрямую использовать квадратурный метод для его численного решения. Причиной этих особенностей является знаменатель ^'(х - 5) в ядре интегрального уравнения. Для применения квадратурного метода численного интегрирования интеграл в (1) разобьём на три части и в первом интеграле сделаем замену переменной х = у2, а в третьем х = В - у2 [4]. Тогда будем иметь:

х ^гцку (х, у2) щк л/х-и2к Г (х,х - у2) . ,

(К(х,5)Ж(s)ds■ 2 ( . =Ж(у2\<Су + ( К(х,5)Ж+ 2 ( —. Ж(х-у2), (5)

0 0 -\/х - у2 щк 0 а/х-у2

где к - величина шага интегрирования; щ, п2 - число точек, по которым вычисляется первый и последний интегралы соответственно. Числа щ и п2 должны удовлетворять условию 1 < щ < п2 < N;

г (x, 5)=* (x, х - 5) - та часть ядра интегрального уравнения, в которой нет указанных особых

точек.

Подставляя (5) в (1), получим интегральное уравнение без особых точек на краях интервала интегрирования и сможем использовать формулу (4) для его численного решения. Теперь остается только выбрать квадратурную формулу для вычисления каждого из интегралов и определиться с величиной целочисленных переменных п1 и п2.

Для краткости введем обозначения для подынтегральных функций первого и третьего интегралов из (5) без учета искомой функции ослабления Ж(х):

К1 (х,?) = 2 г (х,?)/у[х-1, (6)

К3 (х,?) = 2г(х,х- ?)1'1х-, (7)

где К](х, ?) относится к первому интегралу; К3(х, ?) - к последнему, а параметр ? = у2 в обоих случаях.

Для численного решения интегрального уравнения Хаффорда (1) может быть использована квадратурная формула Ньютона-Котеса первой степени (формула трапеций) [4, 9].

Формула Ньютона-Котеса первой степени обладает низкой точностью в сравнении с другими квадратурными формулами более высокого порядка. Рассмотрим квадратурную формулу Ньютона-Котеса второй степени (формулу Симпсона). Для её использования необходимо иметь нечетное количество точек, но при вычислении каждого следующего значения искомой функции в интегральном уравнении число интегрируемых точек будет то четным, то нечетным. Поэтому при четном числе точек будем комбинировать квадратурные формулы Ньютона-Котеса второй и третьей степеней, причем последнюю будем использовать только для последних четырех точек. Следует заметить, что для крайних интегралов в (5) использовать квадратурные формулы с постоянным шагом интег-

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

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

1-^-2 (2,0)

(8)

где г = 4, 5, ... N а первые четыре точки могут быть вычислены с использованием формулы (10). Переменная Щв функции Я(Щ, т1, т2, К) представляет собой массив значений функции ослабления, вычисленных на предыдущих шагах. Функция Я(Щ, т1, т2, К) вычисляет значение среднего интеграла в (5), комбинируя квадратурные формулы Ньютона-Котеса второй и третьей степеней:

К т1 +пр - шоё(и-1,2)-1 Г 0,шоа(и -1,2) = 0,

S (7 ,ть т2,Щ = - X (у + (Ъ] +1 + Ъ] +2 )

J =т

3К (9)

— ( + 37к+1 + 3¥к+2 + ¥к+3 ),

8

где 7 - массив со значениями подынтегральной функции; т\ и т2 - начальный и конечный индексы в массиве 7, в пределах которых требуется вычислить интеграл; п = т2 - т\ + 1 - число точек интегрирования; Пр =[(п -1)/2] - число промежутков интегрирования, в котором квадратными скобками

обозначена целая часть от деления; К - шаг интегрирования; к = т\+2(пр - шо^п - 1, 2)), шо^х, у) -обозначает остаток от деления х на у. Разница т2 - т\ должна быть не меньше 2, поскольку в противном случае использовать формулу Ньютона-Котеса второй степени нельзя.

Формула (8) должна обладать большей точностью, чем решение с использованием формулы Ньютона-Котеса первой степени, но её можно повысить, если крайние интегралы в (5) вычислять с использованием квадратурной формулы более высокого порядка. Для этого необходимо использовать квадратурную формулу с переменным шагом. Для её получения воспользуемся методикой, приведенной в [13], в результате чего получим выражения коэффициентов для квадратурной формулы Ньютона-Котеса второй степени с переменным шагом:

^2,0 (х0,Х1,Х2) = (х2 - Х0)(2х0 -3x1 + Х2)/[б(х0 - Х1)] , (10.1)

С2,1 (,Х1,Х2) = (Х0 -Х2)3/[б(х1 -Х0)( -Х2)] , (10.2)

С2,2 (Х0,Х1,Х2) = (х0 - Х2)(2Х0 -3x1 + 2Х2)/[б(Х1 - Х2)], (10.3)

где х0, XI, х2 - узлы интегрирования, в которых известна подынтегральная функция. При равноотстоящих узлах эти коэффициенты равны коэффициентам формулы Ньютона-Котеса второй степени с равномерным шагом.

Вычисляя крайние интегралы в (5) с помощью формулы Ньютона-Котеса второй степени с неравномерным шагом, а средний - по формуле (9), получим следующее решение интегрального уравнения (1):

" 1 ]

/4 + Sp(X1,71()+ X С2,кК3(х(,(2-к)))+к

Щ4 =

к=0

/[1-С2ДК3 (х(,0)], (11.1)

Щ5 =

К

3

1

/4 + Яр (X1,715) + 2 X К(Х5,X,)Щ + X СиКъ (,(2-к)) /[1 -С2,2К32,0)], (11.2)

2к=2 к=0 /

1

/г + Яр(X 1,71/) + Я(7,2,/-2,К)+ X С^Кэ(,(2-к)й))-2+к /1 -С2,2К3(,0)] , (11.3)

к=0

где Х1 71, - массивы из трех элементов; X1 = ^[]Ь , 71, = Кх(хг, Х;)Щ-;] = 0, 1, 2; коэффициенты С2, к

определяются один раз по формулам (10) для хк = ^/kй, к = 0, 1, 2. Функция Яр(х, /) определяется формулой

Яр (хУ) = X с2,г (х0,х1,х2)/г. (12)

г=0

Проведем сравнение результатов численного решения интегрального уравнения Хаффорда с помощью формул Ньютона-Котеса второй степени, (8) и (11) с результатами расчетов по формулам для однородной сферической Земли. Для вычисления функции ослабления на больших расстояниях (более 10 км) мы будем использовать формулу Фока [5], а для расстояний, меньших 10 км, будем использовать приближенное решение интегрального уравнения Хаффорда для малых расстояний (в сравнении с радиусом Земли) [12]:

(13)

где у(х) - функция ослабления для однородной плоской Земли; х - расстояние в метрах; Я =

Ях - численное расстояние.

Суша

10'

10

10

10

10

x, км

h = 1000 м — — h = 100 м ................... h = 10 м

комбинация ф. Фока и (13)

Рис. 1. Решение интегрального уравнения Хаффорда с помощью формулы Ньютона-Котеса второй степени (11) с шагом интегрирования h = 10, 100, 1000 м в сравнении формулами для однородной Земли

На рис. 1, а и б приведены результаты численного решения интегрального уравнения Хаффорда с помощью формулы Ньютона-Котеса второй степени над поверхностью суши для различных шагов интегрирования. Жирными точками на рисунке обозначены результаты расчета функции ослабления с использованием формул Фока [5] и (13). По результатам расчетов видно, что при плохой проводимости поверхности необходимо выбирать мелкий шаг (не более 100 м в нашем случае), иначе решение становится нестабильным и дает большую погрешность (шаг h = 1000 м), хотя с возрастанием расстояния решение может стабилизироваться. При высокой проводимости, которой обладает поверхность моря, гладкое и точное решение можно получить даже при шаге в 1000 м (рис. 1, в, г).

Оценить точность решений, которую дают различные квадратурные формулы, можно по расчётам, представленным на рис. 2. Численное решение проводилось для поверхности моря с шагом h = 1000 м. На рис. 2, а и б представлены модуль и фаза функции ослабления, найденные с помощью численного решения интегрального уравнения Хаффорда, а также вычисленные по формулам Фока [5] и (13) (жирные точки). На рис. 2, в и г приведены графики модулей относительных ошибок вычисления функции ослабления для абсолютных значений и фаз. Для определения относительной ошибки решения интегрального уравнения Хаффорда использовалась функция ослабления, рассчитанная по формуле Фока [5] и (13). По приведенным графикам видно, что использование формулы Ньютона-Котеса второй степени для среднего интеграла в (5) уменьшает погрешность численного решения приблизительно на один порядок, а замена формулы первой степени на формулу второй

2

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

---- Формула (10) — - Формула (11) ------Формула (14) • • • комбинация ф. Фока и (13)

Рис. 2. Численные решения интегрального уравнения Хаффорда при использовании различных квадратурных методов интегрирования с шагом h = 1000 м

Из приведенных выше расчетов можно сделать вывод, что для вычисления функции ослабления над слабопроводящими трассами необходимо разрабатывать численные методы с переменным шагом интегрирования. Их использование позволит сократить объем расчётов при допустимой величине ошибки. Кроме того, на точность полученного таким методом решения в сильной степени влияет величина погрешности вычисления крайних интегралов в (5). Повышая точность их вычисления, можно существенно повысить суммарную точность всего решения (рис. 2, в, г). Из рассмотренных здесь квадратурных методов численного решения интегрального уравнения Хаффорда наибольшую точность дает квадратурная формула Ньютона-Котеса второй степени. Тем не менее для расстояний свыше нескольких сотен километров даже для морских поверхностей она может давать большую погрешность, что делает такое решение интегрального уравнения непригодным в практических целях. Поэтому необходимо по возможности использовать квадратурные формулы более высокого порядка. Это можно легко сделать для среднего интеграла в формуле (11), если использовать квадратурную формулу Ньютона-Котеса третьей степени для всех точек, кроме четырех конечных, а формулу второй степени - только для одного элементарного интервала, когда общее число точек интегрирования является нечетным. Полученные результаты численного решения интегрального уравнения Хаффорда могут быть использованы для решения интегрального уравнения Фейнберга и других подобных уравнений [8, 10].

Литература

1. Кашпровский В.Е. Экспериментальное исследование распространения радиоволн / В.Е. Каш-провский. - М.: Наука, 1980. - 151 с.

2. Агафонников А.М. Фазовые радиогеодезические системы для морских исследований / А.М. Агафонников. - М.: Наука, 1979. - 164 с.

3. Нагуслаева И.Б. Исследование электрических свойств подстилающей среды и пространственно-временных характеристик электромагнитного поля по данным радиоизмерений и моделирования: автореф. дис. ... канд. физ.-мат. наук. - Иркутск: Иркут, гос. ун-т, 2010. - 26 с.

4. Саломатов Ю.П. Инженерные методы расчета распространения радиоволн вдоль поверхности Земли: учеб. пособие. - Красноярск: ИПК СФУ, 2009. - 164 с.

5. Фок В.А. Проблемы дифракции и распространения электромагнитных волн. - М.: Сов. радио, 2011. - 520 с.

6. Горшнев А.М. Решение интегрального уравнения для функции ослабления над сферической землей // Проблемы дифракции и распространения волн. - 1981. - Вып. 18. - С. 165-170.

7. Фейнберг Е.Л. Распространение радиоволн вдоль земной поверхности. - М.: Наука, 1999. -496 с.

8. Дембелов М.Г. Метод численного решения обобщенного интегрального уравнения Фейнбер-га для геометрически и электрически неоднородных трасс / М.Г Дембелов, Ю.Б. Башкуев // Физика волновых процессов и радиотехнические системы. - 2008. - Т. 11, № 1. - С. 89-94.

9. Распространение радиоволн над электрически и геометрически неоднородными трассами // Проблемы дифракции и распространения волн / Е.П. Проскурин, А.А. Пылаев, Н.П. Тихомиров, А.А. Штейнберг. - 1981. - Вып. 18.- С. 171-183.

10. Братцева А.В. Поле земной волны на высоте над электрически неоднородными трассами / А.В. Братцева, О.П. Буяло, Н.П. Тихомиров // Проблемы дифракции и распространения волн. -1987. - Вып. 21. - С. 179-190.

11. Hufford G. A. An Integral equation approach to the problem of wave propagation over on irregular suface // Quart. Appl. Math. - 1952. - Vol. 9. - P. 391-404.

12. Сенченко А.А. Решение интегрального уравнения Хаффорда при малых расстояниях от антенны / А.А. Сенченко, Ю.П Саломатов // Изв. вузов. Физика. - 2010. - №9/2. - С. 114-115.

13. Вержбицкий В. М. Численные методы (математический анализ и обыкновенные дифференциальные уравнения): учеб. пособие для вузов / В.М. Вержбицкий. - 2-е изд. испр. - М.: Изд. дом «ОНИКС 21 век», 2005. - 400 с.

Сенченко Александр Андреевич

Аспирант каф. радиотехники Сибирского федерального университета (СФУ)

Тел.: 8-902-941-91-67

Эл. почта: alsenchenko@mail.ru

Саломатов Юрий Петрович

Канд. техн. наук, доцент, зав. каф. радиотехники

Института инженерной физики и радиоэлектроники (ИИФиРЭ) СФУ

Тел.: 8-902-965-31-72

Эл. почта: YSalomatov@sfu-kras.ru

Senchenko A.A., Salomatov Yu.P.

A comparison of quadrature methods for a solution of Hufford integral equation

In this research paper there are given numerical solutions of Hufford integral equation for different Newton-Cotes formula. We calculated an attenuation function using numerical solutions of Hufford integral equation. The results were compared to the results obtained by Fock formula.

Keywords: numerical solution of integral equation, Hufford integral equations, Feinberg integral equations, Newton-Cotes formula, attenuation function, radio wave propagation.