Вычислительные технологии
Том 15, № 2, 2010
Регулярный алгоритм аппроксимации негладких решений для интегральных уравнений Фредгольма
первого рода*
В. В. Васин, Т. И. Серёжникова Институт математики и механики УрО РАН, Екатеринбург, Россия e-mail: [email protected], [email protected]
Построен и исследован регулярнзующнй алгоритм для восстановления негладкого решения уравнений Фредгольма первого рода. Алгоритм основан на использовании тихоновской регуляризации с недифференцируемым стабилизатором с привлечением проксимального метода и субградиентного процесса для решения задачи негладкой минимизации. Приводятся результаты численных экспериментов.
Ключевые слова: интегральное уравнение Фредгольма, негладкое решение, тихоновская регуляризация, проксимальный метод, субградиентный процесс.
Введение
Вычислительная практика решения интегральных уравнений первого рода убедительно показала, что применение в тихоновской регуляризации стабилизаторов в форме вариации различных типов позволяет значительно улучшить качество аппроксимации разрывных решений по сравнению с классической квадратичной регуляризацией, когда используется гильбертова норма L^h (см. [1-4]). При применении вариации функции в качестве стабилизатора удается обосновать кусочно-равномерную сходимость регуляризованных решений, что, по-видимому, и объясняет преимущество такого рода регуляризующих алгоритмов при восстановлении негладких решений.
Однако при численной реализации этих методов, т. е. при построении регуляризо-ванного семейства приближенных решений, приходится иметь дело с задачей негладкой минимизации. Здесь наметились два подхода: а) негладкий функционал предварительно аппроксимируется семейством дифференцируемых функционалов, а затем применяются традиционные методы гладкого анализа, например, методы градиентного или ньютоновского типа [1, 2]; б) для регуляризованной негладкой выпуклой задачи минимизации тихоновского функционала применяются субградиентные процессы [3, 4].
В работе [5] был предложен другой параметрический класс недифференцируемых стабилизирующих функционалов на основе нормы пространства Липшица и обоснована равномерная сходимость тихоновских аппроксимаций к непрерывному необязательно дифференцируемому решению исходного уравнения. Оказалось, что при подходящем выборе управляющих параметров регулярные алгоритмы, построенные в рамках упомянутой регуляризации с привлечением ргох-метода, вполне пригодны для г,осп а ног,-
* Работа выполнена при финансовой поддержке РФФИ (грант № 09-01-00053) и Междисциплинарного проекта УрО РАН. © ИВТ СО РАН, 2010.
ления негладких (разрывных) решений интегральных уравнений Фредгольма первого рода.
Заметим, что идея о возможности использования нормы пространства Липшица в качестве стабилизатора была высказана в обзорной статье [6].
В настоящей работе дается описание и обоснование всех этапов алгоритма и приводятся результаты численных экспериментов по решению интегрального уравнения, возникающего при продолжении геофизических полей,
1. Теоремы сходимости приближенных решений
Пусть П — компакт в Яп. Определим множество функций, удовлетворяющих на П условию Липшица
3 d> 0: |п(х^ — п(х2)| < d |я1 — х2|м V х\,х2 € П, (1)
/ п \ 1/2
где 0 < л < 1, |х1 — х2| = ^ = |х1г — , Обозначим через Нм = НМ[П] множе-
ство функций, удовлетворяющих на П условию (1) с фиксированным параметром По аналогии со случаем функции одной переменной определим норму
...... |п(х1) — п(х2) |
и(х) № = тах ф) + вир —:-:-,
Х6П Х1,Х2&П |х1 — х2|^
относительно которой Нм является банаховым пространством (см, [7]), Рассмотрим линейное уравнение
Ап = / (2)
с оператором А, действующим из пространства Нм в пространство непрерывных функций С(П), и единственным решением П € ННепрерывность обратного оператора А-1 не предполагается, поэтому (2) относится к классу существенно некорректных задач. При приближенно заданной правой части /, || / — /1| < 8, для построения регуляризо-ванного семейства приближенных решений используем метод Тихонова в виде [5, 7]
ш1п(||Ап — /||с[п] + а||п||ям : п € Н^[П]}. (3)
Теорема 1. Пусть А — линейный непрерывный оператор в С(П), Кег(А) = {0},
|| / — /б || < 8. Тогда задача имеет решение па, возможно неединственное, и при связи,
8
параметров а(8) —> 0, 0, 8 —> 0 имеет место равномерная сходим,ость
а(8)
б^о
кроме того,
б^о
Иш ||па(б) — П||с(п) = 0, Иш ||па(б)||Ям = ||П||Ям.
Замечание. Если вместо липшицевой нормы в (3) использовать стабилизатор
= |Н|Як + ^1Ы|£р, 1<р<оо, (4)
то решение па задачи (3) единственно.
Используя метод Рптца, можно построить аппроксимации регуляризованного приближенного решения ua элементами конечномерного пространства кусочно-линейных функций. Предположим, что компакт П допускает разбиение симплекеами n-ro порядка с максимальным диаметром h = maxAm. Обозначим через Uh конечномерное подпро-
m
странство кусочно-линейных (линейных на каждом симплексе) функций.
Теорема 2. Пусть выполнены предположения теоремы, 1. Тогда конечномерная задача
min{||Au - fs||с[п] + «ОЫЦрП + ||и^||ям] : uh G Uh} (5)
имеет единственное решение и^, для, которого выполнены соотношения,
lim ||< - иа||с[П] = 0, lim ||<||ям[П] = ||иа||ям[П],
h ^U h ^U
где иа — решение, полученное методом Тихонова со стабилизатором, (4)-
Как показали численные эксперименты, выполненные для интегральных уравнений Фредгольма первого рода с негладкими (разрывными) решениями, при малых значениях параметра а субградиентные методы работают неудовлетворительно, в результате чего не удается восстановить, например, хорошо выраженный разрыв или излом, С дру-
а
решения при данном методе регуляризации,
В значительной степени эти трудности удается преодолеть, если для устойчивой аппроксимации экстремального элемента иа использовать ргох-метод (итерированный вариант метода Тихонова), а именно: если обозначить целевой функционал в (3) через Фа(и), то итерационный процесс для аппроксимации иа строится следующим образом:
uk = argmin {Фа(и) + в ||и - uk-1||H : и G U}, (6)
где в > 0, || ■ ||h — некоторая гильбертова норма,
В отличие от выпуклого функционала Фа(и) функционал Фа'в = Фа(и) + в 11и — v||H является сильно выпуклым, что гарантирует устойчивый счет при нахождении uk субградиентным методом. По-видимому, это обстоятельство объясняет факт повышения качества аппроксимации негладкого решения, однако за счет увеличения трудоемкости алгоритма,
2. Алгоритм решения интегрального уравнения Фредгольма первого рода и численные эксперименты
Рассмотрим одномерное интегральное уравнение Фредгольма первого рода
1
Au = У K(x,y)u(x)dx = f (y), 0 < y < 1 (7)
U
K(x, y), f(y)
тально алгоритм, предложенный в предыдущем разделе. Прежде всего заметим, что конечномерная аппроксимация (см (5)), описанная в теореме 2, не ведет к полной дискретизации регуляризованной задачи, поэтому используем иную, более удобную для численной реализации схему [7].
Зададим равномерную сетку по х и по у с шагом к = 1 /М и аппроксимируем интеграл формулой правых прямоугольников.
Введем обозначения: п(х^) = п, /б(уг) = /г. Тогда дискретный аналог регуляризо-ванной задачи со стабилизатором (4) принимает вид
N
min < max ^ hK(yh Xj)uj — fi| +
+а
i<j<N j=1
1 N |u - u-I
X—^ 7 I I 2 I I i j
- > fiiUj + max Uj +max--—
2 j= jl i<j<N' jl i=j |x, - xjI1
{Uj }N e RN} . (8)
Теорема 3. Экстремальная задача, (8) имеет единственное решение Пм = (п1, п2,... и^), и последовательность кусочно-линейных восполнений пм (х), построенных по вектору пм, равномерно сходится к решению задачи,
-1,
min j || Au — f ||C + а
+ INI*"
: u e H1
Теперь дискретный аналог ргох-метода (6) принимает форму
12
uk ^^min {ФN(u) + ß ||u — uk-1||R : u e RN} = ФN k, (9)
где ФN — целевая функция в задаче (8),
Так как целевая функция в (9) является сильно выпуклой, то ргох-отображение
(см, [8]), определяющее итерационный процесс (9), принадлежит классу псевдосжима-
ющих отображений, и, следовательно, ргох-метод (9) сходится к решению задачи (8)
(см, [9, 10]), Поскольку выпуклая функция субднфференцнруема, то при фиксировап-uk-1 uk
Наиболее простой вариант субградиентного метода имеет вид [11]
Uk,u+1 = uk,u_ и = 0,1,2, ...,nk, (10)
где vk,V e ЗФ^к (uk,V), _ целевая функция в задаче (9), 5Ф означает субдифференциал функционала Ф. Как известно, в общем случае для выпуклой функции при
те те
Ak > 0, £Ak = оо, £ Ak < оо установлена сходимость итераций по функционалу
i=0 i=0
[11], но поскольку Ф^^ ~ сильно выпуклая функция, то имеет место сходимость по аргументу, т.е. lim ||uk,V — uk|| = 0.
Если известна некоторая оценка ФN, k для оптимального значения ФN k в задаче (9),
т, е,
ФN,k < ФN)k + £k, £k > 0,
то целесообразно перейти от задачи (9) к решению выпуклого неравенства
ФЙ (u) — ФN, k < 0. (11)
Теперь для приближенного решения неравенства (11) применим более эффективный
[12]
+ vk, v
uk>v+1 = - А (Ф%%(ик>») - Ф„,к) (12)
где vk'v £ (uk'v), 0 < Л < 2, v = 0,1, 2,..., nk, а величина Фможет уточняться
в процессе счета.
Важно отметить, что итерации uk ргох-метода (9) могут вычисляться субградиентными методами (10), (12) с некоторой погрешностью с сохранением сходимости, а именно, справедлива следующая теорема, обоснование которой следует из результатов [13, 14].
Теорема 4. Пусть Т — ргох-отображение, т. е.
T : v ^ arg min (Ф^(u) + в ||u - v||2 : u £ }. Пусть ргох-метод реализуется с погрешностью
ж
II/+1 - T(zk)||< 7fc, z0 £ , ]T7fc <
i=0
Тогда, lim ||zk — UN|| = 0, где UN — решение задачи, (8).
Замечание. Если известна априорная информация о решении уравнения (2) вида u £ Q, где Q — выпуклое замкнутое множество, заданное, например, системой линейных или выпуклых неравенств, то вместо ргох-метода можно использовать итерационный процесс [14]
uk+1 = Pq(T (uk)),
где T — ргох-отображение, Pq — фейеровекое отображение, отвечающее за ограничение u £ Q.
Таким образом, все описанные этапы алгоритма, построенного на основе тихоновской регуляризации со стабилизатором в форме нормы Липшица, ргох-метода, дискретной аппроксимации квадратурным методом и субградиентного процесса негладкой оптимизации, являются обоснованными, следовательно, мы имеем дело с регуляризую-щим алгоритмом,
В формулах (12) субдифференциалы каждого из слагаемых целевой функции вычисляются по известным правилам субдифференциального исчисления. Приведем, например, формулу для субдифференциала третьего слагаемого
ф(и) = тах т—--=7"- (13)
0<ij<N |x — x
Обозначим через J(u) множество пар индексов i0,j0, для которых в (13) достигается максимум. Определим вектор
11
-—, 0,..., 0, —--
/у . _ /у ^ /у . _ /у . ^
Кад ^Jo I КЗД Jo I
Vioh = 0,...,о,--—, о,...,о,---—, о,...,о
где ненулевые элементы занимают позиции с номерами го и тогда
М = с^„иеЛ«), если - > 0, = ^ со{Ми(—М)}, если — и70 = 0,
-М = Со{-г^}г0,^0€^(и), еСЛИ «¿0 — «70 < 0,
где co(M} — выпуклая оболочка множества M, u = (u0, u2,..., uN),
Для проведения численных экспериментов рассматривается решение интегрального уравнения Фредгольма первого рода [15]
которое моделирует ситуацию, когда необходимо найти продолженное на глубину х = Н гравитационное поле и(х) по измеренному на земной поверхности полю / (у) (плоский случай). Были проведены расчеты по восстановлению трех модельных решений. Использовалось нулевое начальное приближение. В формуле (8)для нормы Липшица полагалось ^ = 0.001, в (9) задавалось в = 300. Количество итераций по формулам (12) принималось равным 60. Для расчета вариантов с точной правой частью для задания параметра регуляризации использовались а = 10-4, 10-5, а для расчетов с возмущенной правой частью — а = 10-3, 10-4. Понятно, что с ростом величины параметра Н, 0 < Н < 2, процесс восстановления приближенного решения задачи (14) усложняется, поскольку свойства матрицы ухудшаются. Для моделей 1 и 2 полагалось Н = 2, для модели 3 — Н =1.
Для каждой модели рассчитывался вариант с точной правой частью и вариант с возмущенной правой частью, величина возмущения полагалась ||/ — /||с = 5 < 0.0108.
Точная правая часть для каждой модели рассчитывалась подстановкой точного решения и в сеточный вариант операторного уравнения Аи = /.
(14)
1 1 о о
0 20 40 60 80 100 120 х
а
1 1 0 о
0 20 40 60 80 100 120 х
б
0 20 40 60 80 100 120 *
в
Рис. 1. Модель 1. а — точное решение; б — восстановленное решение при точной правой части (Д1 = 1.7 х 10-2, Д2 = 1.0 х 10-6); в — восстановленное решение при возмущенной правой части (Д1 = 2.1 х 10-2, Д2 = 8.4 х 10-3)
Точное решение и(ж) (0 < ж, и < 2) для каждой модели рассчитывалось на равномерной сетке по ж, состоящей из 129 узлов с шагом Н = 2/128 ~ 0.016, по следующим формулам:
Модель 1 « =
Модель 2 « = <
Модель 3 « = <
1/32 х г,
г = 0,..., 64,
1/32 х (128 - г), г = 65,..., 128.
' 0, г = 0,..., 42, 2, г = 43,..., 85, 0, г = 86,..., 128.
0, г = 0, ... , 15,
1/8 х (г - 16), г = 16,..., 32,
2 - 1/32(г - 32), г = 33,..., 64,
, «128-г, г = 65, ..., 128.
В каждом варианте расчетов для полученного приближения вычислялись относительная погрешность по решению Д1 и относительная погрешность по невязке Д2:
Д1 =
- х\\ь2 \\х\\ь2
Д2 =
1|Ах - /||ь2
II/II
¿2
И -, 1.6" 1.20.80.40
20 40 60 80 100 120 х
1.61.20.80.4-
0 20 40 60 80 100 120 х
б
0 20 40 60 80 100 120 х
в
Рис. 2. Модель 2. а — точное решение; б — восстановленное решение при точной правой части (Д1 = 7.6 х 10-4, Д2 = 6.9 х 10-4); в — восстановленное решение при возмущенной правой части (Д1 = 3.9 х 10-2, Д2 = 2.7 х 10-2)
Рис. .3. Модель .3. а — точное решение; б — восстановленное решение при точной правой части (Д1 = 1.6 х 10-4, Д2 = 9.5 х 10-5); в — восстановленное решение при возмущенной правой части (Д1 = 7.1 х 10-3, Д2 = 5.1 х 10-3)
На рис. 1-3 приведены графики восстановленных решений как для точных, так и для возмущенных данных для трех модельных решений. На всех рисунках горизонтальная ось — номера узлов равномерной сетки по аргументу х.
Анализ численного моделирования показывает, что предложенный алгоритм достаточно хорошо восстанавливает решение с изломом (разрыв в производной) и разрывом первого рода самой функции, сохраняя структуру исходной модели. Несколько хуже алгоритм работает для более сложной модели 3. По-видимому, в данном случае для повышения качества решения необходимо использовать вместо и0 = 0 начальное приближение, учитывающее некоторые особенности решения, а также более тщательно выбирать управляющие параметры а,
Заключение
В работе предложен регулярный алгоритм, который позволяет восстанавливать негладкие решения интегральных уравнений первого рода. Алгоритм основан на тихоновской регуляризации с недифференцируемым стабилизатором (в форме нормы пространства Липшица) с привлечением ргох-метода и субградиентных процессов для решения задачи негладкой минимизации. Численные эксперименты подтверждают, что использование ргох-метода существенно повышает качество приближенного решения, поскольку дает возможность применять субградиентные процессы для сильно выпуклой функции, что обеспечивает устойчивость счета.
Список литературы
[1] Леонов A.C. Кусочно-равномерная регуляризация двумерных некорректных задач с разрывными решениями. Численный анализ // Журн. вычисл. математики и мат. физики. 1999. Т. 39, № 12. С. 1939-1944.
[2] Vogel C.R. Computational Methods for Inverse Problems. Philadelphia: SIAM, 2002. 183 p.
[3] Васин B.B., Серёжникова Т.И. Двухэтапный метод аппроксимации негладких решений и восстановление зашумленного изображения // Автоматика и телемеханика. 2004. № 2. С. 126-135.
[4] Vasin V.V., Korotkii M.A. Tikhonov regularization with non-differentiable stabilizing functionals // J. Inv. Ill-Posed Problems. 2007. Vol. 15, No. 8. P. 853-865.
[5] Васин B.B. Устойчивая аппроксимация негладких решений некорректно поставленных задач // Докл. РАН. 2005. Т. 402, № 5. С. 586-589.
[6] Тихонов А.Н., Васильев Ф.П. Методы решения некорректных экстремальных задач // Banach Center Publ. 1976. Vol. 3. P. 297-342.
[7] Васин B.B. Аппроксимация негладких решений линейных некорректных задач // Тр. Ин-та математики и механики УрО РАН. 2006. Т. 12, № 1. С. 64-77.
[8] Moreau J.J. Proximité et dualité dans un espace Hilbertien // Bull. Soc. Math. France. 1965. Vol. 93, No. 2. P. 273-299.
[9] Martinet B. Determination approachee d'un point fixe d'ne applications pseudo-contractante // C. R. Acad. Sei. Paris. 1972. Vol. 274. P. 163-165.
[10] Васин В.В., Агеев А.Л. Некорректные задачи с априорной информацией. Екатеринбург: УНФ "Наука", 1993.
[11] Поляк Б.Т. Введение в оптимизацию. М.: Наука, 1983.
[12] Васин В.В., Еремин И.И. Операторы и итерационные процессы фейеровского типа. Теория и приложения. Москва; Ижевск: НИЦ РХД, 2005.
[13] Rockape llar R.T. Monotone operators and the proximal point algorithm // SIAM J. Control and Optimizat. 1976. Vol. 14, No. 5. P. 871-898.
[14] Васин B.B. Проксимальный алгоритм с проектированием в задачах выпуклого программирования. Свердловск, 1982 (Препр. Ин-та математики и механики УНЦ АН СССР).
[15] Лаврентьев М.М., Романов В.Г., Шишатский С.П. Некорректные задачи математической физики и анализа. М.: Наука, 1980.
Поступила в редакцию 30 апреля 2009 г., с доработки — 20 ноября 2009 г.