Решение задачи Дирихле для уравнения колебаний струны методом сопряженных градиентов Текст научной статьи по специальности «Математика»

УДК 519.63

В. И. Васильев, А. М. Кардашевский, В. В. Попов

РЕШЕНИЕ ЗАДАЧИ ДИРИХЛЕ ДЛЯ УРАВНЕНИЯ КОЛЕБАНИЙ СТРУНЫ МЕТОДОМ СОПРЯЖЕННЫХ ГРАДИЕНТОВ

Рассматривается неклассическая задача Дирихле для уравнения колебания струны, описываемое начально-краевой задачей для гиперболического типа второго порядка. Для ее численного решения используется итерационное уточнение начального условия для производной по времени, предложенное в работе [1] для решения ретроспективной обратной задачи теплопроводности. В данной работе предлагается наиболее быстро сходящийся итерационный метод сопряженных градиентов. Приведены примеры расчетов для модельных задач, в том числе и со случайными погрешностями во входных данных.

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

V. I. Vasilyev, A M. Kardashevsiy, V V Popov

The Dirichlet Problem for Oscillation Equation of String by Conjugate Gradient Method

The article considers a non-classical Dirichlet problem for the oscillation equation of string described by an initial-boundary value problem for the second-order hyperbolic equation. For numerical solution there is used the iterative improvement of the initial condition for the time derivative, which proposed in [6] for the solution of a retrospective inverse heat conduction problem. The author proposes the most rapidly converging iterative conjugate gradient method. Examples of calculations for model problems including one with random errors in input data are presented.

Keywords: inverse problems, second-order hyperbolic equation, Dirichlet problem, finite-difference method, explicit difference scheme, conjugate gradient method, random errors.

ВАСИЛЬЕВ Василий Иванович - д. ф.-м. н., профессор, зав. каф. вычислительных технологий ИМИ СВФУ им. М. К. Аммосова.

E-mail: [email protected]

VASILYEV Vasiliy Ivanovich - D. Sc., Professor, Head of the department «Computational Technologies» of Institute of Mathematics and Informatics, North-Eastern Federal University n. a. M. K. Ammosov.

E-mail: [email protected]

КАРДАШЕВСКИЙ Анатолий Михайлович - зав. кабинетом каф. «Вычислительные технологии» ИМИ СВФУ им. М. К. Аммосова.

E-mail: [email protected]

KARDASHEVSKIY Anatoliy Michaylovich - Head of Cabinet of the department «Computational Technologies» of Institute of Mathematics and Informatics, North-Eastern Federal University n. a. M. K. Ammosov.

E-mail: [email protected]

ПОПОВ Василий Васильевич - доцент-исследователь каф. вычислительных технологий ИМИ СВФУ им. М. К. Аммосова.

E-mail: [email protected]

POPOV Vasily Vasilievich - Associate professor and researcher of the department «Computational Technologies» of Institute of Mathematics and Informatics, North-Eastern Federal University n. a. M. K. Ammosov.

E-mail: [email protected]

Введение

Задача Дирихле для уравнения колебаний струны относится к классу условно корректных задач математической физики [2-3]. Доказательство корректности постановки задачи Дирихле для уравнений гиперболического типа второго порядка можно найти в книгах [24]. Общий подход к решению таких некорректных задач базируется на использовании градиентных итерационных методов. Основу этих методов составляет нахождение приближенного численного решения некорректных задач из итерационной минимизации соответствующего функционала. В цикле работ С. И. Кабанихина и его учеников [5-6] численно исследовано решение задачи Дирихле для волнового уравнения методами итераций Ландвебера и наискорейшего спуска, проведено теоретическое исследование их устойчивости и сходимости приближенного решения к решению исходной дифференциальной задачи для одномерного и двумерного уравнений гиперболического типа второго порядка.

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

ционныи метод решения неклассическои задачи для уравнения колебании струны, в котороИ по заданному распределению смещения струны в конечный момент времени определяется скорость струны в начальный момент времени.

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

Решение прямой задачи для уравнения колебаний струны

Колебание тонкой струны описывается дифференциальным уравнением с частными производными гиперболического типа второго порядка

yí

J+1

2 yi

-yi 1 _ y,+1 - 2 y,

f -f? = »• x ^ 0 < < *

где

^ {xi = ih,/=0,1,...,n;h = l/n},

i = 1, n -1,j = 1, m — 1,

(4)

(5)

(6)

(1)

у0+1 = 0, уП+1 = 0, ] = 1,т -1, у0 = р(х,) у) = у0 +™(х,), I = 0,n,

здесь использовано обозначение

У, = °уГ +(1- 2Я)У1

На множестве сеточных функций у е Н таких, что у(х)=0, х ^ тк, определим сеточный оператор А соотношением

АУ = -Ухх, х ^ . Хорошо известно [9], что для решения разностной схемы (4)-(6) при выполнении условия устойчивости

^ 1 1

" й 4 "ТЛА1 (7)

справедлива следующая априорная оценка:

11?Ц<1т1<...<11^(г)1, г е ч, (8)

норма || У ||* определяется следующим образом:

В предположении, что концы колеблющейся струны х=0, х=1 закреплены, выполняются однородные граничные условия первого рода

и (0, г) = 0, и (I, t) = 0, 0 < t < Т. (2)

В начальный момент времени заданы смещение и скорость точек струны

и(х,0) = ф(х), (х,0) = у(х), 0 < х < I. (3)

При достаточной гладкости входных данных ф(х), v(x) и выполнении условий согласования начальных и граничных условий задача (1)-(3) поставлена корректно [8].

Для ее численного решения используем метод конечных разностей. В области определения задачи (1)-(3), представляющей собой прямоугольник [0,1]х[0,Г], введем равномерную прямоугольную сетку

^2 1 - 2 2 П = Т11У + У |

+Т

yt II2 1 , R =1E + aA. ' r—A 4

4

Тем самым норма решения задачи Коши со временем не возрастает. Отметим, что для явной разностной схемы (а=0) условие устойчивости (7) принимает простейший вид т<к. Таким образом, на решении исходной дифференциальной задачи (1)-(3) разностная схема (4)-(6) обладает вторым порядком аппроксимации, и при выполнении условия

а <--" 4

i h

4т

юг = { = Д, 7=0,1,..., т;т = Т / т}.

Задаче (1)-(3) на пространственно-временной сетке юы поставим в соответствие симметричную разностную схему с весовым множителем а:

ее решение при достаточной гладкости входных данных сходится к решению дифференциальной задачи (1)-(3) со скоростью порядка 0(r2+h2). Приведем результаты вычислительного эксперимента. Сначала проведем расчеты, показывающие влияние шагов сетки на точность определения приближенного решения прямой задачи для уравнения колебаний струны. Для этого рассмотрим задачу с точным решением

u(x,t) = sin(x)cos(t), x e[0,l], t е[0,Г],

где l=n, T=n/2.

На рис. 1 приведен график зависимости от шагов пространственно-временной сетки точности решения явной разностной схемы на решении исходной задачи в финальный момент времени

2

2

h

т

Рис. 1. Решение прямой задачи для уравнения колебаний струны

||2™ ||=^В ут

о)2 а.

Расчеты проводились на последовательности сеток, на графике количество узлов пространственной сетки через к вычисляется по формуле я=20-2к, т. е. при п=20, 40, 80, 160, 320, 640. Сплошная линия (1) соответствует расчету с максимально допустимым временным шагом, обеспечивающим устойчивость разностной схемы т=h, следовательно, т принимает те же значения, что и п. Прерывистой линией (2) представлена та же зависимость при т=Ы2, а штрих-пунктирной линией (3) - при т=Ы4. То, что графики представляют собой прямые линии, подтверждает теоретическое утверждение о втором порядке сходимости по h и т решения разностной схемы к решению исходной дифференциальной задачи. Представленные графики также показывают, что при решении задачи Дирихле критерий выхода из итерационного цикла на данной конкретной пространственно-временной сетке следует согласовывать с точностью решения прямой задачи.

Численную реализацию разностной схемы (4)-(6) проведем для задачи (1)-(3) с однородными граничными условиями, начальным условием для скорости смещения струны 1(0=0, ^2(/)=0, 0</<Г, ф(х)=0, 0<х<1 и при неоднородном начальном условии для смещения

ф) = 12е-50( /2)2

пространственно-временная сетка построена таким образом, что условие устойчивости Л выполняется с минимальным запасом. Справа на рис. 2 представлено решение у/, I = 0,п, j = 0,т разностной схемы (п=200, т=200).

Результаты численного эксперимента показывают, что достаточно хорошую точность получаем на достаточно грубых сетках с п=50, т=50. В дальнейшем для задания значения решения в конечный момент времени для неклассической задачи будем использовать решение прямой задачи в финальный момент времени, полученное на достаточно подробной сетке.

Постановка задачи Дирихле для уравнения колебаний струны

В прямоугольнике ищем функцию и(х,0 - решение гиперболического уравнения второго порядка

д2ы

д

дГ дх

k (х)ди

У дх

= о, х е(0, I), о < t < т.

(9)

0 < X < 1.

На рис. 2 приведены результаты решения прямой задачи по явной разностной схеме (о=0) при /=1, 7=0,99 на последовательности сгущающихся сеток с (п, т)=(25, 25); (50, 50); (100, 100); (200, 200). Здесь слева даны графики функции ф(х), а в середине - график решения в финальный момент времени у™, г = 0,п . Отметим, что

Предположим, что концы колеблющейся струны х=0, х=1 закреплены, т. е. выполняются граничные условия первого рода

и (0, t)= 0, и (I, t) = 0, 0 < t < Т. (10)

В начальный и финальный моменты времени задаются смещения струны

и (х,0) = ф(х), и (х, Т) = ф(х), 0 < х < I. (11)

Пусть коэффициент к(х) дифференциального уравнения (7) удовлетворяет условиям 0<к1<к(х)<к2. Задача Дирихле относится к классу неклассических задач математической физики, является условно корректной [3]. Она может иметь не единственное решение. Как показывает нижеследующий простейший пример. Пусть

Рис. 2. Решение прямой задачи для уравнения колебаний струны

фх)=0, ^(x)=0, тогда неклассическая задача (9)-(11) при коэффициенте k(x)=0, помимо тривиального решения u(x,t)=0, имеет семейство решений, задаваемых формулой

u (x, t) = c • sin (kpx)- sin (npt), c € R, k, l, T € N. Разностная задача Дирихле

Обозначим через ты множество внутренних узлов пространственно-временной сетки

®Ы = ®h х®

где

ah = {x¡ = ih,/=1,2,..., n-1},

юг ={tj = jz, j=1,2,..., m-1}.

На множестве сеточных функций y е H таких, что _y(x)=0, х ^ ah, определим сеточный оператор A соотношением

Ay = -(ci (x) yx )x, x € mh, положив, например, a(x)=k(x-0,5h).

В сеточном гильбертовом пространстве Н скалярное произведение и норму введем соотношениями

(у>= |у ||= (у,у)1/2.

В Н оператор А является самосопряженным, положительно определенным А=А*>0 при к(х)>кр к>0 и ограниченным ||А||<4к2^2.

Сначала от задачи Дирихле перейдем к задаче определения решения дифференциально-операторного уравнения

d 2

—y + Ay = 0, X е а, 0 < t < T dt

(12)

при заданных дополнительных условиях, обеспечивающих единственность ее решения

У(х,0) = ф(х), у(х,Т) = ф(х), х е юь. (13)

При использовании симметричной трехслойной разностной схемы с весовым множителем а дискретный

аналог задачи (10)-(11) имеет вид:

Л V

У-2у + У + А\ау + (1-2а)у + ^1 = 0, (х,()е (14)

* (0)

У (п

Ь(х), х £ тк, (15)

где 0<ст<1 и использована безиндексная система обозначений, введенная А. А. Самарским [9]:

^+1

V

У:

У = У % У = У > у = У % * е ah.

Для численной реализации системы линейных алгебраических уравнений (12)-(13) воспользуемся итерационным методом сопряженных градиентов, связанным с уточнением производной по времени в начальный момент времени. На каждой итерации решается прямая задача для уравнения колебаний струны с использованием стандартной симметричной трехслойной разностной схемы, обладающей вторым порядком аппроксимации по h и т [9].

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

д(x,0) = v(x), х € ю.

Обозначим через у разностное решение на момент времени Р=/х, где т>0 - шаг по времени, причем тт=Т.

Дискретный аналог прямой задачи имеет вид:

У - 2 У + У

(x, I )е юк

А

а

уч V

у + (1-2а) у + а у

0,

(16)

В силу самосопряженности оператора А самосопряженным является и оператор Л,. Однако мы ничего не можем сказать об однозначной разрешимости данного сеточного уравнения.

Для решения этого уравнения используем трехслойный итерационный метод сопряженных градиентов, записанный в каноническом виде [9]

vk+1 = ак+Ук + (1-ак+1И-1 -ак+Л+Л, к = 0,1,...,

где гк = Лvk —ф- невязка.

Для вычисления итерационных параметров ^ используются следующие рекуррентные формулы

(г,, г,/

(18)

вк+1 =

а

к +1 _ к = 1,2

(Агк , г

вк+1

к = 0,1,...

(гк, гк) 1

(19)

(20)

Рк (Г'к-и гк-1 )ак , а1 = 1.

В этом итерационном методе вычисление уП осуществляется численной реализацией прямой задачи для очередного приближения искомой сеточной функции vk:

У к - 2 У к + У к

(х г )е юк

°Ук +(1-2°)Ук Ук

0,

с начальными условиями

Ук = 9 Ук = 9 + ™к, * € .

(21)

(22)

У (0) = ф(х,), У (т) = у (0) + ^(х), X е юк. (17) Итерационный метод

Для нахождения решения дискретной задачи Дирихле (12)-(13) используем наиболее быстро сходящийся итерационный метод вариационного типа - метод сопряженных градиентов [10], основанный на последовательном уточнении искомого начального условия v(x) с дальнейшим решением на каждой итерации дискретного аналога прямой задачи (16)-(17). Придадим этой задаче соответствующую операторную формулировку. Последовательно исключая промежуточные значения сеточной функции у (), t € при заданных ф и V на конечный момент времени, получим

ум = Av + Бф,

где Л, Б - соответственно операторные полиномы от положительно определенного и самосопряженного оператора А. С учетом этого факта приближенному решению обратной задачи естественно сопоставить решение следующего операторного уравнения

Av = ф(х) — Бф = ф.

Следовательно, невязка вычисляется по формуле гк = ук — ф, х € юк. Аналогично организуется вычисление вектора ъ=Атк:

V ^к

Zк -

-2 z■

_2

т

V

аzk + (1 — 2а) zk + а zk

(X, ? ) € Щ

кт

с начальными условиями

4 = °> 4 = тгк, х е Щ.

(23)

(24)

Таким образом, итерационный процесс реализуется в следующем порядке:

1. Полагаем к=0 и задаем начальное приближение искомой функции vk, х е юк.

2. Запускаем счетчик итераций к=к+1, последовательно решая прямые задачи (21)-(22) и (23)-(24), определяем, соответственно, невязку гк = ук — ф, х € юк и вектор zk = Лгк.

3. По рекуррентным формулам (19)-(20) определяем значения итерационных параметров вк+1, ак+1.

4. По формуле (18) находим очередное приближение искомого начального условия vk+1(х), х € юк.

2

Т

2

т

5. Этот процесс продолжается до тех пор, пока не выполнится критерий остановки итераций г<г, иначе возвращаемся к пункту 2.

Примеры расчетов

Вычислительную эффективность предложенного метода удобно проиллюстрировать на примере численного решения простейшей задачи Дирихле для уравнения колебания струны. Будем искать приближенное решение задачи Дирихле для уравнения

д2ы д2ы

дt дх

= 0, 0 < х < 1, 0 < I < Т,

с однородными граничными условиями

и (0, t) = 0, и (1, t) = 0, 0 < t < Т,

и условиями в начальный и финальный моменты времени

и (х, 0) = 0, и (х, Т) = ф(х), 0 < х < 1.

Пусть h - шаг равномерной сетки по пространственной переменной

юк = {х | х = хг = гк, г = 0,1,...,п; пк = 1}.

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

, , ди (х,0) , .

и(х,0) = 0, —^—!- = v(x), 0 < х < 1.

Из решения прямой задачи определяется функция ф (х), которая и присутствует в постановке задачи Дирихле:

и(х,Т) = ф(х), 0 < х < 1.

В качестве начального приближения искомого начального условия возьмем функцию, отдаленно напоминающую искомую функцию

v0 = 5sin (лх), 0 < х < 1.

Приведем результаты по решению обратной задачи Дирихле в условиях, когда искомая функция ф(х) имеет разный вид, на разных сетках, при разных Т. В расчетах задавали е=0,0001.

На рис. 3 приведены результаты в разные моменты времени т=100, 1=1, 7=0,24, 0,48, 0,72, 0,96, 1,2; п=25, 50, 100, 200. Решение в разные моменты времени. Соответственно, итераций 9, 8, 9, 5, 13, а время счета 0,332, 0,563, 0,899, 0,445, 1,436 с.

На рис. 4 приведены аналогичные результаты, когда искомая функция задана в виде

у(Х) = 7(-е-10°(х-1 /4)2 + е-20(х-1 /2)2 _ е-100(х-31 /4)2 ,

0 < х < 1.

В этом случае итерации сходятся значительно медленнее, соответственно, количество итераций 55, 48, 43, 39, 63, а время счета составляет, соответственно, 1,736, 2,012, 2,692, 3,261, 6,592 с.

Наконец, в численных экспериментах сеточная функция ф(х), х € соь возмущалась, как в работе [1], следующим образом:

Фз (х) = ф + 3а(х), X £ юь,

где ст(х) - равномерно распределенные на отрезке [-1, 1] случайные величины. Итерационный процесс сопряженных градиентов обрывался по достижению невязки значения д, т. е. при

45)1

< 5.

Проведены расчеты для уровня погрешности во входных данных, определяемых величиной д=0,001 (25). Здесь приведены данные расчетов при использовании предлагаемого итерационного метода без сглаживающего оператора. Для выделения более гладкого решения, принадлежащего вместо Ь2(ю) сеточному пространству Ж (ю), зададим сглаживающий оператор в виде

^ = -5Ухх + У, Х € ю сеточных функций, обращающихся в нуль в граничных узлах.

Графики функций v(x) и ф(х) показаны слева на рис. 5. Здесь приведены результаты счета на пространственной сетке М=200, h=0,005 при помощи явной схемы для волнового уравнения а=0 с временным шагом, удовлетворяющим условию устойчивости: г=0,00495^ 7=0,99, Ж=200.

В данном случае итерационный процесс сопряженных градиентов сходится достаточно быстро, всего за 5-15 итераций, поскольку каждый раз генератор случайных чисел выдает новый набор случайных чисел. В расчетах брали е=0,007. Очевидно, рассчитывать на точное восстановление начального v(x) нельзя в силу пониженной гладкости входных данных. С хорошей точностью приближенное решение на конечный момент времени мы имеем от более гладкого начального условия, которое находится при итерационном решении обратной задачи при сглаживании входных данных. В двумерном случае количество итераций значительно растет. В этом случае для существенного ускорения скорости сходимости итерационного процесса возникает необходимость использования в методе сопряженных градиентов предобуславливателя.

Заключение

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

п-25

Рис. 3. Решение в разные моменты времени

Рис. 4. Решение в разные моменты времени

Рис. 5. Восстановленное начальное условие

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

Авторы благодарны профессору П. Н. Вабищевичу за конструктивные замечания и советы.

Работа выполнена при финансовой поддержке РФФИ (Грант № 13-01-00719).

Л и т е р а т у р а

1. Самарский А. А., Вабищевич П. Н., Васильев В. И. Итерационное решение ретроспективной обратной задачи теплопроводности // Математическое моделирование. - 1997. - Т. 9, № 5. - C. 119-127.

2. Samarskii A. A., Vabishchevich P. N. Numerical Methods for Solving Inverse Problems of Mathematical Physics. De Gruyter, 2007.

3. Kabanikhin S. I. Inverse and ill - posed problems: theory and applications. De Gruyter, 2012.

4. Денисов А. М. Введение в теорию обратных задач. - М.: Изд-во Московского ун-та, 1994.

5. Kabanikhin S. I., Bektemesov M. A., Nurseitov D. B., Krivorotko O. I., Alimova A. N. An optimization method in the Dirichlet problems for the wave equation // Journal of inverse and ill - posed problems. 2012. - № 2 (20). - C. 193-211.

6. Кабанихин С. И., Криворотько О. И. Численный метод решения задачи Дирихле для волнового уравнения // Сиб. журн. индустр. матем. - 2012, - Т. 15, № 4. - С. 90-101

7. Васильев В. И., Попов В. В., Еремеева М. С., Кардашев-ский А. М. Итерационное решение одной неклассичесой задачи для уравнения колебаний струны // Вестник МГТУ им. Н. Э. Баумана. Сер. Естественные науки. - 2015. - № 3. - С. 77-87.

8. Тихонов А. Н., Самарский А. А. Уравнения математичес-

кой физики. - М: Наука, 1972.

9. Самарский А. А. Теория разностных схем. - М.: Наука, 1989.

10. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. - М.: Наука, 1978.

R e f e r e n c e s

1. SamarskiiA.A., Vabishchevich P. N., Vasil'ev V. I. Iteratsionnoe reshenie retrospektivnoi obratnoi zadachi teploprovodnosti // Matematicheskoe modelirovanie. - 1997. - T. 9, № 5. - C. 119-127.

2. Samarskii A. A., Vabishchevich P. N. Numerical Methods for Solving Inverse Problems of Mathematical Physics. De Gruyter, 2007.

3. Kabanikhin S. I. Inverse and ill - posed problems: theory and applications. De Gruyter, 2012.

4. Denisov A. M. Vvedenie v teoriiu obratnykh zadach. - M.: Izd-vo Moskovskogo un-ta, 1994.

5. Kabanikhin S. I., Bektemesov M. A., Nurseitov D. B., Krivorotko O. I., Alimova A. N. An optimization method in the Dirichlet problems for the wave equation // Journal of inverse and ill - posed problems. 2012. - № 2 (20). - C. 193-211.

6. Kabanikhin S. I., Krivorot'ko O. I. Chislennyi metod resheniia zadachi Dirikhle dlia volnovogo uravneniia // Sib. zhurn. industr. matem. - 2012, - T. 15, № 4. - S. 90-101

7. Vasil'ev V. I., Popov V. V., Eremeeva M. S., Kardashevskii A. M. Iteratsionnoe reshenie odnoi neklassichesoi zadachi dlia uravneniia kolebanii struny // Vestnik MGTU im. N. E. Baumana. Ser. Estestvennye nauki. - 2015. - № 3. - S. 77-87.

8. Tikhonov A. N., Samarskii A. A. Uravneniia matematicheskoi fiziki. - M: Nauka, 1972.

9. Samarskii A. A. Teoriia raznostnykh skhem. - M.: Nauka, 1989.

10. Samarskii A. A., Nikolaev E. S. Metody resheniia setochnykh uravnenii. - M.: Nauka, 1978.

^■Hir^ir