On one method for solving the initial-boundary value problem for the Gardner equation

Full Text

1. Введение 1. Начально-краевая задача для уравнения Гарднера. Будем рассматривать начальнокраевую задачу для квазилинейного эволюционного уравнения в прямоугольнике x ∈ [0, L], t ∈ [0,T ], заданного нелинейным оператором Ξ[u](x, t), следующего вида: xxx Ξ[u] := ∂tu + (αu + β u2) ∂xu + γ ∂3 u = f (x, t), u = u(x, t), x ∈ [0, L], t ∈ [0,T ], (1.1) u(x, 0) = u0(x), (1.2) u(0, t)= ϕ0(t), u(L, t)= ψ0(t), ∂xu(L, t)= ψ1(t), (1.3) где α, β ∈ R,γ > 0. При β =0 уравнение (1.1) превращается в уравнение Кортевега-де-Фриза [5]. В общем виде уравнение (1.1) может быть сведено к уравнению Кортевега-де-Фриза при помощи преобразования Миуры [20, 21]. Вопросы разрешимости начально-краевых задач для уравнения Кортевега-де-Фриза изучены, например, в работе [17]. © С. И. Безродных, С. В. Пикулин, 2025 This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4.0/legalcode 353 Уравнение (1.1), называемое уравнением Гарднера, описывает волновые процессы в многокомпонентных стратифицированных жидкостях, в том числе внутренние волны в океане [15, 16], волны перепада уровня на мелкой воде - т. н. боры, или дисперсионные ударные волны [19], волны в конденсированной квантовой жидкости [14], а также ионно-звуковые волны в плазме [22, 24]. Отметим, что различные постановки задач для уравнения (1.1) привлекают внимание специалистов по численному анализу как содержательный пример для верификации различных вычислительных алгоритмов, см. [10], а также [11, 13, 23, 25]. Как правило, при численном решении нелинейных эволюционных уравнений отыскание искомой аппроксимации сводится к серии вспомогательных линейных задач (см., например, [8]). Поэтому при построении таких численных алгоритмов следует использовать адекватную решаемой задаче процедуру приближения нелинейного дифференциального оператора, а также опираться на достаточно точный и эффективный метод решения получаемых вспомогательных линейных задач. Настоящая работа является продолжением исследований [3, 4] по разработке вычислительных алгоритмов решения начально-краевых нелинейных параболических задач, применимых, в том числе, для уравнений типа Колмогорова-Петровского-Пискунова и Бюргерса. В разделе 2 развит аналитико-численный алгоритм решения начально-краевой задачи (1.1)-(1.3), в основе которого лежит сочетание явной экстраполяции по времени нелинейного слагаемого (αu + β u2) ∂xu в левой части уравнения (1.1) и неявной аппроксимации его линейных членов в соответствии с модифицированной схемой Кранка-Николсон [12]. Новизна разработанного в разделе 2 численного алгоритма для уравнения (1.1) заключается в предложенном методе решения вспомогательных линейных задач. Решения этих задач мы строим в полуаналитическом виде с использованием явного вида частных решений линейного уравнения и его фундаментальной системы решений. Применяемый метод позволяет эффективно решать возникающие линейные краевые задачи с алгоритмической сложностью O(K), где K - число отрезков дискретизации по пространственной переменной. Отметим также, что близкий подход использовался в работе [2] при решении сингулярно возмущенной краевой задачи, возникающей в теории полупроводниковых приборов [1]. В разделе 3 предложенный метод проиллюстрирован с помощью численной реализации для начально-краевой задачи о распространении солитона «большой амплитуды» (см. [7]) при отрицательном значении коэффициента β в уравнении (1.1). Полученные при сопоставлении численного результата с явным аналитическим решением данные подтверждают высокую точность и эффективность предложенного метода. 2. Общая схема вычислительного алгоритма. Зададим равномерную дискретизацию временн´ого отрезка [0,T ] с шагом τ = T /N, N ∈ N и введем следующие обозначения: un(x) ≈ u(x, tn), x ∈ [0, L], tn = nτ, n = 0,... , N, ( ) d Fn(x) := α un(x)+ β u2 (x) un(x), (1.4) n dx где un(x) - искомое приближенное решение задачи (1.1)-(1.3) в момент времени t = tn, причем функция u0(x) задана начальным условием (1.2). Будем вычислять приближения un(x) последовательно при n = 1,... , N. Для получения функции u1(x) из u0(x) можно воспользоваться любой двухслойной вычислительной схемой, например, явно-неявной схемой на основе приближения Кранка-Николсон [8]. В настоящей работе при сопоставлении численных результатов с явным решением в разделе 3 функция u1(x) и ее производные по x задаются в качестве известных величин наряду с начальными данными (1.2). Предполагая, что приближенные значения u(x, t) известны на первых n слоях, т. е. заданы функции u0(x), u1(x), ... , un(x), найдем приближение на (n + 1)-м слое, т. е. функцию un+1(x). Для этого будем использовать следующую аппроксимацию нелинейного дифференциального оператора Ξ[u] из (1.1), центрированную в точке t = tn + τ /2, на основе трехслойной явно-неявной схемы [12, 18]: xxx Ξ[u]= ∂tu + (αu + β u2) ∂xu + γ ∂3 u ≈ un+1(x) - un(x) τ 3 + 2 Fn(x) - 1 2 Fn -1(x)+ (1.5) 9 d3 + 16 γ dx3 3 un+1(x)+ 8 γ d3 dx3 1 un(x)+ 16 d3 γ dx3 un-1(x); здесь для нелинейного члена применена явная экстраполяция Адамса-Бэшфорта, а линейное слагаемое со старшей производной приближено по модифицированной схеме Кранка- Николсон [12]. Подставляя приближение (1.5) в уравнение (1.1), получаем уравнение относительно значений искомой аппроксимации un+1(x) на (n + 1)-м слое: 9 d3 16 τγ dx3 un+1(x)+ un+1(x)= gn+1(x), (1.6) где правая часть определяется по формуле 3 d3 1 d3 3 1 gn+1(x)= un(x) - 8 τγ dx3 un - 16 τγ dx3 un-1 - τ ( 2 Fn(x) - 2 Fn-1(x)). (1.7) Используя уравнения вида (1.6), отвечающие индексам n, n - 1, перепишем выражение (1.7) следующим образом: 2 1 3 1 gn+1(x)= un(x)+ 3 (un(x) - gn(x)) + 9 (un-1(x) - gn-1(x)) + τ ( 2 Fn(x) - 2 Fn-1(x)) = 5 1 = 3 un(x)+ 9 un-1(x) - 2 3 gn(x) - 1 3 9 gn-1(x)+ τ ( 2 Fn(x) - 1 2 Fn-1(x)). (1.8) Отметим, что правая часть (1.7) уравнения (1.6) зависит только от функций, уже вычисленных на предыдущих шагах алгоритма. Кроме того, форма (1.8) представления gn+1(x) исключает явную зависимость этой правой части от старших производных функций uj (x), j = 0,... ,n - 1, хотя их первые производные входят в выражения (1.4) для нелинейных членов Fn(x), Fn-1(x). Уравнение (1.6), (1.8) дополняется краевыми условиями un+1(0) = ϕ0(tn+1), un+1(L)= ψ0(tn+1), в соответствии с условиями (1.3) исходной задачи. du (L)= ψ1(tn+1) (1.9) dx Таким образом, решение исходной нелинейной задачи (1.1)-(1.3) сведено к последовательному решению линейной задачи (1.6), (1.8), (1.9) с изменяющимися на каждом шаге алгоритма правой частью gn+1(x) и граничными значениями в (1.9). В следующем разделе 2 представлен численно-аналитический метод, позволяющий эффективно вычислять решение основной линейной задачи (1.6), (1.8), (1.9). 2. Численно-аналитический метод решения основной линейной краевой задачи Для удобства изложения предлагаемого метода решения краевой задачи (1.6), (1.8), (1.9) перепишем уравнение (1.6) в новых обозначениях: 1 d3y L [y](x) ≡ ω3 dx3 + y(x)= g(x), ω > 0, x ∈ [0, L], (2.1) а краевые условия (1.9) - в виде y(0) = Φ0, y(L)= Ψ0. dy (L)= Ψ1. (2.2) dx Нетрудно убедиться в том, что следующий набор функций составляет фундаментальную систему решений уравнения (2.1) F1(x) := eμx, F2(x) := eρx, F3(x) := eρx (2.3) где μ, ρ, ρ - корни характеристического уравнения λ3/ω3 +1 = 0: λ1 =: μ = -ω, λ2 =: ρ = ω β, λ3 =: ρ = ω β-1, β := eπi/3 = √3 -1; (2.4) при этом F1(x) является вещественной, а F2(x) и F3(x) - комплексными функциями. Нам потребуются следующие соотношения между величинами μ и ρ: ρ ρ = ω2, ρ2 = -ωρ, ρ 2 = -ωρ, ρ - ρ = ω √3 i, ρ - μ = ω √3 eπi/6 = i √3 ρ. (2.5) μ + ρ + ρ = ω (-1+ β + β)=0 = ω2 (1 - β - β)= μ2 + ρ2 + ρ 2. (2.6) Введем на отрезке [0, L] сетку, образованную узлами 0= ξ0 < ξ1 < ··· < ξK = L (2.7) и обозначим через Pm пространство многочленов от переменной x степени не выше m 0 (dim Pm = m + 1). Будем предполагать, что правая часть g(x) уравнения (2.1) является кусочно полиномиальной функцией, совпадающей на каждом из интервалов [ξk-1, ξk ] с некоторым полиномом fk(x), т. е. справедливы соотношения g(x) [ξk-1,ξk ] =: fk(x) ∈ Pm, k = 1,... , K. (2.8) Функция (2.8), вообще говоря, может иметь разрывы в точках ξk. Решение y(x) уравнения (2.1) с такой правой частью будем искать в классе C2[0, L], понимая решение уравнения в смысле выполнения интегрального тождества, эквивалентного (2.1) (см. [9]). Предлагаемый численный алгоритм решения задачи (2.1), (2.2) состоит из двух этапов. На первом этапе (см. п. 2.1) построим частные решения уравнения (2.1) отдельно на каждом из промежутков пространственного разбиения (2.7). На втором этапе (см. п. 2.2) с помощью этих частных решений получим решение краевой задачи (2.1), (2.2) на всем отрезке интегрирования x ∈ [0, L]. 1. Построение частных решений на отрезках разбиения. Построим частное решение уравнения (2.1) с правой частью вида (2.8) на k-м отрезке разбиения (2.7) при помощи интегрального представления ξk r k y (x)= fk(ξ) E(x, ξ) dξ, x ∈ [ξk-1, ξk], k = 1,... , K, (2.9) ξk-1 где функция fk (ξ) задана равенством (2.8), а через E(x, ξ) обозначено фундаментальное решение уравнения (2.1), т. е. решение уравнения L E(x, ξ) = δ(x - ξ); (2.10) здесь ξ ∈ R - параметр, а δ( · ) - функция Дирака. Утверждение 2.1. Функция 1 ( -μ eμ (x-ξ), ξ x, E(x, ξ)= 3 ρ eρ (x-ξ) + ρ eρ (x-ξ), ξ > x (2.11) является фундаментальным решением уравнения (2.1), т. е. удовлетворяет равенству (2.10) при всех ξ ∈ R. Отметим, что функция (2.11) является вещественной и экспоненциально убывает с ростом расстояния |x - ξ|. Доказательство. Фундаментальное решение (2.11) может быть получено следующим образом (см. [4, 6]). Образуем из функций фундаментальной системы решений (2.3) и их производных матрицу ⎛ D(x) := ⎝ eμx eρx eρx μ eμx ρ eρx ρ eρx μ2 eμx ρ2 eρx ρ 2 eρx ⎞ ⎠ = M × diag(eμx, eρx, eρx), (2.12) где постоянная матрица M имеет вид ⎛ 1 1 1 ⎞ ⎛ 1 1 1 ⎞ M := ⎝ μ ρ ρ μ2 ρ2 ρ 2 ⎠ = diag( 1, ω, ω2 ) × ⎝ -1 β β ⎠ . (2.13) 1 -β -β Первую строку матрицы D(x) составляют функции (2.3), вторую и третью строки - соответственно первые и вторые производные этих функций. Здесь и далее diag( ·, ·, · ) означает диагональную матрицу с указанными элементами на главной диагонали. Обратим матрицу (2.12) и запишем результат, заменяя независимую переменную x на ξ: D-1(ξ)= diag( e-μξ, e-ρξ , e-ρξ ) × M-1, где M-1 = 1 ⎛ 1 -1 1 ⎞ 1 β -β × diag( 1, ω-1, ω-2 ). (2.14) Тогда произведение - 3 ⎝ 1 β β ⎠ G(x, ξ) := D(x) × D-1(ξ)= M × diag(eμ (x-ξ), eρ (x-ξ), eρ (x-ξ) \ × M-1 (2.15) обращается в единичную матрицу E при x = ξ, т. е. G(x, ξ) x=ξ = E, x, ξ ∈ R, (2.16) причем вторая и третья строки матрицы-функции G(x, ξ) являются в силу построения (2.15) соответственно первой и второй производной по x от ее первой строки, поскольку тем же свойством обладает матрица-функция G(x). Из (2.16) следует, что матричный элемент G (3) 1 μ (x-ξ) ρ (x-ξ) ρ (x-ξ)) (1) (x, ξ)= 3 ω2 (e § βe § βe , (2.17) занимающий правый верхний угол матрицы (2.15), удовлетворяет равенствам 2 (1) G (3) d (3) d (3) (x, ξ) = x=ξ (1) G (x, ξ) dx x=ξ = 0, dx2 G(1) (x, ξ) x=ξ = 1, откуда получаем, что функция (1) E0(x, ξ) := ω3 G (3)(x, ξ) θ(x - ξ) (2.18) является фундаментальным решением уравнения (2.1); здесь θ( · ) - функция Хевисайда. Покажем, что функция (2.11) отличается от E0(x, ξ) прибавлением линейной комбинации функций F2 и F3, принадлежащих фундаментальной системе (2.3) и, следовательно, также удовлетворяет условию (2.10). В самом деле, из (2.17), (2.18) с учетом равенств (2.5) имеем ωβ ω β ω ρ (x-ξ) ρ (x-ξ)) E0(x, ξ)+ 3 F2 + 3 F3 = E0(x, ξ)+ 3 (βe + βe = Утверждение доказано. ω ( eμ (x-ξ), x ξ = 3 β eρ (x-ξ) + β eρ (x-ξ), x < ξ = E(x, ξ). тем Частное решение yk(x) уравнения (2.1) на отрезке [ξk-1, ξk] может быть вычислено по формулам (2.9), (2.11) пу численного интегрирования на основе квадратурных формул. Для случая кусочно полиномиальной правой части g(x) этого уравнения, имеющей вид (2.8), проведем такое вычисление в явном виде. Подставляя фундаментальное решение E(x, ξ) в виде (2.11) в интегральное представление (2.9), находим ξk r k y (x)= fk(ξ) E(x, ξ) dξ = ξk-1 1 ( = 3 -μ x r ξk-1 eμ (x-ξ) fk (ξ) dξ + ρ ξk ξk r r eρ (x-ξ) fk(ξ) dξ + ρ x x eρ (x-ξ) fk(ξ) dξ\ = ξk-1 1 ( r = μ eμ (x-ξ) fk(ξ) dξ + ρ 3 x ξk ξk r r eρ (x-ξ) fk(ξ) dξ + ρ x x eρ (x-ξ) fk (ξ) dξ\. (2.19) Для того, чтобы вычислить интегралы, стоящие в правой части равенства (2.19), рассмотрим следующую более общую ситуацию. Предположим, что заданы некоторый многочлен f (x) ∈ Pm степени m 0 и число ν ∈ C, ν ±= 0, и найдем первообразную функции e-νx f (x). Применяя несколько раз интегрирование по частям, получим r -νx m j -νx f m - e-νx f (x) dx = e ν \ 1 νj j=0 d dxj e f (x)+ C = - ν D(m))j \( 1 ν j=0 f (x)+ C, (2.20) где через 1m+1 обозначен тождественный оператор в пространстве Pm, а через Dm - оператор дифференцирования по x, действующий в этом же линейном пространстве, т. е. (m) df df D(m): Pm → Pm, f D , f (x), ∈ Pm. (2.21) ⊕-→ dx dx Отметим, что оператор (2.21) является нильпотентным: dm+1f m+1 ∀f ∈ Pm : dξm+1 =0 =⇒ (D(m)) = 0. (2.22) Используя свойство (2.22), перепишем формулу (2.20) в виде r -νx f -1 - e-νx f (x) dx = e ν 1m+1 1 - ν D(m) f (x)+ C, C = const. (2.23) Отметим, что входящее в правую часть формулы (2.23) выражение F (x; ν)[ f ] := f 1 -1 - D(m) f (x) (2.24) 1m ν также задает многочлен из пространства Pm. Вычислим, исходя из (2.23), определенный интеграл от функции e-νξ f (ξ) в пределах от x до Ξ, где x, Ξ ∈ R - некоторые числа, пользуясь формулой Ньютона-Лейбница: Ξ r -νξ lξ=Ξ ( \ - e-νξ f (ξ) dξ = r e ν x 1 ( F (ξ; ν) 1 = ξ=x ν e-νx F (x; ν) - e-νΞ F (Ξ; ν) = \ = ν eνx F (x; ν) - F (Ξ; ν) eν(x-Ξ) ; (2.25) здесь в выражениях F ( ·; ν)[ f ] для простоты опущен параметр f. Введем следующее обозначение: J (x; ν, Ξ)[ f ] := ν eνx Тогда с учетом формулы (2.25) имеем Ξ r e-νξ f (ξ) dξ. (2.26) x J (x; ν, Ξ)[ f ]= F (x; ν)[ f ] - F (Ξ; ν)[ f ] eν(x-Ξ). (2.27) Введем также обозначения для выражений (2.24) и (2.26) при f = fk(ξ): f 1 -1 Fk (x; ν) := F (x; ν)[fk ] ≡ 1m+1 - ν D(m) fk(x), Fk (x; ν) ∈ Pm, (2.28) Jk(x; ν, Ξ) := J (x; ν, Ξ)[fk]= Fk (x; ν) - Fk (Ξ; ν) eν(x-Ξ). (2.29) Вычислим интегралы в правой части равенства (2.19) по формуле (2.25) при f = fk и ν = μ или ν = ρ. Используя введенные обозначения (2.24) (2.26), (2.28), (2.29), получим 1 k y (x)= 3 (Jk(x; μ, ξk-1)+ Jk(x; ρ, ξk )+ Jk(x; ρ, ξk )), x ∈ [ξk-1, ξk ]. (2.30) Найдем первую и вторую производные полученного частного решения (2.30) уравнения (2.1) на отрезке [ξk-1, ξk ]. Для этого сначала продифференцируем по x выражение (2.26): Ξ d 2 dx J (x; ν, Ξ)[ f ]= ν eνx r x e-νξ f (ξ) dξ - ν eνx e-νx f (x)= ν J (x; ν, Ξ)[ f ] - ν f (x). (2.31) Дифференцируя по x равенство (2.30), с учетом формул (2.31) и (2.6) находим d 1 yk(x)= ( d Jk(x; μ, ξk-1)+ d Jk (x; ρ, ξk )+ ) d Jk(x; ρ, ξk ) = dx 3 dx 1 dx dx 1 = 3 (μ Jk (x; μ, ξk-1)+ ρ Jk(x; ρ, ξk )+ ρ Jk(x; ρ, ξk )) - 1 3 (μ + ρ + ρ) fk(x)= = 3 (μ Jk (x; μ, ξk-1)+ ρ Jk(x; ρ, ξk )+ ρ Jk(x; ρ, ξk )). (2.32) Повторно дифференцируя по x полученное выражение (2.32), аналогичным образом находим d2 1 yk(x)= (μ2 Jk (x; μ, ξk 1)+ ρ2 J (x; ρ, ξ )+ ρ2 J (x; ρ, ξ )) - 1 (μ2 + ρ2 + ρ2) f (x)= dx2 3 1 - k k k k 3 k 3 = (μ2 Jk (x; μ, ξk - 1)+ ρ2 Jk (x; ρ, ξk )+ ρ2 Jk (x; ρ, ξk )). (2.33) Введем в рассмотрение вектор-функцию dy y d2 \T yk yk (x) := ( (x), k (x), k (x) ; (2.34) dx dx2 здесь и далее ( · )T означает матричное транспонирование. Тогда равенства (2.30), (2.32), (2.33) принимают вид одного векторного соотношения 1 f 2 )T 2 )T ( 2 )T k y (x)= 3 Jk(x; μ, ξk-1)( 1, μ, μ + Jk(x; ρ, ξk )( 1, ρ, ρ + Jk(x; ρ, ξk ) 1, ρ, ρ = 1 T = 3 M × ( Jk(x; μ, ξk-1), Jk(x; ρ, ξk ), Jk(x; ρ, ξk ) ) , x ∈ [ξk-1, ξk ], (2.35) где матрица M имеет вид (2.13). Таким образом, искомое частное решение (2.9) уравнения (2.1) с полиномиальной правой частью на отрезке [ξk-1, ξk ] построено в виде (2.30) и, вместе с первыми двумя производными, - в векторном виде (2.35). Для дальнейшего изложения вычислительного алгоритма нам потребуется найти значения вектора (2.35) в концах отрезка [ξk-1, ξk ]. Для того, чтобы их получить, воспользуемся тождеством J (Ξ; ν, Ξ)[ f ]= 0, Ξ ∈ R, ν ±= 0, f ∈ Pm, (2.36) которое непосредственно вытекает из определения (2.26). Подставляя последовательно x = ξk-1 и x = ξk в формулу (2.35), с учетом (2.36) находим 1 yk (ξk )= 3 T M × ( Jk(ξk ; μ, ξk-1), 0, 0 ) , k = 1,... , K, (2.37) yk (ξk-1 × ( J 1 )= M 0, k 3 (ξk-1 ; ρ, ξk ), Jk (ξk-1 ; ρ, ξk ) )T. Последнее равенство запишем, производя подстановку k ⊕→ (k + 1), в следующем виде: yk+1 (ξk × ( J 1 )= M 0, 3 k+1 (ξk ; ρ, ξ k+1 ), J k+1 (ξk ; ρ, ξ k+1 ) )T . k = 0,... ,K - 1. (2.38) Таким образом, приходим к следующему результату. Утверждение 2.2. Уравнение (2.1) с кусочно-полиномиальной правой частью (2.8) имеет на каждом отрезке [ξk-1, ξk ], k = 1,... , K, частное решение yk(x) вида (2.30), где выражения Jk (x; ·, · ) даны формулами (2.29), (2.28). Вектор (2.34), составленный из производных функции yk (x), может быть найден по формуле (2.35). Значения таких векторов в узлах ξk находятся по формулам (2.37), (2.38). 2. Решение линейной краевой задачи на отрезке. Построим решение краевой задачи (2.1), (2.2), (2.8), предполагая, что на отрезках [ξk-1, ξk ], разбиения (2.7) заданы какие-либо вещественнозначные частные решения yk(x) ∈ C2[ξk -1, ξk] уравнения (2.1), (2.8): L [ yk ](x)= fk(x), x ∈ [xk-1, xk ], k = 1,... , K. (2.39) k Для этих частных решений мы сохраняем обозначение y (x) из предыдущего пункта, не полагая k при этом, что они обязательно заданы формулой (2.9). Как и выше, будем обозначать через y (x) вектор-функцию (2.34). Определим для каждого отрезка [ξk-1, ξk ], k = 1,... , K, разбиения (2.7) фундаментальную систему решений уравнения (2.1), элементы которой запишем в виде строки dk (x) := (eμ (x-ξk-1), eρ(x-ξk ), eρ(x-ξk )). (2.40) Дополняя строку (2.40) ее первой и второй производной по x, образуем по аналогии с матрицей (2.12) матрицу ⎛ Dk (x) := ⎝ eμ (x-ξk-1 ) eρ(x-ξk ) eρ(x-ξk ) μ eμ (x-ξk-1 ) ρ eρ(x-ξk ) ρ eρ(x-ξk ) μ2 eμ (x-ξk-1 ) ρ2 eρ(x-ξk ) ρ 2 eρ(x-ξk ) ⎞ ⎠ = (2.41) = M × diag(eμ (x-ξk-1), eρ(x-ξk ), eρ(x-ξk )), (2.42) где постоянная матрица M дана формулой (2.13); здесь и далее знак × обозначает матричное умножение. y(x) искомое решение краевой задачи (2.1), (2.2), (2.8), а через (x) - вектор Обозначим через y y(x) по аналогии с (2.34). На каждом из отрезков [ξk-1, ξk] это производных порядка 0, 1, 2 от решение имеет вид где yk(x)+ dk (x) × hk, (2.43) yk (x)+ Dk (x) × hk, x ∈ [xk-1, xk ], k = 1,... , K, (2.44) hk ≡ (ak, bk , ck )T, k = 1,... , K, (2.45) - некоторый постоянный вектор-столбец коэффициентов, элементы которого будут найдены далее, а строка dk (x) и матрица Dk (x) заданы равенствами (2.40), (2.41). ∈ Условие гладкости y(x) C2[0, L] эквивалентно требованию непрерывности вектор-функции лах k y(x) во внутренних уз ξ отрезка [0, L]. В силу представления (2.44) это приводит к следующим условиям: yk (ξk)+ Dk (ξk) × hk = (ξ )+ D (ξ ) × h , k = 1,... ,K - 1. (2.46) yk+1 k k+1 k k+1 Соотношения (2.46) представляют собой систему (K - 1) линейных алгебраических уравнений относительно K векторов (2.45). Потребуем в дополнение к этому выполнения краевых условий (2.2), что обеспечит совпадение числа уравнений и числа неизвестных. Для того чтобы сформулировать указанное требование, введем фундаментальные системы решений d(0)(x) и d(L)(x) уравнения (2.1), связанные с концевыми точками отрезка [0, L] x =0= ξ0 и x = L = ξK соответственно: d(0)(x) := ( (x), F (x), F (x)), d (x) := d (x - L), (2.47) F1 2 3 (L) 0 где функции Fj (x), j = 1, 2, 3, являются экспонентами (2.3). Определим также соответствующие матрицы, составленные из производных фундаментальных систем (2.47): D(0)(x) := D(x)= D1(x) × diag( 1, e-ρ ξ1 , e-ρ ξ1 ), (2.48) D(L)(x) := D(x - L)= DK (x) × diag( eμ (L-ξk-1 ), 1, 1 ), (2.49) где матрица D(x) имеет вид (2.12), а D1(x) и DK (x) даны формулой (2.41) при k =1 и k = K. Запишем равенство (2.44) при k = 1, x = ξ0 =0 с учетом (2.48) в следующем виде: y1(0) = D1(0) × h1 = D(0)(0) × diag( 1, e- ρ ξ1 , e-ρ ξ1 ) × h1 = = M × ( a1, e-ρ ξ1 b1, e-ρ ξ1 c1 )T. (2.50) Аналогичным образом равенство (2.44) при k = K, x = ξK = L запишем в виде yK (L)= DK (L) × hK = D(L)(L) × diag( eμ (L-ξk-1 ) , 1, 1 ) × hK = T = M × ( eμ (L-ξK-1 ) aK, bK , cK ) . (2.51) Рассмотрим векторные равенства (2.50), (2.51) покомпонентно и выпишем те три уравнения из получившихся шести, в которые входят данные краевых условий (2.2): Φ0 - y1(0) = a1 + e-ρ ξ1 b1 + e-ρ ξ1 c1, (2.52) Ψ0 - yK yK (L)= eμ (L-ξK-1 ) aK + bK + cK, (2.53) 1 - Ψ (L)= μ eμ (L-ξK-1 ) aK + ρ bK + ρ cK , (2.54) dx где уравнение (2.52) отвечает первой компоненте равенства (2.50), а уравнение (2.53), (2.54) - первым двум компонентам равенства (2.51). Уравнения (2.46), (2.52)-(2.54) образуют систему линейных алгебраических уравнений, в которой число уравнений 3K совпадает с числом искомых скалярных величин (2.45). Для того, чтобы сформулировать результат о разрешимости этой системы и дать эффективные формулы вычисления коэффициентов ak, bk , ck , введем ряд обозначений. Определим числа uk, vk , wk , k = 1,... ,K - 1, характеризующие размер скачков частных решений yk(x) и их производных во внутренних узлах сетки (2.7), следующим образом: (uk, vk, wk )T ≡ zk := M-1 × (yk+1(ξk ) - yk(ξk )), k = 1,... ,K - 1, (2.55) где матрица M-1 имеет вид (2.14). Рассмотрим также следующие линейные комбинации этих чисел: K-1 K-1 K-1 U := \ eμ (L-ξk ) uk, V := \ e-ρ ξk vk, W := \ e-ρ ξk wk. (2.56) k=1 k=1 k=1 Определим величины Φ 0, Ψ 0, Φ 1 и R1, R2, связанные с краевыми условиями (2.2): yK y1(0), 1 Ψ 0 := Ψ0 - yK (L), μL dx Ψ 1 := Ψ1 - (L), (2.57) R1 := Ψ 0 + ω Ψ 1, R2 := Ψ 0 + U - e (Φ 0 - V - W ). (2.58) yk (x) и функций, Отметим, что из вида матрицы (2.14) и вещественности компонент векторов задающих краевые условия (2.2), вытекают следующие соотношения: uk ∈ R, wk = vk, k = 1,... ,K - 1; U, Φ 0, Ψ 0, Φ 1, R1, R2 ∈ R, W = V. (2.59) Введем в рассмотрение число B, зависящее только от произведения (ωL): B := Im (β + 1) (1 - e-( +1) ωL) , (2.60) где β определено последней из формул (2.4). Заметим, что поскольку Re{β + 1} = 3/2 > 0, то при достаточно больших значениях (ωL) величина B отлична от нуля. Теорема 2.1. Предположим, что величина B, определяемая по формуле (2.60), отлична от нуля. Тогда система уравнений (2.45), (2.46), (2.52)-(2.54) однозначным образом разрешима относительно 3K искомых величин ak, bk, ck , k = 1,... ,K - коэффициентов представления (2.43)-(2.45) - при любых значениях параметров Φ0, Φ1, Ψ0 ∈ R и любых заданных функциях yk (x) ∈ C2[ξk -1, ξk ], k = 1,... , K. Справедливы следующие формулы для вычисления коэффициентов bK и a1: (1 - e-( +1) ωL) R1 - (β + 1) R2 bK = , (2.61) 2i B a1 = Φ 0 - 2 Re{V + e-ρL bK }, (2.62) где Φ 0 определено формулой (2.57). Остальные коэффициенты из серий {ak }, {bk } определяются рекуррентным способом: ak+1 = eμ (ξk-ξk-1) ak - uk, k = 1,... ,K - 1, (2.63) bk = e-ρ (ξk+1-ξk ) bk+1 + vk, k = K - 1,... , 1. (2.64) Для коэффициентов {ck } справедливо соотношение ck = bk, k = 1,... , K. (2.65) Доказательство. Умножая равенство (2.46) слева на матрицу M-1, перепишем его с учетом формулы (2.42) и обозначения (2.55) в виде diag( eμ (ξk-ξk-1), 1, 1, ) × hk - diag(1, eρ (ξk-ξk+1), eρ (ξk-ξk+1)) × hk+1 = = M-1 × (yk+1(ξk ) - yk (ξk)) = zk, k = 1,... ,K - 1. (2.66) Покомпонентная запись векторного равенства (2.66) приводит к соотношениям (2.63), (2.64) и следующей формуле: ck = ck+1 e-ρ (ξk+1 -ξk ) + wk, k = 1,... ,K - 1. (2.67) Перейдем к выводу равенств (2.61), (2.62). Последовательно применяя K-1 раз формулу (2.63), выразим через a1 величину eμ (L-ξK-1 ) aK, входящую в правые части уравнений (2.53) и (2.54): eμ (L-ξK-1 ) aK = eμ (L-ξK-1 ) (eμ (ξK-1 -ξK-2 ) aK -1 - uK-1) = = eμ (L-ξK-2 ) (eμ (ξK-2 -ξK-3 ) aK -2 - uK-2) - - eμ (L-ξK-1 ) uK 1 = ··· = = eμ (L-ξ1 ) (eμ (ξ1 -ξ0) a1 - u1) - eμ (L-ξ2 ) u2 - ··· - eμ (L-ξK-2 ) uK - - 2 - eμ (L-ξK-1 ) uK 1 = = eμL a1 - U, (2.68) где U задано первым из равенств (2.56). Аналогичным образом выразим через bK и cK слагаемые в правой части уравнения (2.52), содержащие b1 и c1, используя рекуррентные соотношения (2.64) и (2.67) соответственно: e-ρ ξ1 b1 = e-ρ ξ1 (e-ρ (ξ2 -ξ1 ) b2 + v1) = ··· = e-ρL bK + V, (2.69) e-ρ ξ1 c1 = e-ρ ξ1 (e-ρ (ξ2 -ξ1 ) c2 + w1) = ··· = e-ρL cK + W, (2.70) где величины V и W также заданы формулами (2.56). Подставляя выражения (2.68)-(2.70) в уравнения (2.52)-(2.54), исключим величины aK, b1 и c1 из этой системы. Тогда с использованием обозначений (2.57) получаем следующую систему трех линейных алгебраических уравнений относительно коэффициентов a1, bK и cK : a1 + (e-ρL bK + V ) + (e-ρL cK + W ) = Φ 0, (2.71) (eμL a1 - U ) + bK + cK = Ψ 0, (2.72) μ (eμL a1 - U ) + ρ bK + ρ cK = Ψ 1, (2.73) Исключим a1 из уравнений (2.72), (2.73). Для этого умножим первое из них (2.72) на μ и вычтем результат из второго (2.73), получим (ρ - μ) bK + (ρ - μ) cK = Ψ 1 - μ Ψ 0. (2.74) Затем, подставляя в уравнение (2.72) коэффициент a1, выраженный из уравнения (2.71) в виде a1 = Φ 0 - V - W - e-ρL bK - e-ρL cK, (2.75) получаем (1+ e(μ-ρ) L) bK + (1+ e(μ-ρ) L) cK = Ψ 0 + U - eμL (Φ 0 - V - W ). (2.76) Запишем уравнение (2.74), разделив его на ω = -μ, вместе с уравнением (2.76) с учетом равенств (2.4), (2.5) в следующем векторном виде: β +1 β +1 1 - e-( +1) ωL 1 - e-( +1) ωL bK × cK = R1 R2 , (2.77) где величины R1, R2 заданы формулами (2.58). Решая полученную систему (2.77), находим (1 - e-( +1) ωL) R1 - (β + 1) R2 -(1 - e-( +1) ωL) R1 + (β + 1) R2 bK = , cK = , 2i Im (β + 1) (1 - e-( +1) ωL) 2i Im (β + 1) (1 - e-( +1) ωL) откуда с учетом обозначения (2.60) получаем формулу (2.61), а также соотношение cK = bK. (2.78) Формула (2.62) вытекает из (2.75), (2.78) и (2.59), а формула (2.65) следует по индукции из (2.78), (2.67), (2.64) и (2.59). Теорема доказана. k 3. Вычислительный алгоритм решения линейной краевой задачи. Дадим уточнение формул (2.55) для случая, когда частные решения y (x) уравнения (2.1), (2.8) заданы в виде (2.30), (2.35). Подставляя равенства (2.37), (2.38) в формулу (2.55) при k = 1,... ,K - 1, находим 1 zk = M-1 × (yk+1(ξk ) - yk (ξk )) = M-1 × ( M × ( 0, Jk+1(ξk ; ρ, ξk+1), Jk+1(ξk ; ρ, ξk+1) )T- 3 1 T\ 1 T = - 3 M × ( Jk(ξk ; μ, ξk-1), 0, 0 ) ( 3 -Jk (ξk ; μ, ξ k-1 ), J k+1 (ξk ; ρ, ξ k+1 ), J k+1 (ξk ; ρ, ξ k+1) ) , т. е. для коэффициентов {uk }, {vk }, {wk } имеем следующие выражения: 1 1 uk = - 3 Jk (ξk; μ, ξk-1), vk = 3 Jk+1(ξk ; ρ, ξk+1), wk = vk, k = 1,... ,K - 1, (2.79) где величины Jk(x; ·, · ) определяются по формулам (2.29), (2.28). Остановимся подробнее на вычислении полиномов Fk (x; ν) вида (2.28), которые составляют основу построения искомого решения задачи (2.1), (2.2), (2.8) согласно теореме 2.1 и формулам (2.29), (2.79). Предположим, что многочлен fk(x) ∈ Pm, задающий правую часть уравнения (2.1) на k-м отрезке разбиения (2.7), записан в виде линейной комбинации m j fk(x)= \ fk,j T (k)(x), fk,j ∈ R, (2.80) базисных векторов j=0 T (k) ( j (x) := Tj 2 x - ξk-1 Δk 1 - \, j = 0,... , m, (2.81) где Tj (X) - j-й многочлен Чебышева первого рода, X ∈ [-1, 1], а через Δk обозначена длина рассматриваемого отрезка: Δk := ξk - ξk-1, k = 1,... , K. (2.82) Обозначим через T(m) матрицу оператора дифференцирования многочленов от X в базисе m (k) m {Tj (X)}j=0. Тогда матрица оператора дифференцирования D(m) в базисе {Tj (x)}j=0 равна 2 T(m)/Δk, а матрица оператора (1 - D(m)/ν), стоящего в правой части определения (2.28), имеет вид 2 k R(k)(m, ν) := Em+1 - ν Δ T(m), (2.83) где Em+1 = diag(1,... , 1) - единичная матрица порядка m + 1. Матрица T(m) является верхнетреугольной и имеет нули на главной диагонали. Следовательно, матрица R(k)(m, ν), как и обратная к ней, также является верхне-треугольной и имеет единицы на главной диагонали. Вычисляя эту обратную матрицу и применяя ее к вектору (2.80), получаем вектор-столбец 1 (R(k)(m, ν))- T × (fk,0,... , fk,m) коэффициентов разложения искомого многочлена Fk (x; ν) по базису (2.81). Сформулируем построенный вычислительный алгоритм решения задачи (2.1), (2.2), (2.8) с кусочно полиномиальной правой частью. Алгоритм 2.1. 1. Предполагаем, что многочлены (2.8), задающие правую часть уравнения (2.1), представлены коэффициентами fk,j разложения (2.80) по базису (2.81), j = 0,... , m, k = 1,... , K. 1 1 2. Вычисляем матрицы (R(k)(μ, m))- и (R(k)(ρ, m))- , k = 1,... , K, обратные к матрицам (2.83) при ν = μ и ν = ρ соответственно. j=0 3. Умножая векторы (fk,j )m слева на найденные матрицы для каждого значения индекса k = 1,... , K, вычисляем по формуле (2.28) многочлены Fk (x; μ) и Fk (x; ρ) в виде разложений по базису (2.81). 4. Подставляя в формулу (2.29) последовательно x = ξk при k = 1,... ,K и ν = ρ, ν = μ, находим значения Jk+1(ξk ; μ, ξk+1), Jk(ξk ; ρ, ξk ),k = 1,... ,K - 1, из которых по формулам (2.79) получаем числа {uk }, {vk = wk }. 5. По формулам (2.56)-(2.58), (2.60)-(2.65) находим коэффициенты {ak }, {bk = ck }. 6. Получаем искомое решение y = y(x) вместе с его первыми двумя производными в ви- k де (2.43)-(2.45), где частные решения y (x) и их производные представлены формулой (2.35). При заданном m алгоритм 2.1 имеет сложность O(K), где K - число отрезков разбиения (2.7). Отметим, что если разбиения (2.7) является равномерным, то матрицы (2.83) совпадают между собой при всех k = 1,... , K, поскольку одинаковы длины (2.82) отрезков [ξk-1, ξk ]. В этом случае все шаги алгоритма 2.1, кроме шага 5, сводятся к перемножению матриц, что позволяет достичь высокой эффективности при его численной реализации. 3. Численные результаты Рассмотрим пример численного решения начально-краевой задачи для уравнения Гарднера xxx ∂tu + (αu + β u2) ∂xu + γ ∂3 u = 0, (3.1) т. е. для уравнения (1.1) при f (x)= 0. Выпишем известное [7] солитонное решение, т. е. решение типа уединенной бегущей волны 2 uex(x, t)= 6 γ Γ , Γ > 0, Ω := γ Γ2, (3.2) α + α2 +6 βγ Γ2 ch(Γ (x - Ω t)) где Ω - скорость распространения волны вдоль оси x, а Γ - произвольный параметр. При β < 0 допустимые значения величины Γ имеют верхнюю грань Γlim: α | | Γ ∈ (0, Γlim), Γlim = 6 β , (3.3) γ так что при Γ → Γlim амплитуда A := max (x,t)uex(x, t) решения (3.2) приближается снизу к своему предельному значению Alim = α , |β| а верхняя часть профиля волны обретает характерную плоскую форму. Зададим начальное и краевые условия (1.2), (1.3) для уравнения (3.1) в следующем виде: ∂u ex u0(x)= u ex(0, x), ϕ0(t)= u ex(0, t), ψ0(t)= u ex(L, t), ψ1(t)= (L, t). (3.4) ∂x Рис. 1. Решение задачи (3.1), (3.4), (3.5) - распространение солитона (3.2); а, б, в - график решения в различные моменты времени; г, д, е - график абсолютной погрешности в те же моменты времени. Fig. 1. Solution of the problem (3.1), (3.4), (3.5): propagation of the soliton (3.2); а, б, в: graphs of the solution at different moments of time; г, д, е: graphs of the absolute error at the same moments of time. τ m K lun - uex(·, tn)l L ∞ l∂xun - ∂xuex(·, tn )l L∞ CPU t, sec 1 e-3 4 100 1.2 e-2 2.7 e-1 0.24 1 e-4 5 200 1.3 e-4 3.2 e-3 3.48 1 e-5 6 300 1.4 e-6 3.2 e-5 50.36 1 e-6 7 400 2.9 e-8 5.1 e-7 658.11 Таб. 1. Сопоставление результатов численного алгоритма для задачи (3.1), (3.4), (3.5) с ее точным решением (3.2) при различных значениях параметров τ, m, K. Tab. 1. Comparison of the results of the numerical algorithm for the problem (3.1), (3.4), (3.5) with its exact solution (3.2) for different values of the parameters τ, m, K. Приведены результаты численного решения начально-краевой задачи (3.1), (3.4), проведенного по схеме из раздела 1.2 с использованием алгоритма 2.1 для решения вспомогательной линейной задачи (1.6), (1.8), (1.9), при следующих значениях параметров: α = 1, β = -0.1, γ = 10-3, L = 3, T = 1.6, x0 = 0.2, Alim = 10, Γ= (1 - 10-5) Γlim ≈ 40.82, Ω ≈ 1.67. (3.5) Рис. 2. Найденная численно производная решения (3.2) задачи (3.1), (3.4), (3.5); а, б, в - график производной в различные моменты времени; г, д, е - график абсолютной погрешности этого вычисления в те же моменты времени. Fig. 2. The numerically found derivative of the solution (3.2) of the problem (3.1), (3.4), (3.5); а, б, в: graphs of the derivative at different moments of time; г, д, е: graphs of the absolute error of this calculation at the same moments of time. В таб. 1 для различных сочетаний шага по времени τ, степени m кусочно-полиномиальной аппроксимации и количества K отрезков разбиения (2.7) приведена погрешность численного реex шения задачи (3.1), (3.4), (3.5), вычисленная как максимум L∞-норм lun -u (·, tn)l L∞ на отрезке [0, L] по всем n = 0,... ,N разности аппроксимации un(x) с точным решением (3.2) при t = tn. ∞ Также в таблице дана погрешность l∂xun - ∂xuex(·, tn)l L вычисления первой производной решения и затраченное процессорное время. Численная реализация проводилась при помощи библиотеки Numeric Python на ноутбуке с процессором 11th Gen Intel Core i7. Отметим, что шаг τ по времени может быть выбран очень малым без потери устойчивости численного результата. Это связано с тем, что предложенный в разделе 2 алгоритм 2.1 сохраняет свою эффективность при решении сингулярно возмущенных краевых задач для линейного уравнения (2.1), т. е. при больших значениях параметра ω, что и реализуется для вспомогательного уравнения (1.6) при малом τ. Отметим также, что в тестовых примерах, приведенных в таб. 1, достигаемая точность определяется явно-неявной схемой аппроксимации (1.5), т. е. увеличение указанных в таблице параметров m, K при заданном значении τ не приводит к дальнейшему повышению точности результата. На рис. 1 приведены графики численного решения задачи (3.1), (3.4), (3.5) при τ = 10-6,m = 7, K = 400 в различные моменты времени, а также графики погрешности, найденной как разность приближенного и точного решения (3.2). На рис. 2 даны графики вычисленной производной решения и соответствующие графики абсолютной погрешности.

About the authors

S. I. Bezrodnykh

Federal Research Center “Computer Science and Control” of the RAS

Author for correspondence.
Email: sbezrodnykh@mail.ru
Moscow, Russia

S. V. Pikulin

Federal Research Center “Computer Science and Control” of the RAS

Email: spikulin@gmail.com
Moscow, Russia

References

Supplementary files