О сходимости разностной схемы, аппроксимирующей одну краевую задачу гиперболического типа

Челябинский физико-математический журнал. 2019. Т. 4, вып. 3. С. 333-344.

УДК 519.633 Б01: 10.24411/2500-0101-2019-14306

О СХОДИМОСТИ РАЗНОСТНОЙ СХЕМЫ, АППРОКСИМИРУЮЩЕЙ ОДНУ КРАЕВУЮ ЗАДАЧУ ГИПЕРБОЛИЧЕСКОГО ТИПА

А. С. Сушков

Челябинский государственный университет, Челябинск, Россия mesocyclon@yandex.ru

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

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

Введение

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

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

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

Заметим, что рассматриваемая в данной работе задача при е = 1 может описывать колебания круглой мембраны в среде, оказывающей сопротивление, пропорциональное скорости [2].

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

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

1. Постановка задачи и разностная схема

В области < = (0,1) х (0,Т) рассматривается следующая краевая задача с однородными граничными условиями для уравнения гиперболического типа:

д2м + дм д2м +1 дм о < ^ < т дг2 дг дх2 х дх' '

0 < х < 1, е > 0;

м(х, 0) = 0, м(х, 0) = /(х), 0 < х < 1;

м(0, г) = 0, м(м) = 0, 0 < г < т.

При заданной функции /(х) € С[0,1] ищется такое решение м(х, г) € С ([0,1] х [0, Т]) задачи, что %, §, 0, ^ € С ((0,1) х (0,Т)).

Перейдём к построению разностной схемы на равномерной сетке в прямоугольной области Q =

(xm, гп) € ^Кт.

Запишем задачу (1)-(3) в операторной форме Ьм = ^, где

м

Ьм

д2м дм д2м 1 дм

е--1-------,

дг2 дг дх2 х дх'

0 < г < Т, 0 < х < 1; м(х, 0), 0 < х < 1; м4(х, 0), 0 < х < 1; м(0,г), 0 < г < Т; м(м), 0 < г < Т;

0, 0 < г < Т, 0 < х < 1; 0, 0 < х < 1; /(х), 0 < х < 1;

0, 0 < г < Т; 0, 0 < г < Т.

В разностных производных задача (1)-(3) имеет вид уравнения

мп+1 _ 2мп + мп-1 мп+1 _ мп

е мт_2мт + "т + мт__т

м

п+1 _ 2мп+1 + мп+1

т+1 2мт + мт-1

т

2

т

к2

+

1 мп+1 _ мп+1 1 мт+1 мт-1

хт

2к

где т =1, 2,... , М — 1, п =1, 2,... , N — 1, с граничными условиями

мп = мМ = 0, п = 0,1,..., N,

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

м0 = 0,

мт мт

т

/т

т

0, 1, . . . , М.

2. Исследование аппроксимации разностной схемы

Невязкой метода (4)-(6) назовём сеточную функцию f(h) = Lh[u]h — f(h), где [u]h — точное решение задачи (1)-(3) в узлах сетки, f(h) — приближённое решение разностной схемы (4), Lh — разностная аппроксимация оператора L. В нашем случае невязка примет вид

u(xm,tra+l) 2u(xm, ^га) + u(xm,tra-1) . u(xm,tra+l) u(xmo tn) £-ö--1--

f

u(xm+1,tn+1) — 2u(xm,tn+1) + u(xm-1,tn+1) h2

1 u

T

x m 2 h

u

, m = 1, 2,..., M — 1, n = 1,...,N — 1,

— f (xm),m = 0,1,...,M,

и(х0, £п), п = 0,1,..., Ж, и(жм, £п), п = 0,1,..., N.

Будем говорить, что разностная схема (4) аппроксимирует задачу (1)-(3) с порядком /, если существует такая константа С, не зависящая от Л, что ||/< СЛ/, и

при этом ||//г~8'0> 0, где ^ — линейное нормированное пространство функций /с нормой ||  .

Лемма 1. Пусть точное решение и(х,£) задачи (1)-(3) таково, что функция и(х,£) трижды непрерывно дифференцируема по £ и четырежды непрерывно дифференцируема по х. Тогда разностная схема (4) аппроксимирует задачу (1)-(3) с первым порядком по времени и со вторым по координате.

Доказательство. Пользуясь разложением функции и(х,£) по формуле Тейлора в окрестности точки (хт,£га+1), получаем

д3и(хт,£) т2 д3и(хт,£) т2 д2и(хтД) —ет-—--+

f

dt3

6

dt3

2

dt2

h2 d4u(x,tn+1) h2 d3u(x,tn+1) ~ „

— M---—-, i G [tn-1,tn+1] ,x G [Xn-1, Xn+1]

12

dx4

6

dx3

< 0;

T д2u(xm,¿3) 2 dt2 : 0; 0.

Для

f

' f1(xm,tn),m = 1, 2,..., M — 1,n = 1, 2,...,N — 1; f2(xm),m = 0,1,...,M; f3(xm),m = 0,1,...,M; f4(t„),n = 0,1,...,N; L f5(tn),n = 0, 1,...,N,

введём норму

lf (h)l

max | f 1 (xm, tn) | + max | f2 (xm) | + max | f3 (xm) | + max | f4 (t„) | + max | f5 (tn

Если

то

д3м

дг3

< М

д 4м

дх4

< м14),

д3м

дх3

к2

< м!3),

д2м

дг2

к2

< М4(2) у(х,г) € <3,

II/(К)II* < етМ,(3) + — М4(3) + — МХ4) + М-М<3) + -М2) + —М4( . \\ J \\рк < * 1 6 4 12 х 6 х 2 4 2 4

Отсюда следует, что при т, к ^ 0 имеем /(К) ^ 0. В итоге получаем, что порядок аппроксимации задачи (4)-(6) равен 0(т + к2).

3. Исследование разностной схемы на устойчивость

Существование и единственность решения разностной задачи (4)-(6) следуют из того, что она представляет собой систему линейных алгебраических уравнений с невырожденной трёхдиагональной квадратной матрицей, имеющей диагональное преобладание.

Теорема 1. Пусть для параметров разностной задачи (4)-(6) выполнены неравенства к2 < 16е, к > 0, е > 0. Тогда решение задачи (4)-(6) м(хт,гп) € шКт при т, к ^ 0 сходится к точному решению м(х,г) задачи (1)-(3).

Доказательство. Обозначим через иК линейное нормированное пространство функций м(К) с нормой || • . Для устойчивости (4) необходимо доказать неравенство ||м(К)||иь < С1|/(К)||^ при некотором выборе норм. Перепишем начальные условия (6) в виде мт = 0, мт = т/т или же м0(х) = 0, м1(х) = т/(х), где х € ¿>К.

Воспользуемся принципом замороженных коэффициентов [5], выбирая на сетке произвольную точку х € Сс>к : 0 < х < 1. Решение (4)-(6) будем искать в виде

м

ХтТп ф 0, подставляя его в (4), имеем

ХтТп+1 2ХтТп + ХтТп— 1 , ХтТп+1 ХтТп.

е-о--1--

т2

т

= Хт+1Тп+1 — 2ХтТп + Хт-1Тп+1 + 1 Хт+1Тп+1 — Хт-1Тп+1

= к2 + х 2к '

Разделяя обе части равенства на ХтТп+1, получим

Тп.+1 2Тп + Тп— 1 . Тп+1 Тп

е-—--+ —

т 2Т

п+1

тТ

п+1

Хт+1 2Хт + Хт- 1 . 1 Хт+1 — Хт-1

к2 Хт х 2кХт

-Л.

Принимая во внимание условия (5), для нахождения Хт воспользуемся вспомогательной краевой задачей

(2х + к) Хт+1 + (2Лк2х — 4х) Хт + (2х — к) Хт-1 = 0;

Хо = X

М

0.

Частные решения ищем в виде Хт = дт. После подстановки в уравнение (8) получим (2х + к) д2 + (2Лк2х — 4х) д + (2х — к) = 0. Дискриминант этого квадратного уравнения есть Д = 4Л2к4х2 + 4к2 — 16Лк2х2. Пусть Д < 0, тогда

91,2 =

— (2Лк2х — 4х) ± 2гУ—Л2к4х2 — к2 + 4Лк2х2 2 (2х + к) '

2

Общее решение (8) имеет вид Xm = ciq^ + c2q™. Имеем ф = arg qi,

. , V-A2h4X2 - h2 + 4Ah2X2 , 2X - Ah2X . . /2Х - h sin ф =---, cos ф = , |q1,2| =

V2X - hV2Х + h ' V2X - hV2X + h' V 2x + h'

q™ = |qi|m (cos тф + i sinтф), qm = |q2|™ (cos тф — i sinтф).

Таким образом, Xm = c1 cos тф + c2 sin тф. В силу граничных условий в (8) получаем с1 = 0, с2 sin Мф = 0. Поскольку Xm ^ 0, то с2 = 0, sin Мф = 0, откуда Мф = nk, k = 1,... , М — 1. В итоге фк = nk/M = nkh. Поскольку с2 = 0 произвольно, то для простоты запишем X^ = sin rnnkh.

Обозначая x = rnh, можно принять X(k) (x) = sin nkx G H, где H — гильбертово пространство размерности M — 1, и при всех u, v G H (u, v) = ^f-1 u^h. Поэтому система функций X(k) (x) является в пространстве H ортогональной (см. [6]), причём ||X(k)||H = 1/л/2, поэтому система функций X(k) (x) = v^sinnkx в H является ортонормированной и образует базис.

Заметим, что при D > 0 задача (8) имеет лишь тривиальное решение. Найдём собственные значения A&, используя равенство

2 ,, —A2h4x2 — h2 + 4Ah2x2

sin nkh =---—-—-,

(2x — h) (2x + h) '

которое можно переписать в виде квадратного уравнения

A2h4x2 — 4Akh2x2 + h2 + (2x + h) (2x — h) sin2 nkh = 0.

Принимая во внимание, что 2x > h, получим D = 4x2h4 (4x2 — h2) cos2 nkh > 0 и, кроме того, по теореме Виета

2x — V4x2 — h2 cos nkh

Afc = -ГТ2- > 0.

xh2

Решение разностной схемы (4) представимо в виде ряда Фурье по системе X(k) (x) = л/^sin nkx в форме

м-1

un(x)= J] Xlfc)X(*>(*).

й=1

Для нахождения , используя (7), запишем разностное уравнение второго порядка (е + т + т2А&) — (2е + т) + еТ^— = 0. Его общее решение будет иметь вид т^ = + .

Частное решение этого уравнения будем искать с помощью подстановки Т,!^ = которая приводит к квадратному уравнению (е + т + т2А&) д2 — (2е + т) д+е = 0, откуда получаем

_ 2е + т ± тVI — 4А*е _ 2е + т ± гте — 1 91,2 = 2 (е + т + т2А^) = 2 (е + т + т^) .

Как и раньше, аналогичным образом находим, что

, , I е . те — 1 (2е + т) V

|91,21 = \ . . 2, , ф = а^д1, втф = =, сое ф = =,

V е + т + т2АЛ л/е + т + т2Ай ^е + т + т2Ай

q™ = |qi|n (cos пф + i sin пф) , qn = |q2|n (cos пф — i sinпф). При этом тПй) = cTi(fc) cos пф& + íí2(fc) sin пф&.

Поскольку u0(x) = ^l—1 T0(k)X(k)(x), то T0(k) — коэффициенты Фурье разложения u0(x) = 0 по системе собственных функций X(k)(x), откуда T0(k) = 0, поэтому di(k) = 0.

Аналогично, для u!(x) = M=11 Tik)X(k)(x) = г/(x) коэффициенты Фурье

^fc

1 ^(x)) = т (/(x),X (k)(x

Поэтому

Ti(k) = (т/(x), X (k)(x)) = т (/(x), X (k)(x)) = T/(k).

/ /(k)Ve + т + т2Afc _(fc) /(fc)Ve + т + т2Afc .

u2 =-— =-^ , ==—, ll ) =- , ==— sinmpk.

sin V^V4Afce - 1 V^V 4Afce - 1

В итоге решение разностной схемы имеет вид

wra

Поскольку ura(x) представлена рядом Фурье по ортонормированной системе собственных функций, то в ||-||h = ||-|| имеет место равенство Парсеваля

М-1

£ С (4Ak в — 1) sin2 < g е (4Akе — 1) '

Верно, что 2x > \J(2x + h) (2x — h) , поэтому

:2 — h2 cos nkh 2\

x h2 xh2

2x — V4x2 — h2 cos nkh 2^4x2 — h2 2 nkh Ak =-- > -ггт-sin

Обозначая

2 . o nkh

Ak = — sin2

h2 2

имеем, что

M—i /А

KÍ < --г-= £V 7

eAfc (4e - ¿) t! eW4X2 - h2 (4eAfc - 1) '

Поскольку Ak > 0, то

, . 2X — V4X2 — h2 cos nkh 4

Ak = |Ak|S -X-- 5 fti •

откуда следует, что

,,„„^ v1 if!!! x+т2) " " - Ak i^P^ft2

После преобразований получаем

2< í 4x (gh2 + тh2 + 4т2) \ ^1 (/(k)) < V h2£V4x2 — h2 (16е — h2V kL Ak

Здесь сумма

м-i f(fc)y

Е

At

к=1

задаёт норму, которая определяется следующим образом (см. [7]). Пусть и(х) =

: г>0 = ^м = 0 — сеточная функция на (А^, и

=-^-•

Тогда А-у = — Л^ — линейный ограниченный положительно определённый самосопряжённый сеточный оператор [8], который имеет обратный. Определим норму

v

A

-i

\J(A 1v(x),v(x)). Краевая задача Av = —Лv, v0 = vM = 0 имеет решения X(k) (x) = sinnkx при A = Ak, k = 1, 2,..., M — 1 [8].

Представим v

Поскольку АХ(к)(х) = АкX(к)(х), то X(к)(х) = АкА 1Х(к)(х), следовательно, А-1Х(к)(х) = X(к)(х)А-1. Тогда получим, что А-1^(х) = ^М- ^А-1Х(к)(х), поэтому

м 1 м 1

MlA-1 = (A-1v, v) = ( £ YkX(fc)(x)^ VsX(s)(x) ) =

M-1 M-1 „ M-1 s 2

Akvs (X(k)(x),X(s)(x)) = .

,Л .Л Ak ,л Ak

k=1 s=1 k=1

В итоге получаем

llu" II2 <

nn^

A-

-i

Далее, поскольку ||Л|| = AM-1 [9], то

A1 2 sin2 Пт '

В силу монотонного возрастания функции

h2

2 sin2 ^

2 nh

-1

как функции переменной 0 < h < 1 имеем, что ||A 1| < 1/2. Это позволяет согласовать н°рмы ц/ ¡A-i и ц/ Цс = max |/m| .

m=1,2,...,M

Используя неравенство Коши — Буняковского и тот факт, что /m = 0, получаем

M-1 M

m=1 m=1

M / \2

maxM /2h = ЦАЦ maxM /m hM = ПАп ma* M i/mi = ПАИ И/Mc .

m=1 2 ... M m=1 2 ... M m=1 2 ... M

m=1

Разморозим x G £Ah, заменяя его на максимальное значение hM =1. В результате получаем, что

П 4 (ftf + Wi2 + 4т2) ^

1 П < h2g^4—h2 (16g — h2r П

с •

Принимая во внимание найденную ранее оценку ||А 1| < 1/2, окончательно получим, что

2 (ек2 + тк2 + 4т2) м " к2е^4—к (16е — к2) Н/|с .

Из последнего неравенства и теоремы Лакса [6] получаем утверждение теоремы.

4. Численный эксперимент

Выполнена программная реализация метода (4)-(6) для е = 0.1, /(х) = х+1 при 0 < х < 1, 0 < г < 0.5. Решение будем находить на равномерной сетке. На рис. 1 представлена визуализация решения задачи при М = 2N = 2000. Визуализация производится при помощи программного компонента 1П;егае1луеВа1аВ18р1ау у.1.0.0 [10].

Значения функции и

Рис. 1. Визуализация приближённого решения при М = 2М = 2000

Заметим, что точное решение задачи (1)-(3) невозможно представить в виде совокупности элементарных функций [2], поэтому ограничимся рассмотрением того, как будут изменяться значения

^(М) = тах м(хт,гп), М = 100 + 10к, к = 0,1,..., 290,

М

функции, определяющей приближённое решение в узлах сетки при изменении числа узлов в указанных пределах (рис. 2). Из графика можно предположить, что при увеличении числа узлов метод (4)-(6) сходится к некому точному решению. Аналогичные результаты можно получить, зафиксировав некоторую точку на сетке

СКт.

Рис. 2. График зависимости ф(Ы), M = 100 + 10k, k = 0,1,..., 290

Список литературы

1. Кутузов, А. С. Оценка приближённого решения одной двумерной граничной обратной задачи тепловой диагностики методом квазиобращения / А.С.Кутузов // Вестн. Юж.-Урал. гос. ун-та. Сер.: Математика, физика, химия. — 2009. — № 10 (143). — С. 14-21.

2. Полянин, А. Д. Справочник по линейным уравнениям математической физики / А. Д. Полянин. — М. : Физматлит, 2001. — 576 с.

3. Кутузов, А. С. Моделирование процессов тепломассопереноса, описываемых дифференциальными уравнениями в частных производных на одно- и двумерных областях: дисс. ... канд. физ.-мат. наук / А. С. Кутузов. — Челябинск, 2013.

4. Годунов, С. К. Решение систем линейных уравнений / С.К.Годунов. — Новосибирск : Наука, 1980. — 177 с.

5. Годунов, С. К. Разностные схемы (введение в теорию) : учеб. пособие / С. К. Годунов, В. С. Рябенький. — М. : Наука, 1977. — 440 с.

6. Самарский, А. А. Теория разностных схем / А. А. Самарский. — М. : Наука, Гл. ред. физ.-мат. лит., 1989. — 616 с.

7. Крылов, В. И. Начала теории вычислительных методов. Уравнения в частных производных / В. И. Крылов, В. В. Бобков, П. И. Монастырный. — Минск : Наука и техника, 1986. — 311 с.

8. Самарский, А. А. Введение в теорию разностных схем / А. А. Самарский. — М. : Наука, 1971. — 553 с.

9. Колмогоров, А. Н. Элементы теории функций и функционального анализа / А. Н. Колмогоров, С. В. Фомин. — М. : Наука, 1972. — 496 с.

10. URL: https://github.com/Microsoft/InteractiveDataDisplay.WPF (дата обращения: 18.04.2019).

Поступила в редакцию 02.04.2018 После переработки 22.07.2019

Сведения об авторе

Сушков Андрей Сергеевич, аспирант математического факультета, Челябинский государственный университет, Челябинск, Россия; e-mail: mesocyclon@yandex.ru.

Chelyabinsk Physical and Mathematical Journal. 2019. Vol. 4, iss. 3. P. 333-344-

DOI: 10.24411/2500-0101-2019-14306

THE CONVERGENCE OF A DIFFERENCE SCHEME APPROXIMATING A BOUNDARY VALUE PROBLEM OF THE HYPERBOLIC TYPE

A.S. Sushkov

Chelyabinsk State University, Chelyabinsk, Russia mesocyclon@yandex.ru

The difficulty in applying differential problems in practice lies mainly in the impossibility of obtaining their solutions in an analytical form, which makes the development of numerical methods relevant. In this paper, we constructed one implicit difference scheme approximating a boundary value problem of hyperbolic type with homogeneous boundary conditions. The order of approximation of the difference scheme is found. Particular attention is paid to the proof of its stability and convergence. In the proof, an approach similar to the method of separation of variables in mathematical physics was used. The author found the convergence condition imposed on the parameters of the difference scheme. A numerical experiment is carried out. A program has been developed to find and visualize an approximate solution.

Keywords: hyperbolic type equation, boundary value problem, implicit difference scheme, approximation, difference scheme stability, difference scheme convergence.

References

1. KutuzovA.S. Otsenka priblizhyonnogo resheniya odnoy dvumernoy granichnoy obratnoy zadachi teplovoy diagnostiki metodom kvaziobrashcheniya [Estimation of the approximate solution of a two-dimensional boundary inverse problem of the thermal diagnostics by the quasi-reversibility method]. Vestnik Yuzhno-Uralskogo gosudarstvennogo universiteta. Ser.: Matematika. Fizika. Khimiya [Bulletin of South Ural State University. Ser.: Mathematics. Physics. Chemistry], 2009, no. 10 (143), pp. 4551. (In Russ.).

2. PolyaninA.D. Spravochnik po lineynym uravneniyam matematicheskoy fiziki [Handbook of mathematical physics linear equations]. Moscow, Fizmatlit Publ., 2001. 576 p. (In Russ.).

3. KutuzovA.S. Modelirovaniye protsessov teplomassoperenosa, opisyvayemykh differentsial'nymi uravneniyami v chastnykh proizvodnykh na odno- i dvumernykh oblastyakh [Modelling of heat and mass transfer processes, described by partial differential equations on one- and two-dimensional domains. PhD thesis]. Chelyabinsk, 2013. (In Russ.).

4. Godunov S.K. Resheniye sistem lineynykh uravneniy [Solving linear equation systems]. Novosibirsk, Nauka Publ., 1980. 177 p. (In Russ.).

5. Godunov S.K., Ryaben'kii V.S. Raznostnye skhemy (vvedeniye v teoriyu) [The difference schemes (introduction to the theory): a tutorial]. Moscow, Nauka Publ., 1977. 440 p. (In Russ.).

6. Samarskii A.A. Teoriya raznostnykh skhem [The theory of difference schemes]. Moscow, Nauka Publ., 1989. 616 p. (In Russ.).

7. KrylovV.I., BobkovV.V., Monastyrnyi P.I. Nachala teorii vychislitel'nykh metodov. Uravneniya v chastnykh proizvodnykh [The beginning of the theory of computational methods. Partial differential equations]. Minsk, Nauka i Tekhnika Publ., 1986. 311 p. (In Russ.).

8. Samarskii A.A. Vvedeniye v teoriyu raznostnykh skhem [Introduction to the theory of difference schemes]. Moscow, Nauka Publ., 1971. 553 p. (In Russ.).

9. Kolmogorov A.N., FominS.V. Elementy teorii funktsiy i funktsional'nogo analiza [Elements of the functions theory and functional analysis]. Moscow, Nauka Publ., 1972. 496 p. (In Russ.).

10. Available at https://github.com/Microsoft/InteractiveDataDisplay.WPF, accessed 18.04.2019.

Accepted article received 02.04-2018 Corrections received 22.07.2019