Сеточные схемы для решения уравнения гиперболического типа с запаздыванием в производной Е.Е. Таширова linetisa@yandex.ru УрФУ (Екатеринбург) Аннотация Построен численный метод для решения уравнения гиперболиче- ского типа с запаздыванием во второй производной по простран- ственной переменной и функциональным запаздыванием в неод- нородности. Исследован порядок локальной погрешности метода. Получена теорема о порядке сходимости алгоритма с помощью вложения в общую схему численных методов для функционально- дифференциальных уравнений. Приводятся результаты расчетов для тестового примера. 1 Постановка задачи Наличие запаздывания в производной может привести к новым эффектам, например, появлению коле- баний [1]. Численные алгоритмы для уравнений параболического типа с запаздыванием во второй произ- водной изучались в ряде статей, например [2]. Численный метод для уравнения переноса с запаздыванием в производной и функциональным запаздыванием в неоднородности был построен и исследован в работе [3]. Численные алгоритмы для уравнений гиперболического типа с запаздыванием в производной ранее не исследовались. При этом, численные методы для уравнений гиперболического типа с функциональным запаздыванием в неоднородности были сконструированы и исследованы в работе [4]. Рассмотрим уравнение ∂2u 2 2∂ u 2 2 2∂ u (x, t) = a (x, t) + a c (x, t − τ ) + f (x, t, u(x, t), ut (x, ·)) : t ∈ [t0 , T ], x ∈ [0, X] (1) ∂t2 ∂x2 ∂x2 с граничными: u(0, t) = g1 (t), u(X, t) = g2 (t) : t ∈ [t0 , T ] и начальными условиями u(x, t) = ϕ(x, t) : x ∈ [0, X], t ∈ [t0 − τ, t0 ]. Здесь x, t — независимые переменные, u(x, t) — искомая функция, ut (x, ·) = {u(x, t + ξ), −τ 6 ξ < 0)} — функция-предыстория искомой функции к моменту t, τ — величина запаздывания; f (x, t, u(x, t), ut (x, ·)) — функционал, определенный на [0, X] × [t0 , T ] × R × Q, Q = Q[−τ, 0) — множе- ство функций u(ξ), кусочно-непрерывных на [−τ, 0) с конечным числом точек разрыва первого рода, в точках разрыва непрерывных справа, ku(·)kQ = supξ∈[−τ,0) |u(ξ)|. Будем предполагать, что функционал f и функции g1 , g2 , ϕ таковы, что задача имеет единственное решение u(x, t). Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org 335 2 Численный метод Разобьем отрезок [0, X] на части с шагом h = X/N , где N — некоторое целое число. Введем точки xi = ih, i = 0, 1 . . . , N . Разобьем отрезок [t0 , T ] на части с шагом ∆. Будем считать, что m = τ /∆ – целое число. Введем точки tj = t0 + j∆, j = −m, . . . , M . Будем обозначать приближение точного решения u(xi , tj ) через uij . Введем дискретную предысторию к моменту tj при каждом фиксированном i: {uil }j = {uil : j − m 6 l 6 j}. Оператором интерполяции-экстраполяции дискретной предыстории назовем отображение: I : {uil }j → vji (·) ∈ Q[−τ, ∆]. Будем говорить, что оператор интерполяции-экстраполяции имеет порядок погрешности p на точном ре- шении, если существуют константы C1 и C2 , такие, что для всех i = 0, . . . , N, j = 0, . . . , M и t ∈ [tj − τ, tj+1 ] выполняется неравенство |vji (t − tj ) − u(xi , t)| 6 C1 max |uil − u(xl , ti )| + C2 ∆p . j−m6l6j Например, кусочно-линейная интерполяция 1 vji (ξ) = ((tl − tj − ξ)uil−1 + (tj + ξ − tl−1 )uil ), tl−1 6 tj + ξ 6 tl , −τ 6 ξ 6 0 (2) ∆ с экстраполяцией продолжением 1 vji (ξ) = ((−ξ)uij−1 + (∆ + ξ)uik ), tj 6 tj + ξ 6 tj + 1, ξ>0 ∆ имеет второй порядок [5]. Будем предполагать, что оператор интерполяции-экстраполяции липшицев, то есть найдется такая кон- станта LI , что для любых двух предысторий {uil }j и {νli }j выполненяется sup |vji (t − tj ) − ηji (t)| 6 LI max |uil − νli |, tj −τ 6t6tj +∆ j−m6l6j где vji (·) = I({uil }j ), ηji (·) = I({νli }j ). Также будем предполагать, что оператор интерполяции-экстраполяции согласован: vji (tl − tj ) = uil , l = j − m, . . . , j. Рассмотрим метод uij+1 − 2uij + uij−1 2 (ui+1 j − 2uij + uji−1 ) 2 2 (ui+1 i i−1 j−m − 2uj−m + uj−m ) = a + a c + Fji (vji (·)), ∆2 h2 h2 i = 1, . . . N − 1, j = 0, . . . M − 1 (3) с граничными uj0 = g1 (tj ), ujN = g2 (tj ) и начальными условиями uji = ϕ(xi , tj ) : − m 6 j 6 0. Здесь Fji (v(·)) — некоторый функционал, определенный на функциях v(·) = vji (·) = I({uik }j ) ∈ Q[−τ, ∆], связанный с функционалом f (xi , tj , uij , vji (·)) и липшицевый по переменной v(·) с константой LF . Метод явный, при каждом фиксированном j можно выразить uij+1 . Невязкой (без интерполяции) метода назовем: u(xi , tj+1 ) − 2u(xi , tj ) + u(xi , tj−1 ) u(xi−1 , tj ) − 2u(xi , tj ) + u(xi+1 , tj ) ψji = 2 − a2 − ∆ h2 u(xi−1 , tj−m ) − 2u(xi , tj−m ) + u(xi+1 , tj−m ) − a2 c2 − Fji (utj (xi , ·)). (4) h2 336 Будем говорить, что невязка имеет порядок hp1 + ∆p2 , если существует такая константа C, что |ψji | 6 C(hp1 + ∆p2 ) для всех i = 1, . . . N − 1, j = 0, . . . M − 1. Рассмотрим метод (3), в которых функционал Fji определяется следующим образом: Fji (vji (·)) = f (tj , xi , uij , vbji (·)), (5) где vbji (·) ∈ Q[−τ, 0); vbji (ξ) = vji (ξ), −τ 6 ξ < 0. Теорема 1. Если для точного решения задачи (1) существуют и непрерывны все частные производные вплоть до 4-го порядка включительно, Fji определяется соотношением (5), то невязка имеет порядок h2 + ∆2 . Доказательство. Невязка определяется равенством u(xi , tj+1 ) − 2u(xi , tj ) + u(xi , tj−1 ) u(xi−1 , tj ) − 2u(xi , tj ) + u(xi+1 , tj ) ψji = 2 − a2 − ∆ h2 u(xi−1 , tj−m ) − 2u(xi , tj−m ) + u(xi+1 , tj−m ) − a2 c2 − f (tj , xi , u(xi , tj ), utj (xi , ·)). h2 Разложим функцию u(x, t) по формуле Тейлора в окрестности точек (xi , tj ), (xi , tj−m ). Так как существу- ют и непрерывны все частные производные вплоть до 4-го порядка, получаем следующие равенства для значений функции в точках: ∂u 1 ∂2u 2 1 ∂3u u(xi , tj±1 ) = u(xi , tj ) ± (xi , tj )∆ + (x i , tj )∆ ± (xi , tj )∆3 + O(∆4 ), ∂t 2 ∂t2 6 ∂t3 ∂u 1 ∂2u 2 1 ∂3u u(xi±1 , tj ) = u(xi , tj ) ± (xi , tj )h + (x i , tj )h ± (xi , tj )h3 + O(h4 ), ∂x 2 ∂x2 6 ∂x3 ∂u 1 ∂2u 1 ∂3u u(xi±1 , tj−m ) = u(xi , tj−m ) ± (xi , tj−m )h + 2 (xi , tj−m )h2 ± (xi , tj−m )h3 + O(h4 ). ∂x 2 ∂x 6 ∂x3 Подставив эти соотношения в формулу для ψji , получим: ∂2u ∂2u ∂2u ψji = 2 (xi , tj ) − a2 2 (xi , tj ) − a2 c2 2 (xi , tj−m ) − f (tj , xi , u(xi , tj ), utj (xi , ·)) + O(∆2 + h2 ). ∂t ∂x ∂x В силу уравнения (1) ψji = O(∆2 + h2 ).  3 Исследование сходимости Обозначим величину погрешности метода в узлах через εij = u(xi , tj ) − uij . Будем говорить, что метод сходится с порядком hp + ∆q , если существует такая константа C, не зависящая от ∆, h, что выполняется неравенство |εij | 6 C(hp + ∆q ) для всех i = 0, . . . , N, j = 0, . . . M. Будем исследовать сходимость полученной схемы (3), используя общую дискретную схему с эффектом наследственности [5, 6, 7, 8]. Будем рассматривать задачи с однородными краевыми условиями: u(0, t) = u(X, t) = 0, t0 6 t 6 T. При каждом tj определим значения дискретной модели вектором γej = (u0j , u1j , u2j , . . . ujN −1 , uN 0 j ) ∈ Γ, где e 0 — знак транспонирования, Γ e — векторное пространство размерности N + 1 со скалярным произведением: N X −1 (e γ, ω e) = ei ω γ e i h, γ γ0, γ e = (e eN )0 ∈ Γ, e1 , . . . γ e ω ω0 , ω e = (e e N )0 ∈ Γ. e1, . . . ω e (6) i=1 337 В пространстве Γ e введем операторы A и A: e ui−1 j − 2uij + ui+1 j Aγj = ωj , ωj0 , ω ωj = (e ejN )0 , ej1 , . . . ω ωji = −a2 : 1 6 i 6 N − 1, ωj0 = 0, ωjN = 0, (7) h2 e = ∆2 A. A (8) Оператор A самосопряженный и положительный [9] в смысле скалярного произведения (6). Перепишем систему (1) в виде: ej+1 − 2e γ γj + γ ej−1 γj + c2 Ae + Ae γj−m = Fj (v(·)), (9) ∆2 где Fj (v(·)) = (Fj0 (vj0 (·)), Fj1 (vj1 (·)), . . . , FjN −1 (vjN −1 (·)), FjN (vjN (·)))0 , γk }j ) ∈ QN +1 [−τ, ∆], v(·) = I({e QN +1 [−τ, ∆] — пространство вектор-функций, каждая компонента которых принадлежит Q[−τ, ∆]. Приведем уравнение (9) к явной форме ej+1 = (2E − A)e γ ej−1 − c2 Ae γj − γ eγj−m + ∆2 Fj (v(·)). (10) В пространстве Γ e введем норму: p ke γn kΓe = (Ae γn , γ en ). (11) Заметим, что норма оператора A e (8), в отличие от нормы оператора A (7), не зависит от h и ∆, а лишь от σ = a2 ∆2 /h2 . Введем вектор γj = (γj0 , γj1 , . . . , γjm+1 )0 = (e γj−m , γ ej )0 ∈ Γ, где Γ — векторное пространство ej−m+1 , . . . , γ размерности q = (m + 1)(N + 1). Введем норму в пространстве Γ: e m X kγk2Γ = kγ k k2Γe . (12) k=0 В результате получаем разностную схему: γj+1 = Sγj + ∆Φ(tj , I({γk }j ), ∆), (13)   0 1 0 ... 0 0   0  0 0 1 ... 0 0  .. .. .. .. .. ..     где S =   , Φ(tj , I({γk }j ), ∆) =    . .  .. . . .   0   0  0 0 ... 0 1  γk }j ))/∆ Fj (I({e −c2 A e 0 0 . . . −1 2E − A e Определим функцию точных значений для схемы (13) соотношениями zj = (zj0 , zj1 , . . . zjm+1 )0 = (e zj−m , zej−m+1 , . . . , zej )0 , zej = (u(x0 , tj ), u(x1 , tj ), . . . , u(xN −1 , tj ), u(xN , tj ))0 . Стартовые значения модели (13) можно определить следующим образом: zj−m , zej−m+1 , . . . , zej )0 , γj = zj = (e zej = (ϕ(x0 , tj ), ϕ(x1 , tj ), . . . , ϕ(xN , tj ))0 , j = −m, . . . , 0. Будем говорить, что стартовые значения модели имеют порядок ∆p1 + hp2 , если найдется константа C, не зависящая от zj , yj , ∆, h, такая, что kzj − yj kY 6 C(∆p1 + hp2 ), j = −m, . . . , 0. Погрешность аппроксимация (невязка) с интерполяцией в общей разностной схеме [5, 6, 7, 8] определя- ется следующим образом: dn = (zn+1 − Szn )/∆ − Φ(tn , I({zi }n−m ), ∆), n = 0, . . . , M − 1. (14) 338 Будем говорить, что метод имеет порядок погрешности аппроксимации с интерполяцией ∆p1 + hp2 , если существует константа C, не зависящая от dn , ∆, h, такая, что kdn k 6 C(∆p1 + hp2 ), n = 1, . . . , M. Данное определение невязки отличается от введенного ранее определения невязки без интерполяции (4). Однако справедливо следующее утверждение. Теорема 2. Пусть невязка в смысле (4) имеет порядок ∆p1 + hp2 , функции Fji липшицевы, оператор интерполяции-экстраполяции I имеет порядок погрешности p0 на точном решении, σ зафиксировано, тогда невязка с интерполяцией в смысле (14) имеет порядок, равный ∆min{p0 ,p1 } + hp2 . Доказательство. По определению нормы в Γ (12): m X kdn k2Γ = kdkn k2Γe . k=0 Для k = 0, . . . , m − 1 выполнено 2 zen−m+k+1 − zen−m+k+1 kdkn k2Γe = k(zn+1 k − znk+1 )/∆k2Γe = = 0. ∆ Γ e Рассмотрим kdm n kΓ e. 2 kdm 2 n kΓ m 2 0 m−1 e = (zn+1 + c Azn + zn − (2 − A)znm )/∆ − ∆(F n (I({zk2 }n ))) Γe = 2 zen+1 + zen−1 − 2e zn 1 e = zn + c2 Ae + (Ae ezn−m ) − ∆(F n (I({e zk }n ))) = ∆ ∆ Γ e 2 zen+1 − 2e zn + zen−1 = ∆2 zn + c2 Ae + (Ae zn−m ) − (F n (I({e zk }n )) . ∆2 Γ e Отсюда по определению оператора A (7) и нормы (11)   2 zen+1 − 2e zn + zen−1 1 e + zn + c2 Ae (Ae ezn−m ) − F n (I({e zk }n )) 6 ∆2 ∆2 Γ e N −1 1 e X u(xi , tn+1 ) − 2u(xi , tn ) + u(xi , tn−1 ) u(xi+1 , tn ) − 2u(xi , tn ) + u(xi−1 , tn ) 6 2 k AkΓ 2 − a2 − ∆ ∆ h2 e i=1 u(xi+1 , tn−m ) − 2u(xi , tn−m ) + u(xi−1 , tn−m ) 2 − c2 a2 2 − Fni (I({uik }n )) h. (15) h Оценим каждое слагаемое под знаком суммы в (15), используя предположения теоремы: u(xi , tn+1 ) − 2u(xi , tn ) + u(xi , tn−1 ) u(xi+1 , tn ) − 2u(xi , tn ) + u(xi−1 , tn ) − a2 − ∆2 h2 u(xi+1 , tn−m ) − 2u(xi , tn−m ) + u(xi−1 , tn−m ) − c2 a2 ± Fni (utn (xi , ·)) 6 h2 6 |ψni | + |Fni (utn (xi , ·)) − Fni (I({uik }n ))| 6 6 C1 (∆p1 + hp2 ) + LF kutn (xi , ·) − I({uik }n )kQ 6 C1 (∆p1 + hp2 ) + LF C2 ∆p0 . (16) В результате из (15), (16) получаем N −1 1 e X kdm 2 n kΓ e 6∆ 2 kAkΓ (C1 (∆p1 + hp2 ) + LF C2 ∆p0 )2 h = kAk e e (N − 1)(C1 (∆p1 + hp2 ) + LF C2 ∆p0 )2 h 6 Γ ∆2 e i=1 339 6 (C3 ∆min{p1 ,p0 } + C4 hp2 )2 , q q где C3 = kAkΓe X (C1 + LF C2 ), C4 = kAk e e e X C1 . Γ То есть справедлива оценка kdn kΓ 6 C3 ∆min{p1 ,p0 } + C4 hp2 .  По определению устойчивости в общей разностной схеме необходимо, чтобы выполнялось условие: kSkΓ 6 1. К сожалению, в отличие от численных алгоритмов для уравнений параболического типа и уравнения переноса с аналогичными эффектами запаздывания в производных, эффективных критериев проверки устойчивости для данного алгоритма не найдено. Вложение в общую разностную схему с последействием [5, 6, 7, 8] проведено, откуда, используя [7, Теорема 1], получаем следующее утверждение. Теорема 3. Пусть схема устойчива, невязка в смысле (4) имеет порядок ∆p1 + hp2 , функции Fji липши- цевы, оператор интерполяции-экстраполяции I липшицев и имеет порядок погрешности p0 на точном решении, стартовые значения имеют порядок ∆p3 + hp4 , σ зафиксировано, тогда метод сходится с по- рядком ∆min{p0 ,p1 ,p3 } + hmin{p2 ,p4 } . Опираясь на теорему 1, получаем следующее следствие. Следствие. Если для точного решения задачи (1) существуют и непрерывны все частные производ- ные вплоть до 4-го порядка включительно, Fji определяется соотношением (5), применяется кусочно- линейная интерполяция (2), схема устойчива, σ зафиксировано, тогда метод сходится с порядком h2 + ∆2 . 4 Примеры численных расчетов Пример. Рассмотрим уравнение: ∂2u 2 2 π 2 a2 c2   2∂ u 2 2∂ u 2 2 (x, t) = a (x, t) + a c (x, t − 1) + 1 + π a + et sin πx − ∂t2 ∂x2 ∂x2 e √3 − e−t sin πx + u(x, t − τ (t)), 0 6 t 6 3, 0 6 x 6 1 (17) при τ (t) = 4t/3, a = 1, c = 0.5 с начальными: u(x, t) = et sin πx, −1 6 t 6 0, 06x61 и граничными условиями: u(0, t) = 0, u(1, t) = 0, 0 6 t 6 3. Уравнение имеет точное решение u(x, t) = et sin πx. На рис. 1 приведено приближенное решение этого уравнения методом (3), для которого Fji определяется соотношением (5), с кусочно-линейной интерполяцией (2) c числом точек разбиения по x, равным 15, по t – равным 60, на рис. 2 и 3 – разница между точным и приближенным решением. Рассмотренный показывает, что при увеличении количества шагов погрешность метода уменьшается. Это подтвержает выводы теоремы 3 и следствия из нее. Благодарности Исследования поддержаны Программой повышения конкурентоспособности ведущих университетов РФ (соглашение 02.А03.21.0006 от 27 августа 2013 г.) и проектом РНФ №14-35-00005. 340 u(x,t) 25 20 15 10 5 0 3 2 1 0.8 1 0.6 0.4 Time t 0.2 0 0 Distance x Рис. 1: Приближеное решение уравнения при N = 15, M = 60 0.02 0 -0.02 -0.04 -0.06 -0.08 1 3 2 0.5 1 0 0 Distance x Time t Рис. 2: Разница между точным и приближенным решением при N = 15, M = 60 Список литературы [1] L. Luo, Y. Wang. Oscillation for nonlinear hyperbolic equations with influence of impulse and delay. Int. J. Nonlinear Science, 14(1):60–64, 2012. [2] P. Garcia, M.A. Castro, J.A. Martin, A. Sirvent. Numerical solutions of diffusion mathematical models with delay. Mathematical and Computer Modelling, 50(5-6):860–868, 2009. [3] V.G. Pimenov. Numerical method for modeling controlled advection equation with delay in derivative with respect to time. Proceedings of international conference «System dynamics and control processes». Yekaterinburg: IMM UrO RAN. 2015, 272-278 (in Russian). = В.Г. Пименов. Численный метод моде- лирования управляемого уравнения переноса с запаздыванием в производной по времени. Труды международной конференции «Динамика систем и процессы управления (SDCP’14)». Екатеринбург: ИММ УрО РАН. 2015, 272-278. [4] V.G. Pimenov, E.E. Tashirova. Numerical methods for solving a hereditary equation of hyperbolic type. Proc. Steklov Inst. Math. (Suppl.), 281(suppl. 1):126—136, 2013. [5] A.V. Kim, V.G. Pimenov. i-Smooth Analysis and Numerical Methods for Solving Functional Differential Equations. Regulyarn. Khaotichesk. Dinamika, Izhevsk, 2004 (in Russian). = А.В. Ким, В. Г. Пименов i- Гладкий анализ и численные методы решения функционально-дифференциальных уравнений. Ижевск: НИЦ Регулярная и хаотическая динамика, 2004. 341 u(x,t) -3 x10 5 0 -5 -10 -15 -20 1 3 2 0.5 1 Distance x Time t 0 0 Рис. 3: Разница между точным и приближенным решением при N = 30, M = 120 [6] V.G. Pimenov. General Linear Methods for the Numerical Solution of Functional-Differential Equations. Differential Equations, 37(1):116–127, 2001. [7] V.G. Pimenov. Difference methods of solving partial differential equations with delay. Yekaterinburg: UrFU, 2014 (in Russian). = В.Г Пименов. Разностные методы решения уравнений в частных производных с наследственностью. Екатеринбург: УрФУ, 2014. [8] V.G. Pimenov, A. B. Lozhnikov Difference schemes for the numerical solution of the heat conduction equation with aftereffect. Proc. Steklov Inst. Math. (Suppl.), 275(suppl. 1):137—148, 2011. [9] A.A. Samarskii. The Theory of Difference Schemes. Nauka, Moscow, 1977; Marcel Dekker, New York, 2001. 342 Grid schemes for solving hyperbolic equation with delay in derivative Ekaterina E. Tashirova Ural Federal University (Yekaterinburg, Russia) Keywords: numerical methods, hyperbolic equation, time delay, order of convergence. A numerical method is constructed for solving equations of hyperbolic type with time delay in the second derivative with respect to space variable and functional delay in nonhomogeneity. The local error order of the algorithm is investigated. A theorem on the order of convergence of the algorithm is obtained by means of embedding it into a general scheme of difference methods for functional differential equations. Results of calculating a test example are presented. 343