Численный алгоритм прямого решения вариационной задачи А. Н. Тихонова Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ ЦА Г И Т о м VII 197 6

№ 4

УДК 533.6.011.5:532.582.33

ЧИСЛЕННЫЙ АЛГОРИТМ ПРЯМОГО РЕШЕНИЯ ВАРИАЦИОННОЙ

ЗАДАЧИ А. Н. ТИХОНОВА

В. Д. Перминов

Предложен новый эффективный численный алгоритм решения интегральных уравнений 1 рода методом регуляризации А. Н. Тихонова. Эффективность алгоритма проверена на примере решения уравнения Абеля.

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

ь (х)

А[х,г] = ( К (х, я) г (я) — / (х); а0 < а (х) < 5 < Ь (*) < Ь0. (1)

а (х)

Задача о решении такого уравнения поставлена некорректно. Если для заданной /(х) решение уравнения (1) единственно, то решение г(«) может быть получено методом регуляризации, предложенным А. Н. Тихоновым в работах [1], [2]. Согласно этому методу, приближения гл (в) определяются, как функции, доставляющие минимум функционалу

!’о *0

ма

[г, /] = | {А [х, г] - f }^Лх + а|[р (я) + Я (я) .г'2] (2)

где р ($) и <7 («) — заданные непрерывные положительные функции, а а — некоторый параметр. При этом, если правая часть / уравнения (1) задана со среднеквадратичной ошибкой 8, то при а = сб2 (с — любая постоянная величина, не зависящая от В)

вир | га (в) — г (в) | г (8, с) (е ^ О при 5 0),

Йо 5 Ь() *

где (б-) доставляет минимум Ма [г, /].

Как и при решении любой вариационной задачи, для решения (2) можно пользоваться либо прямыми, либо непрямыми методами вариационного исчисления. В известных автору работах (см., например, [3], где приведена обширная библиография) использовалась вторая группа методов, т. е. для функционала (2) выписывалось уравнение Эйлера, которое затем решалось на ЭЦВМ с помощью разностных методов. Этот метод, на наш взгляд, обладает двумя основными недостатками. Во-первых, программу, реализующую этот алгоритм, трудно сделать универсальной, так как точность используемых квадратурных формул зависит от пределов интегрирования а{х) и Ь (х). И, во-вторых, трудно обобщить этот алгоритм на двумерные интегральные уравнения.

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

Разобьем интервал [д0> &о] на N интервалов узлами в у. а0 = в0< «1<С • • • < %=6о. в которых зададим значения функции 2} = Тогда сплайном 5(«) называется

функция, непрерывная вместе со своими первой и второй производной, совпадающая с кубическим полиномом на каждом интервале 5/_1 С $ О; (у = 1, 2, ... ,Ы) и удовлетворяющая условиям

5(5;) = г;- О' = 0, 1...Ы).

Для полного определения сплайна необходимо задать еще два дополнительных краевых условия, в качестве которых можно взять значения на концах первой или второй производной, если эти значения известны, или, в более общем случае, некоторые связи между вторыми производными в двух начальных и двух концевых узлах сетки. При наличии этих условий задача определения коэффициентов кубических полиномов сводится к решению системы линейных алгебраических уравнений с трехдиагональной матрицей, которая легко и устойчиво решается методом прогонки. Свойства сходимости сплайнов и их производных можно найти в работе [4].

Определим теперь N + 1 фундаментальных сплайнов с помощью соотношений

А) (Я;) = 8,у (г, У, = 0, . .., Ы); ^

А) (во) = А) («О, А"/(зя_1)=А’{8к). / (3)

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

. N N

г (в) а 5 (в) = ^ А) («) г у, г' (5) = 5' (я) = 2 А\ (в) г].

1=0 /—О

Подставляя эти выражения для г (я) и г' (э) в (2), после некоторых элемен-

тарных преобразований получим

N N

М* [г, Л = 2 г* г1 + 2 Р/ г1 + сь (4)

к, } = 0 /=0

N

[//£ Л'у Н" аЧ’ь (^**)] ~\~ Р?г'> (5)

<=0

N N ■

?* = - 2 2 /,*, с, =2 в^\; (в)

г=0 /=0

Ь.

I

/«= \к(Х1^)Ак^)(И-, (7)

'а1

здесь Л/ — коэффициенты выбранной квадратурной формулы для внешнего интеграла в (2).

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

известными методами линейной алгебры.

Заметим, что использование фундаментальных сплайнов позволяет для

широкого класса ядер К (х, «) вычислять интегралы (7) в квадратурах. В том случае, когда этого сделать не удается, для вычисления этих интегралов может быть использована любая квадратурная формула. Кроме того, из (5) — (7) видно, что с,, р* и /,* вообще не зависят от параметра регуляризации о, а коэффициенты матрицы распадаются на отдельные слагаемые, каждое из которых либо не зависит от а, либо пропорционально ему. Это позволяет рассчитывать эти величины (или их слагаемые) только один раз, что, естественно, уменьшает общее машинное время.

Остановимся на некоторых возможностях и преимуществах предлагаемого алгоритма:

— алгоритм универсален в классе уравнений типа (1) в том смысле, что с его помощью могут быть решены интегральные уравнения первого рода, ядра которых имеют слабую особенность (типа уравнения Абеля, см., например, [5]);

— дополнительные условия, используемые для образования фундаментальных сплайнов, накладывают довольно слабые ограничения на искомую функцию. В частности, условия (3) означают лишь, что на двух крайних интервалах искомая функция аппроксимируется параболой, а не полиномом третьей степени, как в остальных интервалах. Если граничные условия для искомой функции известны точно, то этот факт легко учесть при построении системы фундаментальных сплайнов;

— так как сплайны имеют непрерывную вторую производную, то имеется возможность без принципиальных переделок увеличить порядок регуляризации п [2] [здесь п — порядок старшей производной во втором слагаемом в (2)];

— если воспользоваться двумерными фундаментальными сплайнами [4], то можно (по крайней мере, для прямоугольных областей определения искомой ■функции) построить алгоритм решения двумерного линейного интегрального уравнения первого рода. Задача сводится к решению системы № линейных алгебраических уравнений (здесь N—число узлов на каждой стороне прямоугольника). Это в свою очередь требует запоминания № коэффициентов матрицы. При не слишком больших N и использовании внешних запоминающих устройств такая задача может быть решена на БЭСМ-6.

Предложеннный алгоритм для одномерного случая был реализован в виде программы на автокоде БЕМШ для БЭСМ-6 и опробован на задаче решения интегрального уравнения Абеля, к решению которого сводится задача определения плотности в осесимметричных течениях методом сдвига интерференционных полос [5, 6]. Это уравнение имеет вид

1

к (Л ^ = Б (и). (8)

3 У V — и

и

где V =и = ; V (^)= ^ ^° ; к — ; здесь г и у — расстояние от

Л?2 К2 Ро А

оси симметрии; р — плотность, 5—сдвиг интерференционных полос, /? —радиус течения, б —постоянная Гладстона — Дейла для данного газа, р0 — некоторое известное значение плотности и X — длина волны источника света.

Сравнение эффективности предлагаемого алгоритма решения уравнения (8) с принятым в настоящее время методом сведения таких уравнений с помощью квадратурной формулы к системе линейных алгебраических уравнений [5,6] проводилось на модельной задаче. Непосредственной подстановкой в (8) легко

проверить, что решение vo = (l—г»)2 соответствует функции 50 (и)= к (1—и)5'2.

15

Ошибки, присутствующие при обработке экспериментальных данных, моделировались наложением на 50 случайного возмущения, т, е. вместо 5“0 задавалась при расчетах следующая функция:

5(и) = .9оИ-Нд. (9)

где к) — случайная величина, равномерно распределенная в интервале [ —1,1]; Д — некоторый постоянный множитель, характеризующий максимальную величину возмущения.

При решении модельной задачи считалось, что £=/; (в) = <?($)= 1, использовалась регуляризация первого и второго порядков (в последнем случаев сглаживающий функционал (2) входит еще и квадрат второй производной [2]) с определением оптимального значения параметра регуляризации а с помощью критериев квазинаилучшего приближения [3] и невязки [7].

Сформулированная задача решалась при N=16 для значений Д = 0; 0,01;

0,05. При Д = 0 правая часть задана точно (по числу разрядов в ячейке машины) и основную роль играют ошибки аппроксимации интегрального оператора в (8) — ошибки квадратурной формулы и ошибки замены искомой функции сплайнами, проведенными через значения функции в узлах. В этом случае ошибки восстановления ч (и) методом квадратур составляют 0,8%, а методом регуляризации первого порядка — около 0,1% от максимального значения. При А = 0,01 и 0,05 (фиг. 1) ошибки восстановления методом квадратур составляют соответственно 2,5—3 и 11—-12%, а методом регуляризации первого порядка 2 и 7—8%. Применение регуляризации второго порядка позволяет не только получить в качестве решения более гладкую кривую, но и несколько улучшить точность восстановления решения (1,5% для Д == 0,01 и 6% для Д = 0,05).

Примерно такая же картина имеет место при увеличении вдвое числа узлов (N = 32).

На фиг. 2 приведен для N = 32 пример восстановления решения модельной задачи, в которой ошибка экспериментальных данных равна ~10% (Д = 0,1).

\ A=J2;A^%f

•л

•V X

X > X X — точное решение х метод р" • регуляризация п-2

X >

\ • \ • \ “ X > X* X

х X ч • 1 • X X X

. . А X

4 < X а V X '

Фиг. 2

Видно, что при таких ошибках решение, полученное методом квадратур без предварительного сглаживания экспериментальных данных, является неудовлетворительным (ошибка восстановления ~32%). Применение регуляризации второго порядка позволяет восстановить искомое гладкое решение вдвое точнее (ошибка ■—15—16%).

В заключение отметим, что расчет по предлагаемому алгоритму для одного значения а требует для N=16 менее 1 с машинного времени.

ЛИТЕРАТУРА

1. Тихонов А. Н. Решение некорректно поставленных задач и метод регуляризации. Докл. АН СССР, т. 151, № 3, 1963.

2. Тихонов А. Н. О регуляризации некорректно поставленных задач. Докл. АН СССР, т. 153, № 1, 1963.

3. Тихонов А. Н. О некорректно поставленных задачах. В сб. „Вычислительные методы и программирование*, вып. 8. М., Изд-во МГУ, 1967.

4. Альберг Дж., Нильсон Э., УолшДж. Теория сплайнов и ее приложения. М., „Мир", 1972.

5. Физические измерения в газовой динамике и при горении. М., Изд. иностр. лит., 1957.

6. Возли некий М. И., Казанджан Э. П. Результаты применения интерференционного метода для расчета поля плотностей около некоторых тел вращения при осесимметричном сверхзвуковом обтекнии. Труды ВВИА, № 1137, 1967.

7. Морозов В. А. О регуляризующих семействах операторов. В сб. „Вычислительные методы и программирование”, вып. 8. М., Изд-во МГУ, 1967.

Рукопись поступила 241Ш 1975 г.