CC BY
35
44 Вестник СамГУ — Естественнонаучная серия. 2006. №2(42).
ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА
УДК 519.999
ВЕЙВЛЕТ-АНАЛИЗ В ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТОНКОПРОВОЛОЧНЫХ АНТЕНН1
© 2006 А.С. Пименов2
Отыскание функции распределения тока и диаграммы направленности для системы тонких круглоцилиндрических проводников сводится к решению интегрального уравнения Фредгольма первого рода.
Для его решения применяется метод вейвлет-Галеркина на базе сплай-новых вейвлет. Использование псевдоразреженности матрицы системы линейных алгебраических уравнений (СЛАУ) позволяет применять алгоритмы технологии разреженных матриц. Изучаются свойства матриц и приводятся результаты численных расчетов.
Вычислительные трудности решения интегральных уравнений в задачах электродинамики во многом обусловлены заполненностью матриц, возникающих при их дискретизации. Это приводит к огромному объему вычислений (особенно в трехмерных задачах). Однако, как было показано в [1], в базисе из вейвлет-функций [2] большинство элементов матрицы СЛАУ оказываются очень малыми по абсолютной величине, т.е. она будет псевдоразреженной [3]. Использование этого факта позволяет существенно уменьшить объем вычислений и требования к оперативной памяти компьютера и открывает путь к созданию эффективных алгоритмов решения существенно трехмерных задач.
В настоящей работе рассматривается применение метода вейвлет-Галеркина к системе электрически тонких проводников.
1. Постановка задачи
Рассматривается система вертикальных электрически тонких круглоцилиндрических проводников (рис. 1).
1 Представлена доктором физико-математических наук профессором И.А. Блатовым
2Пименов Александр Сергеевич ([email protected]), кафедра прикладной математики Cамарского муниципального университета Наяновой, 443001, Россия, г. Самара, ул. Молодогвардейская, 196.
а(1)
Рис. 1. Система проводников
Задача состоит в отыскании функции распределения тока по проводникам при заданном возбуждении системы. Необходимо решить следующее уравнение Поклингтона:
Г , 1 д20(і, і') , ,
Е(1) = і[щі0С(1,Ґ)----------^-^сіГ]І(Г),
ше0 діді'
0(1, і) =
, Я(і, і') = V\Лі) - Кі')| + а2(і),
(1.1)
(1.2)
Я(і, і')
где Ь — контур, образованный совокупностью осей проводников; і, і — координата, отсчитываемая вдоль Ь; Е — заданная функция распределения стороннего поля ^-составляющая), возбуждающего систему; I — искомая функция распределения осевого тока; О — функция Грина; Г(і) — радиус-вектор точки на Ь; а(1)— радиус проводника в сечении /; |3 = — волновое число;
ш =
бл-108
^ —круговая частота; X — длина волны; |1о = 4л- 10 —магнитная
постоянная; ео = 8.854 ■ 10-7 — электрическая постоянная.
Один из проводников рассматривается в качестве активного вибратора. В среднем его сечении на контуре Ь выделяется короткий отрезок, на котором задается сторонний кусочно-синусоидальный ток (физически такой локальный источник соответствует зазору вибратора, к которому подключается фидер). Функция Е определяется по формуле:
Е(1) = -301(0(1, Ь) - 2ео8(|12 - Н\/2)0(11, (72 + к)/2) + 0(1, к)). (1.3)
С целью удобства физической интерпретации при тестировании в качестве окончательной характеристики предлагается использовать диаграмму
направленности (ДН) — функцию угловых сферических координат, описывающую распределение поля излучения на сфере весьма большого радиуса (много большего размеров исследуемой системы).
В данном случае ненормированная комплекснозначная ДН определяется по формуле:
В конечном итоге практический интерес представляет амплитудная ДН Ат(Г).
Для решения уравнения (1) применяется метод вейвлет-Галеркина. Вначале в п. 2 строятся сплайновые вейвлеты. Путем скалярного умножения левой и правой частей уравнения на эти функции мы получаем СЛАУ относительно неизвестных коэффициентов. Формулируется теорема о том, что матрица СЛАУ является псевдоразреженной, т.е. большинство ее элементов по модулю меньше заданного барьера.
Теперь для решения поставленной задачи применяем метод Галеркина, сначала приводя его к уравнению Поклингтона и проводя оценку ядра. Формулируется и доказывается теорема о сходимости в зависимости от малого параметра а. Если п — число точек разбиения контура, то \\иа - м„|| —> 0 при п — то, а — 0 : ап ^ С1, где С1 —некоторая константа.
Приводятся результаты численных экспериментов. Первые правдоподобные результаты получаем за 6 секунд. Решение этой задачи методом кол-локации дает те же результаты за 12 минут.
2. Сплайновые вейвлеты
Пусть [а,Ь] = [0,1], т — натуральное число и П0 — такое целое число, что 2П0 < 2т — 1 < 2П0+1. Рассмотрим семейство Д = (Дп, п = щ, П0 + 1,...} разбиений отрезка [0,1] Дп : 0 = х0 < хЦ <... < х^п = 1 с постоянным шагом к = Нп = 1/2п. На каждом из разбиений Дп рассмотрим пространство сплайнов Ьп = 5 (Дп, т - 1,1). Тогда для каждого к ^ п0 пространство S (Дк, т — 1,1) можно представить в виде прямой суммы Ьк = Ьп0 0 Шщ+\ 0 Жп„+2 0 ...0 Шк, где через Шк обозначено ортогональное дополнение пространства Ьк-\ до пространства Ьк. Искомый wavelet-базис будем строить как объединение базиса в Ьп0 и всех базисов в пространствах Шп, щ < п ^ к.
Вначале построим базис в ортогональном дополнении Шп пространства Ьп-1 до пространства Ьп. Зафиксируем п ^ п0 + 1. В случае необходимости
где V—орт направления излучения.
Нормировка ДН выполняется к своему максимуму:
(1.4)
F (0, ф) = /(0, ф)/ тах Ат(/).
(1.5)
будем считать, что каждое из разбиений Дп продолжено с тем же шагом на всю числовую ось узлами xn, —т < i < +ю. Нормализованные В-сплайны на разбиении Дп будем обозначать Nm-1,j,n.
Зафиксируем некоторое целое i ^ 0 такое, что i + 2m — 1 ^ 2n-1, то есть отрезок [xn_1, xn+21m_ 1] целиком содержится в [0,1]. Будем искать функцию ¥i,n(x) е Wn в виде
2i+3m—2
Vi,n(x) ^ ^ ajNm—1,j,n. I2.1)
j=2i
Для того чтобы Yi,n е Wn, достаточно потребовать выполнения условий
(V;,n> Nm-1,k,n-1) = 0, k = i — m + 1, i — m + 2,..., i + 2m — 2, (2.2)
поскольку остальные условия ортогональности выполняются автоматически в силу дизъюнктности носителей.
Подставляя представление (2.1) в (2.2), получим однородную систему 3m — 2 уравнений с 3m — 1 неизвестными 2i+3m—2
уЛ1 a j(Nm—1, j,n, Nm— 1,k,n—1) = 0, k = i — m + 1, i — m + 2,..., i + 2m — 2, (2.3)
j=2i
которая всегда имеет нетривиальное решение. Находя это нетривиальное решение, получаем искомый набор коэффициентов и функцию y;,n(x) в виде (2.1).
Из представления (2.1) и свойств В-сплайнов вытекает, что носитель sup(^i,n) с [x^., x^i+4 2], т.е. содержит 4m — 2 смежных частичных отрезка.
Возникает вопрос, нельзя ли построить wavelet с меньшей длиной носителя? Следующая теорема показывает, что это не так.
Теорема 2.1. Пусть p < m — 1, и функция Y;,n(x) вида
2i+m+2p
¥i,n(x) = ^ ^ ajNm—1,j,n (2.4)
j=2i
удовлетворяет условиям
(¥i,n, Nm— 1,k,n— 1) = 0, k = i — m + 1, i — m + 2,...,i + m + p — 1. (2.5)
Тогда Yi,n(x) = 0.
Следствие 2.1. Совокупность функций (Yi,n), i = 0,1,...,2n—1 — 2m + 1
линейно независима на каждом отрезке [xn_1, xrn^—jn_ 1], 0 ^ j ^ 2n—1 (т.е. линейно независимы их сужения на каждый такой отрезок).
Замечание 2.1. Справедлива формула у.,п(x) = ^0,П(2П—n0x — i/2n0—1), т.е. совокупность построенных wavelet-функций получается сдвигом одной единственной функции Y0,n(x).
Таким образом, мы построили совокупность полуортогональных линейно независимых wavelet {^,>}, i = 0,1,...,2n—1 — 2m + 1. Однако размерность ортогонального дополнения Wn равна 2n—1, т.е. до базиса в Wn нам не хватает ровно 2(m — 1) функций. Построим недостающие wavelet-функции. Для
этого рассмотрим функции yi,n(x) при -2m + 2 ^ i ^ 2n-1 -1 на расширенном разбиении Дп.
Первую группу из m - 1 недостающих wavelet будем искать в виде
—m
¥i,n(x) = Vi,n(x) - ^ aj-yjn(x), -m + 1 ^ i ^ -1, (2.6)
j=—2m+2
из условий
(Vi,n> Nm-1,k,n-1) = 0, к = -m + 1, -m + 2,...,-1, (2.7)
где скалярное произведение понимается в смысле L2[0,1].
Подставляя (2.6) в (2.7), получим систему линейных алгебраических уравнений
- m
/ i aj(yj,n, Nm- 1,k,n-1) = (yi,n, Nm- 1,k,n-1), /~
j=-2m+2 (2.8)
к = -m + 1, -m + 2,..., -1
для определения aj. Матрица системы (2.8) не вырождена, так как в противном случае существовало бы нетривиальное решение соответствующей однородной системы, что означало бы, что функция £-=-2m+2 aj¥j,n(x) является wavelet-функцией на [0,1] с носителем [x^ x2 2], что невозможно
в силу теоремы 2.1. Решая систему (2.8), получаем, что функция (2.6) является искомой wavelet-функцией, так как ортогональность к В-сплайнам Nm-1,k,n-1 при к ^ 0 имеет место в силу ортогональности им всех wavelet из линейной комбинации (2.6), а при -m + 1 ^ к ^ -1 — в силу условий (2.7).
Тем самым мы построили совокупность m-1 wavelet-функций (2.6). Их линейная независимость с ранее построенными функциями вытекает из вида (2.6) и следствия 2.1.
Вторую группу из m - 1 недостающих wavelet будем искать в виде
2n-1-1
¥i,n(x) = yi,n(x) - ^ ajyj,n(x), 2n-1 - 2m + 2 ^ i ^ 2n-1 - m, (2.9)
j=2n-1-m+1
из условий
(Vi,n, Nm-1,k,n-1) = 0, к = 2n-1 - m + 1,2n-1 - m + 2,...,2n-1 - 1, (2.10)
где скалярное произведение понимается в смысле L2[0,1].
Подставляя (2.9) в (2.10), получим систему линейных алгебраических уравнений
2n-1 -1
^ ^ aj(yj,n(x), Nm- 1,k,n-1) = (yi,n, Nm- 1,k,n-1), (2.11)
j=2n-1-m+1
к = 2n-1 - m + 1,2n-1 - m + 2,..., 2n-1 - 1
для определения aj. Решая систему (2.11), получаем, что функция (2.9) является искомой wavelet-функцией. Тем самым мы построили совокупность
т—1 wavelet-функций (2.9). Вместе с функциями (2.1) и (2.6) они образуют
искомый базис в ^п^ если 2п—1 > 2т — 1.
В качестве базиса в Ьщ выберем совокупность “усеченных” В-сплайнов
{ф* = iVm-u.no> -т + 1 < к < 2й0-1}. (2.12)
Итак, совокупность функций (2.12) и (2.1), (2.6), (2.9) при по +1 ^ п ^ к образует искомый wavelet-базис в пространстве Ьк.
Замечание 2.2. Из рассмотренных ранее аппроксимационных свойств сплайнов вытекает, что совокупность функций
{Ит-и,п0 |Гт+1}и
{{Vi.nl 2П01-2т+1} и {Vi.nl -=1—т— 1} и {Vi.nl 2=П2П-1т-2т+2}}}
^•П=Пг\ + 1 ^
^п=щ + 1
образует базис в Ь2[0,1].
Замечание 2.3. Для функций (2.6) и (2.9) при i е [2п—1 — 2т + 2,2п—1 — т] в силу симметрии справедливы формулы
¥;>(х) = ¥по + 1,2п—1—2т—1+1(2п—Щ —1(1 — х))
3. Алгоритм построения ’^вув^-базисов
Из полученных в предыдущем пункте результатов вытекает следующий алгоритм построения wavelet-базисов в S (Дк, т — 1,1) в случае произвольных т ^ 1, к : 2к > 2т — 1.
1) Для п = по + 1 численно решаем СЛАУ (2.3) при i = 0 следующим образом:
а) для ] = 0,1,...,3т — 2 вычисляем скалярные произведения
(Ит— 1.у,п, Ит— 1.к.п— 1), применяя для вычисления интегралов на каждом частичном отрезке, входящем в носитель В-сплайна, квадратурную формулу Гаусса, точную для многочленов степени 2(т — 1) [5].
б) полагаем азт—2 = 1 и находим аj, 0 ^ ^ 3т — 3, решая систему
с квадратной симметричной матрицей одним из прямых методов (например, методом квадратного корня [5]).
2) Для п = п0 + 1 аналогичным образом решаем системы (2.6) при — т +
+ 1 ^ ^ —1.
3) Определяем для п е [п0 + 1, к], i е [0,2п—1 — 2т + 1] —функции
¥;,п(х) = 2(п—п0—1)/2¥п0+1,0(2п—п0—1 х — г'/2п°), для п е [п0 + 1, к], i е [—т + 1, —1] — функции
Щп(х) = 2(п—п0—1)/2^п0+1,- (2п—п0—1 х),
а для п е [п0 + 1,к], i е [2п 1 — 2т + 2,2п 1 — т] —функции ^,п(х) = 2(п—п0—1)/2Уп0+1,2п—1—2т—;+1(2п—'10—1(1 — х)).
4) Присоединяем к построенным функциям функции (фк,п0, —т+1 ^ к ^ ^ 2 ° — 1}; где ф~ Nт—1,к,Щ'
Замечание 3.1. Таким образом, wavelet-базис в S (Дк, т — 1,1) получается сжатиями и симметричными преобразованиями т — 1-й функции ^о+и при i = —т +1, — т + 2,...,—1, а также сжатиями и сдвигами функции упо+1,о.
Замечание 3.2. Для вычисления функций из замечания 3.1 необходимо уметь вычислять значения В-сплайнов Nm—1,у,п—1. Их вычисление осуществляется с помощью рекуррентных формул (3.1), которые сводят вычисление Ит—1']п—1 к вычислению В-сплайнов первой степени по формулам
^1, j.n— 1 = ‘
х — хп_1, х е [хп—1, хп-1],
^ L j j + 1^
хп—2 — х, х е [хп—1, хп—2],
0, х £ [хп—1, хп—2].
^т- 1,*(0 = 1 Хк_ + Хк+т_ 1 А^ш-2Д+1(0. (3-1)
хк+т 1 хк хк+т хк+1
Замечание 3.3. Построение wavelet-функций в данном разделе проводилось в случае [а, Ь] = [0,1]. Однако очевидно, что на произвольном отрезке [а, Ь] соответствующую wavelet-систему можно получить из построенной с помощью линейной замены х' = , рассматривая \¥ауе1е1;-функции пере-
менной х'. Поэтому в дальнейшем мы будем считать, что wavelet-система построена для произвольного отрезка [а,Ь].
Замечание 3.4. В последующих разделах для краткости обозначений и приведения их в соответствие с общепринятыми в теории wavelet-функций положим
N т— 1, у,п(х) = Ф;,п(х),
считая т фиксированным.
4. Интегральные уравнения и метод Бубнова—Галеркина
Определение 4.1. Определенную на квадрате [а, Ь] X [а, Ь] функцию К(х, у) будем называть асимптотически т-гладкой, если она т раз непрерывно дифференцируема, и для всех ее производных 1-го порядка (1 ^ I ^ т) и некоторой константы а > 0 справедливы оценки
д К (х, у)
д хяду1'
1
I х —у|
а—I'
(4.1)
я
Пусть К(х, у) — асимптотически т-гладкая функция. Рассмотрим интегральное уравнение Фредгольма второго рода
и(х) + Г К(х, у)и(у) йу = /(х) (4.2)
^а
с заданной функцией / и неизвестной функцией и. Сложности численного решения уравнений такого вида (особенно многомерных) традиционными численными методами связаны с тем, что матрицы, получающиеся при их дискретизации, оказываются заполненными, т.е. состоящими из ненулевых элементов.
Рассмотрим для (4.2) метод Бубнова-Галеркина на базе построенных wavelet-функций степени т - 1. Зафиксируем некоторое натуральное к > щ и будем искать решение (4.2) в виде
2по -1 к-п0 2п0+;-1-т
и = й0уфу,по + с1]^у,по+1 (4.3)
у=-т+1 1=1 у=-т+1
из условий
ПЬ . 2п0 -1 к-по 2п0+;-1-т ^
К(х, у)( ^ ^ й0уф-/,ио(у) + С-у^у.по+г(у)) ф/,по(х) йх =
/=-т+1 1=1 /=-т+1
j=-m+1 i=1 j=-
Г f (х)ф/,и0(х) dx, -m + 1 ^ l ^ 2n0 - 1,
da
Uo
,ъ пъ , 2П0 -1 k—no 2n0+l-1-m
ПЪ / ^ " * Ж 0 Ж_ \
K(x, y)( ^ ^ d0j^j,n0(^) + z z cij¥ У,и0+1-Сул ¥/,n(x) dx
j=-m+1 1=1 j=-m+1
Г ъ
= I /У,Лч,г /,А ,J,„ . 1 7 on-1
a
f(x)¥/,n0(x) dx, -m + 1 ^ l ^ 2n - m, щ + 1 ^ n ^ k. (4.4)
Совокупность условий (4.4) представляет собой систему линейных алгебраических уравнений с квадратной матрицей порядка 2k+m-1, элементы которой имеют вид
ХЪ рЪ г*Ъ
ф ^j,n0 (x^/,n0 (x)dx + I K (x, _у)ф j,n0 (у)ф/,п0 (x)dxdy, ъ ^ъ Ja Jarъ rъ
K(x, у)фj,n0(y)¥/,n(x)dxdy, I I K(x, у)ф(y)¥/,n(x)dxdy,
ъ rb ja ja (4.5)
K(x, y)yj,n(y)¥/,n+s (x)dxdy, s * 0,
nh rb nh
I ¥ j,n (x)V/,n (x)dx + I I K (x, y)y j,n (y)¥/,n( x)dxdy.
*Ja da da
Для классических систем функций в методе Галеркина числа (4.5) оказываются, в основном, ненулевыми и недостаточно малыми, чтобы ими можно было пренебречь и рассматривать СЛАУ (4.4) как разреженную. Получим оценки чисел (4.4) в случае wavelet. В целях упрощения технических деталей рассмотрим случай а = 0 в (4.1).
Теорема 4.1. Найдется такая константа С > 0, не зависящая от к, что при —т + 1 ^ І ^ 2п—1 - т, —т + 1 ^ ^ 2п—1 - т, по + 1 ^ п ^ к справедливы
оценки
Пь
‘
K(x, y)yj,n(y)^l,n(x) dx dy n+s-1, 1 < s ^ к — n
< C2-n(1 + jl - jj)-
а при -m + 1 ^ j ^ 2'
-b rb
K(x, y)yj,n+s (y)rn,n(x) dx dy
nb
< C2-(n+s(m+1))(1 + jl - j/2sj)-
(4.6)
(4.7)
Из теоремы 4.1 следует, что матрица системы (4.4) является псевдораз-реженной, т.е. в ней очень много малых по модулю элементов. Учитывая, что СЛАУ для интегральных уравнений Фредгольма второго рода хорошо обусловлены, пренебрегая этими малыми элементами, мы получим хорошую разреженную аппроксимацию этой матрицы, для которой можно применять алгоритмы для разреженных систем [6].
5. Метод решения поставленной задачи
Будем решать уравнение (1.1).
Представим его в следующем виде:
Е(1) = Г К(I, Г)1(Г)М',
JK(l, l')I(l')dl', (5.1)
1 d2G(l, l')
где K(l,l) = I) ---------------) — ядро интегрального уравнения
шє0 dldl’
Фредгольма первого рода (5.1). Функции E(l),K(l, l') имеют разрывы первого рода в точках l = co,l = сі,..., в которых контур L прерывается (т.е. на концах проводников). Поскольку токовая функция на концах проводников должна иметь нули (из физических соображений), она является непрерывной, но негладкой в этих точках.
Проведем оценки ядра для рассматриваемых систем (от l линейно зависит только вертикальная координата радиус-вектора, радиус проводников a(l) = const).
Лемма. Для ядра уравнения (5.1) верна оценка
C
<5‘2>
Доказательство. Пусть z(l) — вертикальная координата, соответствующая положению точки l на контуре L.
д2С{1,Г) _т,\ ф o2{z{T)-z{V)f o.d(z(D-z(0)\
(¥ + ¥2~3 & +р ЇР зф ¥ );
Сначала докажем, что jr(l) - Г(і')| ^ Cjl - l'j. Если l, l' лежат в пределах одного проводника, то |?(Г) — ?(lr)\ = \I — Ґ|. Иначе \?(1) — ?(1Г)\ =
= л/Ш ~ z(V))2 + (х(0 - х(1'))2 ^ \х(1) - х(Г)| ^ C\l - V|, где С =
тт(|х(/) - х(1')|) — минимальное расстояние между двумя проводниками контура.
Так как
= А < А ш-т1_ г_
я2 я3 ^ я3’ я5 ^ /г5 ^ я3’
то учитывая, что |г(0 - z(l')| ^ Ь, Я = - ?(1')\2 + а2, получаем:
1 дЮ1' 1 ^ С|/гз|+|/г2| + :51 \+Р^\к3\+5р\ к4
<--------—--------•
«7 - V)2 + а2)3!2 ’
е-г|Щ/,0 1 С2
\0(1,01 = I
Я(ї, ї') ' ^ |Я| ^ ((ї - І'}2 + а2)3/2’
1 д2С(ї, І') С
\Щ,Ґ)\ < |СО[10С(/, 1)\ + І-------^| < ——2----------------5^2
шє0 дїдї' ((І - І')2 + а2)3/2
Лемма доказана.
Применим метод Бубнова-Галеркина на базе вейвлет-функций степени т - 1 для решения системы (5.1). Зафиксируем натуральное к > по и ток будем искать в виде
2п0 -1 к-по 2п0+і-1-т
1 = ^]<$]т + сі№ ]'по+і (5.3)
]=-т+1 і=1 ]=-т+1
из условий
/у=о К(х ,У)!(У)аУ, ф/, »о(х)) = (£(х) , ф/, «о(х)) ,-т + 1 ^ (54)
^ ^ 2П0 - 1,
(/У=0 К(х, у)/(у)4у, у/,и(х) | = (£(х), ^/,п(х)), -т + 1 ^ (55)
^ I ^ 2п0 - 1, по + 1 ^ и ^ к.
Здесь ф/,п0(х) —сплайны и ^/,п(х) — сплайновые вейвлеты третьей степени [2] .
Совокупность условий (5 . 4)—(5. 5) представляет собой систему линейных алгебраических уравнений относительно коэффициентов dоj, с квадрат-
ной псевдоразреженной матрицей порядка 2к + т -1, элементы которой имеют вид:
К (х, у)ф ^щ(у)ф1,п0( x)dxdy, (5.6)
их=0 J у=0
К (х, у)ф 1щ(у)у1,п(. x)dxdy, (5.7)
их=0 J у=0
Г Г К (х, у)у j,n(y)yl,n(x)dxdy. (5.8)
их=0 ^ у=0
Для решения системы необходимо вычислить элементы ее матрицы. Матрица системы состоит из сплайнового, вейвлетно-сплайновых и вейвлетных блоков.
и применением формулы Гаусса с р узлами на каждом отрезке.
Согласно [3], большинство элементов матрицы пренебрежимо малы по модулю. Поэтому заполнять матрицу мы будем следующим образом. Элемент матрицы будем считать нулевым, если он по модулю меньше заданного барьера 0 < е ^ 1. Такие элементы мы хранить просто не будем.
Матрицу мы вычисляем целиком, применяя быстрое вейвлет-преобразование (FWT). Идея его состоит в том, что все сплайны и вейвлеты уровня п - 1 представляются в виде линейной комбинации сплайнов предыдущего уровня (п). Поэтому мы можем вычислять интегралы по самому мелкому разбиению со сплайнами самого верхнего уровня к. Затем интегралы с вейвлетами уровня к мы вычисляем как линейную комбинацию вычисленных интегралов и таким же образом (со сплайновыми коэффициентами) находим интегралы со сплайнами уровня к -1. Таким образом, мы получаем быстрый алгоритм для вычисления всех элементов матрицы (с учетом, конечно, что интегралы на уровне к уже посчитаны).
Естественно, надо учитывать, что алгоритм FWT весьма прост в одномерном случае, а в двумерном использование его напрямую приведет к большим затратам памяти (больше, чем уйдет на хранение самой матрицы). Поэтому в каждой точке разбиения отрезка [0; Ь] х, мы считаем все г ь
интегралы ^=0 К(х,, у)у(у^у(*), где у обозначает все вейвлеты уровней п0 + + 1..к и сплайны уровня п0 указанным выше способом. При этом в памяти храним последние незаконченные (т.е. недосчитанные) т сплайнов каждого уровня п0..к и 3т - 2 вейвлет каждого уровня п0 + 1..к и прибавляем к сплайнам посчитанные интегралы (*), а досчитанные сплайны добавляем к вейвлетам и сплайнам следующих уровней и т.д. Законченные вейвлет-ные и сплайновые элементы уровня п0 мы добавляем в матрицу, если они больше заданного барьера.
После заполнения матрицы мы вычисляем вектор правой части:
и с помощью метода сопряженных градиентов находим решение искомой СЛАУ. Матрица этой СЛАУ А не является эрмитовой, поэтому для сходимости метода необходимо домножить обе части на комплексно-сопряженную к матрице А матрицу А*. Реализуя этот вариант, нет никакой необходимости перемножать А и А*, т.к. это потребует определенного времени и структура матрицы А нарушится, вместо этого в методе сопряженных градиентов на каждой итерации на вектор решения умножается сначала А, а потом А*.
Интегралы будем вычислять разбиением всего контура на 2к отрезков
(5.9)
6. Обоснование метода
Обоснование метода основывается на следующей теореме об операторных уравнениях.
Рассмотрим операторные уравнения
КаЫ = /, РпКаЫп = Рп/, (6.1)
где Ка — семейство операторов, зависящих от параметра а е (0, а0], действующих в банаховом пространстве Е с нормой || ■ ||, Рп : Е ^ Еп —операторы проектирования на подпространства Еп с Е.
То, что второе уравнение (6.1) следует из первого, можно показать следующим образом:
(Каип, ^п) = (/, Уп);(Каип /, Уп) = 0 ^
^ Каип - I ^ Еп ^ Рп(Каип - /) = °.
Пусть РпК — сужение оператора РпК на Еп.
Теорема. Пусть удовлетворены следующие условия.
1. Первое из уравнений (6.1) при всех а е (0, а0] имеет единственное
решение иа е Е, и может быть найдена константа С1 > 0, для которой
|| иа - Рпиа ||^ 0 при п ^ю, а ^ 0 : ап ^ С1.
2. Оператор РпК обратим, и для всех п ^ п0, а е (0, а0] : па ^ С1 верны
оценки II (РпКа)-1 РпКа |К С2.
Тогда | иа - ип ||^ 0 при п ^ю, а ^ 0 : ап ^ С1 и верны оценки || ип -
Рпиа 11^ С2 У иа Рпиа ||.
Доказательство.
(/=КаЫа )
РпКпип = Рп/ РпКа(ип - Рпиа) = РпКаиа - РпКаРпиа;
РпКа(ип - Рпиа) = РпКа(иа - Рпиа);
ип - Рпиа = (Рп-Ка) РпКа(иа - Рпиа) ||ип - Рпиа\\ ^ С2||иа - Рпиа||;
||иа - ипУ ^ ||иа - Рпиа\\ + ||un - Рпиа\\ ^ (1 + С2) У иа - Рпиа У * °.
Теорема доказана.
Выполнение условий теоремы для рассматриваемого метода Буб-нова-Галеркина в случае Е = Ьто[/], Еп = S(Д, 3,1) следует из аппроксима-ционных свойств сплайновых вейвлет, ограниченности последовательности операторных норм У Рп ||то [6], априорных предположений о гладкости функции тока, представления оператора задачи в виде РпКа = РпКа,1 + РпКа,2, где
К(х, у), |х - у| ^ Сза
Ка,1 —интегральный оператор с ядром К\(х,у) = . 0 \ \ > с
где С3 достаточно большая, но не зависящая от а, п константа. В самом
деле, из вида К(х, у) следует обратимость РпКау.
РпКа = РпКа,1 + РпКа,2;
\К(х,у)\ ^ ------------5-, \х-у\^ С3а,
((х - у)2 + а2)) 2
Ка, 1 ип = К1(х, у)ип (у)йу =
^0
Jr*min(L,x+Cз а)
К1(х, у)ип(у)йу = ип(у')/(х) ^ ип(х)/(х), а ^ 0
тах(0,х-Сз а)
т.е. имеет место сильная сходимость. Т.к. сфера в конечномерном пространстве является компактным множеством, то следует равномерная сходимость на этой сфере, и т^е£п,|МК1 \\РпКаЛи\\ ^ С||/(х)|| ^ Значит, существует обратный оператор ||(РпКад)-1|| ^ а2.
Тогда, т.к. а — малый параметр, имеем следующее.
(рпЖа)-^=(Р^Л+Р^2Т1_= (адГ1(/ + (РЖ,!)-РпКО^))-1 =
= (I + (^Го-1 2)-1( )-1;
\\Ка,2(х,у)\\ < -£-;=> ||М|| = < Со < 1;
аСз
1
||(/+МГ1|| = ||/-М + М2-М3 + ...|| < 1 + ||М|| + ||М||2 + ... = у-щ =>
^ КЛ^)-^ < Са2 ^У (ККаТ1 РпКа |К С
7. Антенна Уда—Яги
Вычислительный эксперимент проводился на системе Уда-Яги с пятью проводниками одинакого радиуса (рис. 2).
Длина волны 4 метра. Вертикальные размеры проводников (слева направо). 2; 2; 1,8; 1,72; 1,64 м.
Абсциссы проводников. 0,4; 1; 1,8; 2,8; 3,7 м.
Радиусы проводников 0,01 м.
Второй проводник слева — активный вибратор с зазором 0,04 м.
Вся система обладает зеркальной симметрией относительно плоскости г = 0.
Было проведено несколько экспериментов, нацеленных на изучение точности решения уравнения (5.1) при различных к. Точность решения определяется двумя параметрами. 1) точность решения интегрального уравнения как сравнение Е(1) с левой частью уравнения (5.1) при подстановке в нее полученного тока 1(1); 2) различия в диаграммах направленности при различных к.
Было выяснено, что для обеспечения приемлемого вида решения и диаграммы направленности необходимо, чтобы интегрирование производилось по разбиению из более чем 3000 точек. Меньшее количество точек не обеспечивает достаточную точность интегрирования ядра, что приводит к слиш-
1 -
0.5
2
3
-0.5
-1 [
Рис. 2. Антенна Уда-Яги
ком большим погрешностям. Для малых к (к < 9) добавляются дополнительные точки разбиения.
Похожая на правильную диаграмма направленности получается при минимальных к = 6, количество отрезков 29 = 512, количество точек в формуле Гаусса (р) 6. Время счета — 8 секунд, размер матрицы — 64 КЬ. Барьер (т.е. число, для которого любое число, меньшее его по модулю считается равным 0) выбран 0.001, количество ненулевых элементов — 4107 (из 4761 элементов, т.е. 86%). Графики ДН приведены на рисунке 3, 4.
При к = 11 абсолютная погрешность решения (как разность между Е(1) и левой частью уравнения (5.1)) достигает 0.05 процентов, но относительная погрешность — 365 процентов на концах проводников. При р = 4 время счета матрицы 1.5 минуты, размер матрицы 4 МЬ. Барьер тот же, количество нулевых элементов — 4017082 (из 4206601, т.е. 95%). На рис. 5-8 представлены графики тока и ДН.
Процент заполненности матрицы можно видеть в табл. 1.
Таблица 1
Процент заполненности матрицы при различных k
к = 6 к = 1 к= 8 к = 9 к = 10 к = 11
85.4% 57.9% 32.2% 16.8% 8.8% 5.1%
Сравнение времени выполнения при различных барьерах представлено в табл. 2.
Таблица 2
Сравнение времени выполнения при различных барьерах е
е к = 6 II оо II II 40 О II к = 11
0.001 7 сек 16 сек 39 сек 1 мин 35 сек 3 мин 52 сек 6 мин 35 сек
0 7 сек 21 сек 1 мин 13 сек 5 мин 40 сек 18 мин 80 мин
Рис. 3. ДН к = 6, 0=f
Рис. 4. ДН к = б, ф = О
Рис. 5. Ток, к = 11
Рис. 6. Пик тока на отрезке [2;4]
Рис. 7. ДН, к =11, 0=| Рис. 8. ДН, к =11, ф = О
Заключение
Таким образом, для уравнений Фредгольма первого рода специального вида предложен и строго математически обоснован метод вейвлет-Галерки-на. Написана программа и проведен численный эксперимент, продемонстрировавший преимущество этого метода перед применяемыми в настоящее время численными методами решения таких задач.
Литература
[1] Blatov, I.A. On estimates of the inverse matrices elements for the Galerkin method for singular integral equation based oh the spline wavelets / I.A. Blatov // Proceedings of the International conference on computational mathematics. Part II. Novosibirsk, 2002. P. 356-361.
[2] Добеши, И. Десять лекций по вейвлетам / И. Добеши. Ижевск: НИЦ ’’Регулярная и хаотическая динамика”, 2001. 464 с.
[3] Блатов, И.А. Об алгебрах операторов с псевдоразреженными матрицами и их приложениях / И.А. Блатов // Сибирский мат. журнал. 1996. Т. 37. №1. С. 36-59.
[4] Glinsky, B.M. Composition of wave fields of vibrators and calibration shots of the OMEGA series / B.M. Glinsky, V.V. Kovalevsky, A.S. Alekseev // Third international conference ’Monitoring of nuclear tests and their consequences”. - Borovoye, 2002. - P. 23-25.
[5] Бахвалов Н.С. Численные методы j Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. М.: Наука, 1987.
[6] Писсанецки, С. Технология разреженных матриц j С. Писсанецки. М.: Мир, 1988.
Поступила в редакцию 18jXIj2005; в окончательном варианте — 19jX//j2005.
APPLICATION OF SPLINE WAVELET-FUNCTIONS TO NUMERICAL MODELLING OF THIN-WIRE ANTENNAS3
© 2005 A.S. Pimenov4
Finding function of distribution of the current and the directional diagram for the system of the slender conductors is reduced to the solution of an integral Fredholm equation of the first type. The Galerkin method on the base of spline wavelets is applied to solve it. The algorithms of the technology of the sparse matrices are researched and the results of the numerical experiments are stated.
The computational difficulties of integral equations solving in the tasks of electrodynamics are caused in many respects by the fullness of the matrices that appear during their discretization. It results in an enormous size of computations (especially in 3-dimensional problems). However, as it is shown in [1], in the basis of wavelet functions [2] the most elements of the SLAE matrix turn out to be very small in absolute value, i.e. it will be pseudosparce [3]. That allows to reduce significally the amount of computations and the memory requirements and opens the way to the creation of the effective solution algorithms of essential 3-d problems.
In this work the application of the wavelet-Galerkin method to the system of electrically thin conductors is considered.
Paper received 18jXIj2005. Paper accepted 19jX//j2005.
3Communicated by Dr. Sci. (Phys. & Math.) Prof. I.A. Blatov
4Pimenov Alexander Sergeyevich ([email protected]), Dept. of Applied Mathematics, Samara Marina Nayanova University, Samara, 443001, Russia.
