Analysis of the difference scheme of wave equation equivalent with fractional differentiation operator

Full Text

Введение. Развитие науки и техники привело к широкому использованию на практике самых различных волновых процессов. В результате этого живые системы подвергаются воздействию волн не только естественного, но и искусственного происхождения. Можно сказать, что живые организмы везде сталкиваются с различного типа волнами. Физические поля самых разных типов в процессе эволюции живых организмов оказывают на них существенное влияние. Сравнительно давно появилось отчетливое понимание того, что живые организмы не только активно используют волновые процессы в своей деятельности, но и способны чутко реагировать на внешние воздействия, имеющие волновую природу. ISSN: 2310-7081 (online), 1991-8615 (print); doi: http://dx.doi.org/10.14498/vsgtu1209 © 2014 Самарский государственный технический университет. Образец цитирования: В. Д. Б е й б а л а е в, А. З. Я к у б о в, “Анализ разностной схемы аналога волнового уравнения с оператором дробного дифференцирования” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2014. № 1 (34). С. 125-133. 125 В. Д. Б е й б а л а е в, А. З. Я к у б о в Ответная реакция живой системы на внешние волновые воздействия может происходить на различных структурных уровнях живого организма - от молекулярного до субклеточного, клеточного уровня. Волновые процессы активно используются и в технической деятельности людей. Несмотря на значительные усилия исследователей, до сих пор задача создания адекватных количественных моделей волновых процессов остаётся актуальной. Особенно это актуально, когда речь идет о системах с фрактальной структурой. При описании свойств систем с фрактальной структурой нельзя использовать представления евклидовой геометрии, поэтому необходимо привлечь представления геометрии дробной размерности. Особенность систем с фрактальной структурой в том, что для них существенны такие эффекты, как память, сложная природа пространственных корреляций и эффекты самоорганизации. Создание адекватных математических моделей для систем, где проявляются свойства самоорганизации, детерминированного хаоса, также требует привлечения нетрадиционных подходов, основанных на применении математического аппарата дифференциальных уравнений дробного порядка [1-7]. В качестве математической модели волнового процесса в нелокальных средах исследуем в области D = {(x, t) : 0 x l, 0 t T} первую краевую задачу для волнового уравнения с производными дробного порядка. Задача. Найти решение уравнения β α D0t u(x, t) = C(x, t)CD0x u(x, t) + f (x, t), C (1) удовлетворяющее начальным условиям u(x, 0) = ϕ1 (x), ut (x, 0) = ϕ2 (x) (2) u(l, t) = µ2 (t). (3) и граничным условиям u(0, t) = µ1 (t), Здесь 1 < α 2, 1 < β C 2, f (x, t) - достаточно гладкая функция, C(x, t) α D0t u(x, t) = 1 Γ(2 - α) t 0 0, u (x, s) ds (t - s)α-1 - частная дробная производная Капуто [10]. Численный метод решения краевой задачи (1)-(3). Для численного решения задачи (1)-(3) в области D введём сетку с шагом h по x и τ по t: hτ = (xm , tn ) : xm = mh, tn = nτ ; h = l/M, τ = T /N ; m = 0, 1, . . . , M, n = 0, 1, . . . , n . Для дробных производных Капуто в случае 1 < α 2, 1 < β 2 на частичных отрезках [tn-1 , tn ], [xm-1 , xm ] имеют место разностные аппроксимации [11] 126 Анализ разностной схемы аналога волнового уравнения . . . 1 Γ(2 - α) tn tn-1 u (x, s) ds = (tn - s)α-1 = 1 Γ(2 - β) xm xm-1 u(x, tn+1 ) - 2u(x, tn ) + u(x, tn-1 ) + O(τ 4-α ), (4) Γ(3 - α)τ α u (s, t) ds = (xm - s)β-1 = u(xm+1 , t) - 2u(xm , t) + u(xm-1 , t) + O(h4-β ). (5) Γ(3 - β)hβ Обозначим через un ≈ u(xm , tn ) приближённое решение в точке (xm , tn ), m n n n = 0, 1, . . . , N , m = 0, 1, . . . , M , аналогично обозначим Cm = C(xm , tn ), fm = = f (xm , tn ). Воспользовавшись (4) и (5), для уравнения (1) получим разностную схему с весами un+1 - 2un + un-1 m m m n = σΛβ un+1 + (1 - σ)Λβ un + fm , m m Γ(3 - α)τ α (6) где Λβ un+1 m = n+1 n um+1 Cm - 2un+1 + un+1 m m-1 , β Γ(3 - β)h n Λβ un = Cm m un - 2un + un m m+1 m-1 . β Γ(3 - β)h В случае σ = 0 получаем явную разностную схему на 5-точечном шаблоне: n - 2un + un un+1 - 2un + un-1 m m-1 m m m n u n = Cm m+1 + fm . α β Γ(3 - α)τ Γ(3 - β)h (7) Теорема 1. Явная разностная схема (7) устойчива, если CΓ(3 - α)τ α < 1, Γ(3 - β)hβ где C = max |C(x, t)|. D Д о к а з а т е л ь с т в о. Решение un+1 можно представить в виде un+1 = Aun + Γ(3 - α)τ α f n . Тогда S = A - оператор перехода с одного временного слоя на другой. Условием устойчивости по начальным данным разностной схемы (6) является S 1. Действительно, y n+1 = Ay n или y n+1 = y n A · yn ... y0 . Следовательно, начальные возмущения затухают. Условие S 1 есть условие того, что спектр оператора S лежит внутри круга единичного радиуса на комплексной плоскости. Это означает, что max |λ| < 1, где λ - собственное число оператора перехода S. 127 В. Д. Б е й б а л а е в, А. З. Я к у б о в Для нахождения собственных значений оператора перехода решение un m представим в виде возмущения: un = λn eiωm , m (8) где i - мнимая единица, λ - собственные числа оператора перехода [12]. В результате подстановки (8) в (7) получим квадратное уравнение λ2 - 2 1 - n ω 2Cm Γ(3 - α)τ α sin2 λ + 1 = 0. β 2 Γ(3 - β)h (9) Произведение корней уравнения (9) равно единице, а его дискриминант D= n n ω Cm Γ(3 - α)τ α ω 16Cm Γ(3 - α)τ α sin2 sin2 - 1 . β β 2 2 Γ(3 - β)h Γ(3 - β)h (10) Если дискриминант (10) отрицательный, то корни λ1 , λ2 - комплексно-сопряжённые и равны единице по модулю. Следовательно, условие |λ| < 1 выполняется, если они выбраны так, что CΓ(3 - α)τ α < 1. Γ(3 - β)hβ Итак, решение un разностной задачи m u0 m n - 2un + un un+1 - 2un + un-1 m m-1 m m m n u n = Cm m+1 + fm , Γ(3 - α)τ α Γ(3 - β)hβ = ϕ1 (xm ), u1 = u0 + τ ϕ2 (xm ), un = µ1 (tn ), un = µ2 (tn ), m m 0 M можно найти явным методом. В случае σ = 1 получим полностью неявную схему n+1 - 2un+1 + un+1 un+1 - 2un + un-1 m m-1 m m m n u n = Cm m+1 + fm . Γ(3 - α)τ α Γ(3 - β)hβ (11) Теорема 2. Неявная разностная схема (11) безусловно устойчива. Д о к а з а т е л ь с т в о. Для нахождения собственных значений оператора перехода решение un представим в виде возмущения (8). В результате подm становки (8) в (11) получим квадратное уравнение [1 + p] · λ2 - 2λ + 1 = 0, где p= n ω 4Cm Γ(3 - α)τ α sin2 . β 2 Γ(3 - β)h Дискриминант этого квадратного уравнения отрицательный. Таким образом, корни комплексно-сопряжённые и имеет место неравенство √ p |λ| = < 1. 1+p 128 Анализ разностной схемы аналога волнового уравнения . . . Следовательно, неявная разностная схема (11) безусловно устойчива. Для n = 1, 2, . . . , N имеем tn 1 Γ(2 - α) α (C D0t u(x, t))n = 0 u (x, s) ds = (tn - s)α-1 = = = = 1 Γ(2 - α) 1 Γ(2 - α) 1 Γ(3 - α) n i=1 n i=1 n i=1 1 Γ(2 - α) n ti i=1 ti-1 ti u(x, ti+1 ) - 2u(x, ti ) + u(x, ti-1 ) + O(τ 2 ) τ2 u(x, ti+1 ) - 2u(x, ti ) + u(x, ti-1 ) +O(τ 2 ) τ2 u (x, s) ds = (tn - s)α-1 ti-1 - ds = (tn - s)α-1 (tn - s)2-α 2-α ti ti-1 = u(x, ti+1 ) - 2u(x, ti ) + u(x, ti-1 ) + O(τ 2 ) × τ2 × (tn - ti-1 )2-α - (tn - ti )2-α = = 1 Γ(3 - α)τ 2 n 2-α u(x, ti+1 ) - 2u(x, ti ) + u(x, ti-1 ) (t2-α - tn-i )+ n-i+1 i=1 + O(τ 2 ) 1 = Γ(3 - α)τ 2 2-α (tn-i+1 - t2-α ) = n-i n u(x, tk+1 ) - 2u(x, tk ) + u(x, tk-1 ) × k=1 × t2-α - t2-α + O(τ 3-α ), (12) n-k n-k+1 так как t2-α - t2-α τ 2-α . n-i+1 n-i Аналогично для m = 1, 2, . . . , M имеем β (C D0x u(x, t))m = 1 Γ(2 - β) xm 0 1 = Γ(2 - β) = 1 Γ(3 - β)h2 u (s, t) ds = (xm - s)β-1 m xk k=1 xk-1 u (s, t) ds = (xm - s)β-1 m u(xk+1 , t) - 2u(xk , t) + u(xk-1 , t) × k=1 2-β 3-β ), (13) × x2-β m-k+1 - xm-k + O(h Воспользовавшись разностными аппроксимациями (12) и (13), получим составную разностную схему с весами: 1 Γ(3 - α)τ 2 n uk+1 - 2uk + uk-1 t2-α - t2-α = m m m n-k+1 n-k k=1 129 В. Д. Б е й б а л а е в, А. З. Я к у б о в n = σ · Λβ un+1 + (1 - σ)Λβ un + fm , m m где n Λβ un+1 = Cm m Λβ un m = n Cm 1 Γ(3 - β)h2 1 Γ(3 - β)h2 m 2-β un+1 - 2un+1 + un+1 xm-k+1 - x2-β , k+1 k k-1 m-k k=1 m 2-β 2-β un - 2un + un k+1 k k-1 xm-k+1 - xm-k . k=1 Пример. Найти решение задачи β α D0t u(x, t) = C D0x u(x, t), 0 < x < 4, 0 < t < 10, u(x, 0) = sin(πx/4), ut (x, 0) = 0, 0 x 4, u(0, t) = 0, u(4, t) = 0, 0 t 10. C (14) На рис. 1-3 приведены графики решений задачи (14) при различных значениях параметров α, β. Как видно из рис. 1, при переходе к дробной производной по времени период колебаний становится меньше, а при переходе к дробной производной по координате, наоборот, период колебаний растёт. При этом длина волны остаётся неизменной (рис. 2 и 3). Следовательно, при переходе к дробным производным в первом случае скорость распространения волны становится больше, чем в стационарном случае, а во втором случае - меньше. Выводы. Значительные успехи при рассмотрении волновых процессов в системах с фрактальной структурой связаны с использованием формализма интегралов и производных дробного порядка. Повышенный интерес к дифференциальным уравнениям дробного порядка обусловлен их физической интерпретацией [7, 13]. Было показано, что переход к производной дробного порядка по времени позволяет учитывать эффекты памяти системы [7]. Это и вызвало их широкое применение в практически во всех областях естествознания. Область применимости решений дифференциальных уравнений дробного порядка значительно шире, чем дифференциальных уравнений с целочисленным дифференцированием, поскольку последние оказываются их частным случаем. Уравнения в производных дробного порядка позволяют учесть обратимые и необратимые процессы [7, 13]. Это позволяет не только получить принципиально новые, но более глубоко осмыслить известные результаты, позволяя при этом создать адекватные количественные модели исследуемых явлений Исследование свойств систем с фрактальной структурой должно быть проведено с учётом эффектов памяти (временная нелокальность) и пространственных корреляций (пространственная нелокальность). Учёт эффектов нелокальности в рамках традиционных подходов приводит к появлению в дифференциальных уравнениях интегрального оператора, где ядро интегрального оператора несёт информацию о природе нелокальности. Для решения таких уравнений интегральные операторы представляются в виде ряда дифференциальных операторов с возрастающим показателем порядка дифференцирования и при наличии малого параметра ограничиваются несколькими членами ряда. В отсутствие малого параметра такой подход оказывается 130 Анализ разностной схемы аналога волнового уравнения . . . Рис. 1. Графики решения задачи (14) при x = 0.7 и различных значениях параметров α и β: 1 - α = 2, β = 2; 2 - α = 2, β = 1.7; 3 - α = 1.7, β = 2 Figure 1. The graphics of solution (14) for x = 0.7 and different values of parameters α and β: 1 - α = 2, β = 2; 2 - α = 2, β = 1.7; 3 - α = 1.7, β = 2 Рис. 2. Графики решения задачи (14) при t = 0.4 и различных значениях параметров α и β: 1 - α = 2, β = 2; 2 - α = 1.7, β = 2; 3 - α = 1.5, β = 2; 4 - α = 1.4, β = 2 Figure 2. The graphics of solution (14) for t = 0.4 and different values of parameters α and β: 1 - α = 2, β = 2; 2 - α = 1.7, β = 2; 3 - α = 1.5, β = 2; 4 - α = 1.4, β=2 Рис. 3. Графики решения задачи (14) при t = 0.9 и различных значениях параметров α и β: 1 - α = 2, β = 2; 2 - α = 2, β = 1.7; 3 - α = 2, β = 1.5; 4 - α = 2, β = 1.4 Figure 3. The graphics of solution (14) for t = 0.9 and different values of parameters α and β: 1 - α = 2, β = 2; 2 - α = 2, β = 1.7; 3 - α = 2, β = 1.5; 4 - α = 2, β = 1.4 131 В. Д. Б е й б а л а е в, А. З. Я к у б о в непродуктивным и, кроме того, полученные уравнения также не всегда удаётся разрешить.

About the authors

Vetlugin D Beybalaev

Dagestan State University

Email: kaspij_03@mail.ru
(Cand. Phys. & Math. Sci.), Associate Professor, Dept. of Applied Mathematics 43a, M. Gadzhiev st., Makhachkala, 367025, Russian Federation

Amuchi Z Yakubov

Daghestan State Institute of National Economy

Email: yakubovaz@mail.ru
(Cand. Phys. & Math. Sci.), Associate Professor, Dept. of Computer Science & Computer Engineering 5. D. Ataev st., Makhachkala, 367008, Russian Federation

References

Supplementary files