CC BY
2
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
2023, Т. 165, кн. 3 С. 190-207
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
ОРИГИНАЛЬНАЯ СТАТЬЯ
УДК 519.63 10.26907/2541-7746.2023.3.190-207
КОНСЕРВАТИВНАЯ ПОЛНОСТЬЮ ДИСКРЕТНАЯ СХЕМА МКЭ ДЛЯ НЕЛИНЕЙНОГО УРАВНЕНИЯ КЛЕЙНА-ГОРДОНА
Р. З. Даутов, Г. Р. Салимзянова
Казанский (Приволжский) федеральный университет, г. Казань, 420008, Россия
Аннотация
Предложено семейство методов Петрова-Галеркина-МКЭ для решения нелинейного уравнения Клейна-Гордона. Дискретные схемы сформулированы в терминах решения задачи и его производной по времени и обеспечивают сохранение полной энергии на дискретном уровне. Численно исследована простейшая двухслойная схема подобного типа. На основе решения тестовых задач с гладкими решениями показано, что схема позволяет определить как решение задачи, так и его производную по времени с погрешностью порядка 0(Я2 + т2) в среднеквадратической норме, где т и Я характеризуют шаги сетки по времени и пространству соответственно.
Ключевые слова: метод Петрова-Галеркина, метод конечных элементов, уравнение Клейна-Гордона, двухслойная схема
Введение
1. Нелинейное уравнение Клейна-Гордона (или нелинейное волновое уравнение)
u''(x,t) — △ u(x,t) + y>(u(x,t)) = 0 (1)
является одним из классических уравнений теории нелинейных волн. Здесь x = (xi, Х2,..., xd) - пространственная переменная, u' = du/dt - производная u(x, t) по времени t, ^ - заданная функция, Д - оператор Лапласа,
д2и д2и
= тт^т + • • • + тт^т • dxi dxd
Многочисленные приложения уравнения (1) обусловлены достаточно типичным сочетанием модельной (квадратичной) дисперсии и нелинейности. Среди наиболее «востребованных» в приложениях нелинейностей выделим нелинейность
¡р(и) = sin u (2)
(уравнение синус-Гордона) и кубическую нелинейность
tf(u) = —u + u3. (3)
Обзор задач, приводящих к уравнению Клейна-Гордона, можно найти в монографиях [1,2]. Среди приложений этого уравнения отметим теории магнетиков, дислокаций, джозефсоновских переходов и целый ряд других областей физики.
Наряду с уравнением (1) большой интерес представляют его модификации вида
где а, а и в являются положительными постоянными или функциями. Слагаемые в правой части (4), а также функция а описывают наличие в моделируемой системе диссипативных эффектов, неоднородностей и других возмущений.
Если уравнение (1) (например, в случаях (2) и (3)) допускает точные аналитические решения, то модифицированное уравнение (4) слишком сложно для точного интегрирования, и его решения могут быть получены только приближенно с использованием численных методов.
2. Уравнение (1) является уравнением гиперболического типа, и для его приближенного решения могут быть использованы соответствующие известные численные методы. В связи с этим отметим следующие работы, посвященные сеточным методам (метод конечных разностей [3-7], метод конечных элементов [8,9]), а также работы, посвященные бессеточным методам (спектральные и псеводоспектральные методы [10], дифференциальный метод квадратур [11], метод радиальных базисных функций [12] и т. д.).
Поскольку уравнение Клейна-Гордона описывает консервативные системы, то следует отметить, что консервативные численные схемы его решения демонстрируют лучшие результаты, чем неконсервативные, и при прочих равных условиях являются более точными. Ключевым моментом является то, что такие схемы сохраняют некоторые инвариантные свойства дифференциального уравнения и более точно воспроизводят особенности решения при вычислениях на больших отрезках времени. Для решения уравнения Клейна-Гордона нам известны только консервативные конечно-разностные схемы на ортогональных сетках (см., например, [6,7]). Консервативные сеточные схемы на основе метода конечных элементов, предназначенные для решения задачи в произвольных областях, нам неизвестны.
3. Основной целью настоящей работы является численное исследование точности предлагаемого нами консервативного сеточного метода решения уравнения (1). Он также применим для решения уравнения (4). Первоначально введением новой переменной V = и' уравнение сводится к системе уравнений, которая далее и аппроксимируется на основе комбинации метода конечных элементов (по пространственным переменным) и метода Петрова-Галеркина (по временной переменной). Метод является т-слойным и имеет формальный порядок аппроксимации 0(тт + Нк), где т и Н характеризуют шаги сетки по времени и пространству соответственно, т > 2, к > 1. Результаты проведенных вычислительных экспериментов для двухслойной схемы на основе линейных треугольных конечных элементов позволяют утверждать, что метод эффективен и пригоден для практического использования.
1. Постановка задачи
1.1. Исходная задача. Пусть Т > 0, О - ограниченная область в Кй с лип-шицевой границей Г, d > 1, (т = О х (0, Т]. Рассмотрим неоднородное уравнение Клейна-Гордона
и" — &у(а gradu) + <р(и) = -аи' + в △и' + 7,
(4)
и''(х,г) — Аи(х,г) + р(и(х, г)) = /(х, г), (х,г) е ,
при краевых условиях Дирихле
1
и(х,г) = 7(х,г), (х,г) е Г х (0,Т].
хДля удобства изложения будем считать, что функция 7 определена в Qт •
В начальный момент времени t = 0 решение уравнения и его производная по времени считаются известными:
u(x, 0) = uo(x), u'(x, 0) = vo(x), x G (6)
Сформулированная выше задача является классической и хорошо исследованной, условиям существования ее обобщенных решений посвящено множество работ. В связи с этим отметим лишь статью [13] и книгу [14, Гл. 1]. Для определения дискретного метода нам понадобится ее обобщенная формулировка, изученная в [14, с. 20] в случае степенной функции ^. С этой целью введем дополнительные обозначения и определим пространства функций.
1.2. Пространства функций. Пусть а • b = ^j «Л обозначает скалярное произведение векторов а и b, |а| = (а • а)1/2 - длина вектора а. Используем стандартные обозначения: Vu - для градиента функции u,
/ du du
vdxi ' ' dxd
Lp(0) и Hm(0) - для обозначения пространств функций Лебега и Соболева, p G [1, те], m > 1, H1(0) = {u G H 1(0) : u|r = 0}.
Положим H = ¿2(0) и для конечного p > 1 определим пространства Vp = H 1(0) П Lp(0) и Vp0 = Hq(0) П Lp(0). Через (•, •) будем обозначать как скалярное произведение в L2(0),
(u,v) = / JQ
(u, v ) = I u v dx, '0
так и отношение двойственности между Ур0 и его сопряженным Н + Ьр (О), где р' определяется как сопряженный к р показатель, 1/р + 1/р' = 1. Пусть также
(Ум, V«) = / Ум • У«вх. Jn
Отметим, что вложение
Н 1(О) с Ьр(П) (7)
непрерывно при любом конечном р при в < 2 и р < 2в/(в — 2) при в > 3. Следовательно, при таких р имеем Ур = Н 1(О), Ур = Нд(О).
Далее для заданной функции м переменных х С О и 4 С [0, Т] под м(4), 4 € [0, Т], будем понимать функцию х ^ м(х, 4), х С О. Подобные функции, определенные на [0, Т] со значениями в банаховом пространстве В функций от х, будем рассматривать как элементы пространства Ьр(0, Т; В), т. е. пространства (классов) измеримых функций [0, Т] ^ В таких, что
( fT \ 1/Р
u||lp(0,t;B) =У J ||u(î)||B dtj < те, 1 < p< те.
При р = те норма в нем определяется функционалом эир е884е(о,т)1|м(£)||.в (подробнее, см., например, [13, с. 469]). Отметим, что Ьр(0,Т; Ьр(О)) = ¿р(ф).
1.3. Обобщенное решение задачи. Уравнение (5) перепишем в виде системы, введя дополнительное неизвестное « = м':
«'(*) — Дм(4) + ¥>(м(*)) = / (4), (8)
м'(4) = «(4), (9)
где г € (0,Т]. Краевые и начальные условия примут вид
и(£)|г = 7 (г), и(0) = ио, и(0) = ио.
Считая, что функция р является непрерывной на всей числовой оси, определим ее потенциал Ф:
Ф(а) = Фо + I АО ¿£, а € К,
о
о
где постоянная Фо может быть выбрана произвольно. Дополнительно будем предполагать, что для некоторого конечного р > 1 функция р удовлетворяет следующим условиям:
(Н\) если и € Ур, то /о Ф(и(х)) ¿х < ж, р(и) € Ьр'(£1).
Относительно других данных задачи будем считать, что (Н2) I € Ь2(), и0 € Ур, у0 € Н, ^ является следом на Г х [0,Т] некоторой функции п из Ьж(0,Т; Ур) такой, что г' € Ьж(0,Т; У*).
Отметим, что функции, указанные в (2), (3), удовлетворяют условиям (Н\) при р = 2 (например) и р = 4 соответственно. В этом нетрудно убедиться, если принять во внимание неравенство Гельдера для интегралов и вложение (7).
Для определения обобщенного решения задачи умножим обе части уравнения (8) на произвольную функцию С € Ь2(0,Т; Ур0), а уравнение (9) - на произвольную функцию г € Ь2(0,Т; Н). Полученные равенства проинтегрируем по области (т и учтем, что
- Аи(г) ((г) ¿х = Уи(г) ■ У((г) ¿х, ■ 'о -'о
/о ио
поскольку ф)1г = 0 .В результате придем к тождествам
Г (у'(г), с (г)) ¿г + Г а(и(г), С (г)) ¿г = Г (I (г), с (г)) ¿г, (10)
о о о
[ (и'(г),п(г)) ¿г = ( (у(г),/(г)) ¿г, (11)
оо где функционал а определен следующим образом:
а(и,( ) = [ (Уи ■УС + р(и) С) ¿х. (12)
о
Определение 1. Пару функций и,у, удовлетворяющую условиям
и - 7 € Ь2(0,Т; Ур0), и',ь € Ь2(0,Т; Н), V' € Ь2(0,Т; У*), и(0) = ио, v(0) = vо,
а также тождествам (10), (11) для всех ( € Ь2(0,Т; У£), / € Ь2(0,Т; Н), назовем обобщенным решением задачи (5)-(6).
Заметим, что при выполнении условий (Н\), (Н2), данное определение корректно в том смысле, что все слагаемые в тождествах (10), (11) являются конечными, а начальные условия имеют смысл (см., например, [14, Замечание 1.1, с. 21]).
Доказательство существования определенного выше обобщенного решения в случае, когда
р(и) = 1и1ри, р> 0, 7 = 0, р = р + 2, (13)
дано в [14, Теорема 1.1, с. 20]. Оно непосредственно обобщается для рассматриваемой нами задачи. Единственность решения имеет место при дополнительных к (Н1) условиях на функцию ^ (например, условии липшиц-непрерывности оператора из Ур в V*, порожденного этой функцией). Эти условия выполнены для функций (2), (3), а в случае (13) они приводят к ограничению р < 2/(п — 2) при в > 3 (см., например, [14, Теорема 1.2, с. 27]).
Отметим следующее важное и хорошо известное свойство решения рассматриваемой задачи, справедливое при / = /(х), 7 = 7(х). Определим функционалы
П(и)= [ (-№и\2 + Ф(и))(1,х, Е(1>,и) = -(г>,1>) + Щи)-(/,и), .¡о ^2 ) 2
которые будем называть функционалами потенциальной и полной энергии соответственно. Они определены на функциях м С Ур и « С Н. Справедливо утверждение (см., например, [14, Теорема 1.6, с. 36]):
#(«(г),м(г)) = Е(«0,м0) для всех 4 С [0, Т], которое принято называть законом сохранения (полной) энергии.
2. Определение схемы Петрова—Галеркина—МКЭ
Укажем способ дискретизации поставленной выше задачи, обеспечивающий сохранение энергии на дискретном уровне. Его можно рассматривать как вариант метода Петрова-Галеркина с возмущениями, основанный на аппроксимациях, принятых в методе конечных элементов. Полностью дискретную схему получим последовательной дискретизацией задачи сначала по пространственным переменным, а затем - по переменной времени.
Для облегчения изложения далее будем предполагать, что О является многогранником в .
2.1. Полудискретная схема МКЭ с численным интегрированием.
Пусть Н > 0 означает максимальный размер конечных элементов, ТН = {К., К2,..., К^} - точное разбиение (триангуляция) области О семейством лагран-жевых конечных элементов К1,..., К какой-либо формы. На основе этого разбиения определим пространство конечных элементов
14 = к е С(П): 1>ь\к е Рк УК е Тй},
где Рк - некоторое конечномерное пространство функций на элементе К. Определим также его подпространство
У0 = {«н С Ун : «н|г = 0}.
Пространство Ун (У°) будем рассматривать как аппроксимацию Н 1(О) (Нд(О)). Ясно, что справедливы вложения Ун С Ур, У° С Ур0.
Определим в Ун базис Лагранжа {^¿(х)}|=1, соответствующий множеству узлов интерполяции шн = , N = N(Н), т.е. функции С У, такие, что ^г(аз) = , .7 = 1,..., N, где - символ Кронекера. Тогда для произвольной функции С Ун справедливо разложение
N
«н(х) = «г ¥>г(х),
г=1
где VI = у-(а^), г = 1, . . . , N - узловые параметры функции V-.
Функцию гч из V/) назовем V/, -интерполянтом функции и (Е С(Г2), если и/(а^) = и(щ), г = 1,..., N, т. е. функцию
N
и/(х) = ^и(а*) (х).
г=1
Для удобства изложения далее будем предполагать, что внутренние узлы сетки имеют номера от 1 до п, а граничные (лежащие на Г) - от п +1 до N.
Для вычисления интеграла от функции I по конечному элементу К зададим некоторую квадратуру БК (I) с Ь узлами:
ь
Кг/
f f (x) dx « £ c?f (d?) = Sk(f ),
Jk i=i
где c?, dK - веса и узлы квадратуры. Квадратурную формулу для вычисления интеграла по области ^ определим как составную:
i f (x) dx = £ i f (x) dx « £ S? (f ) = Sn(f )■ Jq KeTh K ?eTh
Через (•, •)h обозначим аппроксимацию скалярного произведения в Ь2(0.),
(u,v)h = Sq(uv).
Будем предполагать, что оно определяет скалярное произведение в Vh. Введем также следующие обозначения
(Vu, Vv)h = Sn(Vu^Vv), ah(u,Ç) = (Vu, VÇ)h + (f(u),Ç)h. (14)
Далее через u0h и v0h будем обозначать функции из Vh, которые являются некоторыми аппроксимациями функций uo и vo соответственно.
Определение 2. Функции uh,vh из HX(0,T; Vh), удовлетворяющие условиям uh - yi G H\0,T; V0), uh(0)= uoh, vh(0) = voh, а также тождествам
f (vh(t),C(t))h dt + f ah(uh(t)X(t)) dt = f (f (t),C(t))h dt,
o o o
f (uh(t), n(t))h dt = f (vh(t),n(t))h dt
oo
для всех Z G L2(0,T; V°), n G L2(0,T; Vh), назовем приближенным решением задачи (5)-(6) по методу конечных элементов с численным интегрированием.
2.2. Полностью дискретная задача. Пусть т > 0. Определим на отрезке [0, Т] сетку узлов шт = {Ь^, ] = 0,1,..., М} с шагами т^ = Ь^ — Ь^-1, т^ < т, и фиксируем целое т > 1. Обозначим через и-т(Ь) непрерывную на [0, Т] функцию со значениями в V-, которая является алгебраическим многочленом степени т на каждом отрезке ]. Множество всех таких функций обозначим
через (0,Т; Ун). Определим также пространство 5т-1(0, Т; Ун), состоящее из функций £ : [0, Т] ^ Ун, являющихся многочленами степени т — 1 на каждом отрезке [¿¿-1,4]]. Функции из 5т-1(0, Т; Ун) не предполагаются непрерывными на всем [0, Т] (они могут терпеть разрыв в узлах сетки). Отметим, что кусочно-определенная производная мНт принадлежит 5т-1(0,Т; Ун).
Полностью дискретную схему Петрова-Галеркина-МКЭ для задачи (5)-(6) определим следующим образом.
Определение 3. Найти такие функции мНт,«Нт из (0, Т; Ун), что мнт — 7/т С (0, Т; Ун0), мнт(0) = мон, «нт(0) = «он, которые удовлетворяют тождествам
[ Кт(4), С(4))н +/ ан(мнт(4), С(^)) в4 =/ (/(4),Ш)н в4, (15)
0 0 0
I (мнт(*),п(*))н в4 = I («нт(*),п(*))н в4 (16)
для всех С С 5т-1(0, Т; Ун0), п С 5^(0, Т; Ун).
(0, Т; н
Отметим, что условие мнт — 7/т С (0, Т; У°) равносильно условию
мнт(аг,4]) = 7(аг, ) для всех г = п + 1,..., N, = 1,..., М.
2.3. Сохранение энергии на дискретном уровне. Определим на пространстве Ун дискретные аналоги функционалов П и Е:
ПнЫ = + Ф(«)), Ен(г,,и) = ^(г»,г»)л. +Щ(г») - (/,«)/,.
Теорема 1. Пусть / = /(х), 7 = 7(х), мнт, «нт - решения схемы (15), (16). Тогда Ен(«нт(4), мнт(4)) = Ен(«0н,м0н) для всех 4 С Шт .
Доказательство. Поскольку 7 не зависит от 4, то 7/т совпадает с 7/ и также не зависит от 4. Поэтому мнт (4) С У° для всех 4. Это позволяет выбрать в (15) и (16) пробные функции равными
л(/) = / Чт4 < ], п(/) Г «нт 4 < ,
Цг) = \ 0, 4 >4], п(4) = \0, 4 >4],
где ^ > 1. Придем к следующим равенствам
Г 1з Г 1з Г 1з
/ Кт (4)Хт(4))н в4 + / ан(мнт (4)Хт(4)) в4 = (/,мнт(4))н в4,
0 0 0 Г1з Г1з
/ (мнтИХт(4))н в4 = («нт(4),«нт(4)) в4,
00
из которых вытекает соотношение
Г 1з Г 1з Гьз
/ («нт(4),«нт(4))н в4 + / ан(мнт(4), мнт(4)) в4 = (/,мнт(4))н в4. (17)
0 0 0
Преобразуем это равенство. Заметим, что (4) = 1/2в/в4ад2(4),
0
0
и, как следствие, для любых функций u и v справедливы следующие равенства
(v(t),v'(t))h = 1/2 d/dt (v(t),v(t))h,
ah{u{t),u'{t)) = h{u{t)), {f,u'{t))h=d/dt{f,u{t))h.
Поэтому равенство (17) можно записать в следующем виде
Г Ь ^
] -¿¿Ен{1'нт$),инт^))<и:. = 0. Следовательно, = Ен(г'0ь,иоь), 3 = 1, • • •, М. □
3. Реализация простейшей двухслойной схемы
Рассмотрим реализацию определенной выше схемы в случае т = 1. Положим
3) > К
j = uhT (tj ), uJh G Vh, j = 0,1, . . . , M. Тогда справедлива формула
икт(1) = и( 1 + (и( - и( х) Ц^А * € [^-1,^].
тз
Аналогичное представление имеет функция гнт € (0,Т; Vн). Пусть в момент времени Ьз^1 решение схемы, а именно, пара функций и^ 1, г^ 1, известно. Это имеет место при ] = 1, поскольку в силу начальных условий имеем и^ = ион, г" = ион. Выберем в тождествах (15) и (16) функции ((Ь) и п(Ь) равными нулю всюду, кроме отрезка ], на котором примем их равными соответственно
произвольно фиксированным (н € Vи щ € Vн. Тогда придем к следующим уравнениям
Ггз Ггз Ггз
/ (г'нт(Ь),Сн)н ¿Ь + ан(инт(Ь)Хн) ¿Ь = (/(Ь)Хн)н ¿Ь, (18) ГЬ ГЬ
/ (и'нт(Ь),пн)н ¿Ь = / (гнт(Ь),пн)н ¿Ь. (19)
В силу линейности функции Ь ^ инт(Ь) на отрезке \tj-i,Ьз] и определения (14) формы ан имеем
1 г Ь г1
- [ ah(uhT(t)Xh)àt-= [ ah(cr u3h + (I - a)ul 1 ,(h) da = i Jtj-i Jo
Л ' h
2h
где в силу равенства ^(s) = Ф'(в)
Ф(и) — Ф^)
Ф[u,v] = <p(a u + (1 — a)v) da =
Jo
Вычислив остальные интегралы по переменной Ь в (18) и (19), придем к следующим уравнениям, справедливым для всех ] = 1,... ,М,
j j—i
(Vh~Jh (v«ft,va) + (*HA-%Ch)h = (fj,çh)h, (20)
V Tj J h V У h
j j — 1 j j — 1 " (21)
Tj J h V 2 J h
uv
Здесь использованы обозначения
7—1 1 г Ь
йн = ин+^~ , Р = -\ та. (22)
2 Т7 4-1
В силу произвольности п7 € У7 из (21) следует
7 7—1 7 , 7—1 и, — Щ, V? + VÍ , -= И т и—^ ^
т7 2
Отсюда найдем v7 и из (22) и^. Тогда уравнение (20) примет вид
т—2 (ил,Сл)л + (Vил, УСл)л + (Ф[2ил - и7л—1,и7л—1],Сл)л = (/,Сл)л, (24) где С7 - произвольная функция из У^1,
/ = Я + 7 + 1)/г2, Г = Т7/2 .
Уравнение (24) сводится к системе (нелинейных) алгебраических уравнений.
7 = 2ил - и7—
После его решения найдем ил, а потом определим и^ = 2ил — и^ 1 и v7 из (23).
3.1. Реализация схемы в случае линейных элементов. Уточним расчетные формулы, указанные выше, при использовании линейных конечных элементов, когда конечные элементы из Тл являются симплексами в Д (отрезками -при в = 1, треугольниками - при в = 2, тетраэдрами - при в = 3), а функции из У7 на каждом конечном элементе К являются аффинными функциями, т. е. Рк в определении пространства конечных элементов У7 есть множество полиномов вида ао + 01x1 + ... + а^х^. Сетка узлов для указанных элементов совпадает с множеством всех вершин симплексов из Тл .
Напомним, что функция V7 € У7 имеет представление
N
V7(x) = V» уч(х), (25)
¿=1
где v¿ = V7 (a¿), г = 1,..., Ж, - узловые параметры функции V7.
Укажем подходящие для использования квадратурные формулы в одномерном и двумерном случаях. В одномерном случае, т. е. при в = 1 , для вычисления интеграла на отдельном элементе можно использовать квадратуру либо трапеций, либо Симпсона. Эти же квадратуры, а также квадратуру центральных прямоугольников можно использовать для вычисления интеграла, определяющего / . В двумерном случае, т. е. при в = 2 , можно использовать квадратуру вида
1к Пх) с1х « 1|1 (/(<**) + /(4') + /(4')),
где |К | - площадь К € Т7, а узлы вк совпадают либо с вершинами К, либо с серединами сторон К . Первая квадратура точна на многочленах первой степени, вторая - на многочленах второй степени. Легко убедиться в том, что в этих случаях билинейный функционал (•, •) 7 определяет скалярное произведение в У7 .
Реализация схемы, описанной в предыдущем пункте, требует определения функций иол и vо7, а также решения уравнения (24). Рассмотрим эти этапы.
1) Задание начальных значений и07 и v07 . Укажем два способа. Для непрерывных функций ио и vо в качестве иол и vо7 можно выбрать их У7 -интерполянты, т. е. функции
N N
иол(х) = ио^) ¥(х), иол(х) = ^ ио^) ¥(х).
¿=1 ¿=1
В общем случае определим начальные значения как наилучшие приближения функций ио и vо в дискретном ¿2^), т.е. как решения уравнений
(иол, С)л = (ио, С)л V С € Ул, ^ол,п)л = К,п)л Vп € Ул. (26)
Это равносильно решению двух систем алгебраических уравнений. В самом деле, выбрав С = ¥ в (26) и подставив разложение функции иол в первое равенство в (26), придем к системе уравнений
Мио = М = {(^^, ^)л}&=1, ^ = {(ио, >N=1
относительно вектора узловых параметров ио = (иол(а1),..., u07(aN))т функции иол . Если узлы квадратуры на конечном элементе совпадают с его вершинами, то матрица масс будет диагональной. Аналогичный вид имеет система уравнений для определения функции v07 .
2) Представление решений. Рассмотрим теперь представления функций и^, ] = 1,..., М, а также пробных функций из Ул в (20) и (21). Поскольку функции из равны нулю на границе Г, то
п
Сл(х) = £^¿(х), VСл € уо,
¿=1
т.е. система функций {^¿}П=1 образует базис в У^.
Согласно определению 2.2. имеем и^|г = 7/(£/), т.е. и^(а; г = п + 1,..., N. Поэтому в силу (25) справедливо разложение
п N
Ч(х) = Х!и/ ¥(х)+ 7/(x), 7/(х) = ^2 7КЛ') ¥(х).
¿=1 ¿=п+1
Через и7 = ( и 1,..., иП)т обозначим неизвестные узловые параметры и^ . 3) Вычисление решений на слое по времени. Разложение и л имеет вид
п
и л(х) = £ ик (х) + 7/(х), 7/ = (7/ + 7/—1)/2, (27)
к=1
и пусть ио = (и1,..., ип)т - вектор его неизвестных узловых параметров. Определим также вектор-функцию р : Дп ^ Дп с компонентами
п
^(и) = (Ф^Х) и^ (х) + 27/(х) — ил—1, ил—^ , л, 1 < г < п
к=1
Выбрав Сл. = ¥, г = 1,...,п, и подставив разложение (27) в (24), придем к системе нелинейных алгебраических уравнений
7—2 Мои + Коио + Ро(ио) = Фо . (28)
Здесь индекс 0 указывает, что матрицы и векторы имеют длину п, Мо и Ко -подматрицы матрицы масс М и матрицы жесткости К соответственно,
Мо = Ко = {(у*,
Фо = {(/,^)л - Г-2 (^+ (V7/, V*) Л}П=1.
Сделаем следующие замечания.
a) Матрицы М и К стандартной процедурой сборки матриц, принятой в МКЭ, вычисляются лишь один раз в самом начале вычислений.
b) Вектор Фо можно вычислить следующим образом. Сначала стандартной процедурой сборки, принятой в МКЭ, вычислим вектор д = {(/, . Затем определим нулевой вектор-столбец д длины N и его г-й элемент, г = п + 1,..., N, положим равным ) + 7(а^^_1))/2; наконец, вычислим вектор Фо = д — г-2 Мд — Кд.
c) Для вычисления вектора Ро(и) при заданном и = (и1,...,Цп)т сначала вычислим функцию (точнее, ее узловые параметры)
п
м^(х) = 2^ ий ^(х) + 27/(х) — и-1.
к=1
Поскольку Ф[и, V = (Ф(и) — Ф(^))/(и — «), то вектор Ро(и) с компонентами
Р>;.(С7) = (з, 1 < 1' < д =-—
мь —
можно вычислить стандартной процедурой сборки, принятой в МКЭ.
4. Вычислительные эксперименты
В этом разделе мы представляем некоторые результаты вычислений по построенной нами двухслойной схеме для одномерных и двумерных уравнений Клейна-Гордона с известными точными решениями. Вычисления были организованы с целью определения максимального порядка точности схемы на гладких решениях в различных нормах. Для каждой тестовой задачи вычисления были организованы следующим образом.
1. При d =1 отрезок ^ был разбит на п конечных элементов (подотрез-ков) равной длины Н. В двумерном случае разбиение квадратной области ^ на треугольные конечные элементы было получено на основе равномерной ортогональной сетки из (п +1) х (п +1) узлов после деления одним и тем же способом каждого элементарного квадрата со стороной Н на два треугольника. Для вычисления интегралов при реализации схемы МКЭ использовались квадратурные формулы, обеспечивающие диагональность матрицы масс.
2. Функции в начальных условиях дискретной схемы выбирались равными У -интерполянтам мо и Уо. Для решения системы нелинейных алгебраических уравнений (28) на каждом слое времени был использован итерационный метод Ньютона с начальным приближением с предыдущего слоя. Во всех рассмотренных случаях метод сходился квадратично. Итерации завершались следующим образом: после того как равномерная норма разности соседних итераций становилась меньше 10-8, выполнялась еще одна итерация. При таком условии итерационный метод не оказывал влияния на погрешность схемы.
3. Приближенный метод зависит от двух параметров h и т. Учитывая приближения, использованные при построении схемы, естественно ожидать, что погрешность решения в среднеквадратичной норме имеет второй порядок малости как по h, так и по т, т. е. справедлива оценка
L(u) = ||u - uhT< C (h2 + т2). (29)
Для экспериментальной проверки этой оценки были проведены следующие вычисления.
Для набора точек сетки n = n\,U2,... (в одном направлении) и соответствующих им шагов сетки h = h\,h2,... при т = Th были найдены приближенные решения задачи с известным точным решением и вычислены величины L(u) и L(v). Одновременно вычислялись погрешности C(u) и C(v), где
C(и) = max max \u(ai,tj)— uhT(ai)\. (30)
1<i<N 1<j<M J J
Для вычисления интегралов в определениях L(u) и L(v) использовалась составная квадратурная формула Симпсона. Предположив, что Eh к C ha, где Eh -любая из указанных выше норм, и что постоянная C слабо зависит от h , нашли значения ai по формуле
ai = log (Ehi+1 /Ehi)/ log (hi+i/hi) .
Отметим, что значения ai к 2 соответствуют ожидаемому значению a = 2.
4. Вычислениями было подтверждено, что предложенная схема сохраняет полную энергию (с точностью до ошибок округления). Как в одномерном, так и двумерном случаях был решен ряд задач при различных шагах сетки с исходными данными, удовлетворяющими условиям теоремы 1. По найденным решениям была вычислена вариация функционала энергии
AE = max Eh(vhT(t),uhT(t)) - min Eh(vhT(t),uhT(t)).
Во всех случаях она имела порядок 10-16 ~ 10-13 в зависимости от шага h, что свидетельствует о постоянности функционала энергии.
4.1. Решение одномерных задач. Приведем решения двух тестовых задач.
Тест 1. Рассмотрим одномерное уравнение Клейна-Гордона (1) на отрезке [-10,10] при T = 20 с начальными условиями
u(x, 0) = 0, u'(x, 0) = 4 sech(x)
при f(u) = sin u, имеющее точное решение [15]
u(x,t) = 4 arctg(sech(x)t).
На рис. 1 представлен график приближенного решения и его погрешности на грубой сетке узлов. В таблице 1 указана зависимость от h норм погрешности (29), (30), а также соответствующих им значений a , при т = h .
Тест 2. Рассмотрим при d = 1 уравнение Клейна-Гордона (1) на отрезке [-15,45] при f(u) = —u + u3 и T = 60, имеющее точное решение
u(x,t) = tgh((x - ct)/d), с = 0.5, d= >/2(1 - с2).
Рис. 1. Тест 1: приближенное решение (слева) и его погрешность на сетке узлов с шагами Н = 0.25, т = 0.25
Табл. 1
Нормы погрешности и порядки сходимости решения теста 1 при т = Н
h Ци) OíL(u) ад OiL(v) С {и) Ote (и) C{v) <*C(v)
1/100 5.53 е-1 — 1.22 е-1 1.33 е-1 — 2.29 е-2 —
1/200 1.45 е-1 1.93 3.17е-2 1.95 3.54 е-2 1.91 6.13 е-3 1.90
1/400 3.67 е-2 1.98 7.98 е-3 1.99 8.97 е-3 1.98 1.58 е-3 1.96
1/600 1.64 е-2 1.99 3.55 е-3 2.00 3.99 е-3 1.99 6.97 е-4 2.01
1/800 9.20 е-3 2.00 2.00 е-3 2.00 2.25 е-3 2.00 3.98 е-4 2.01
1/1000 5.89 е-3 2.00 1.28 е-3 2.00 1.44 е-3 2.00 2.50 е-4 2.00
Табл. 2
Нормы погрешности и порядки сходимости решения теста 2 при т = Н
h Ци) ад aL(v) С(и) ®C(u) C(v) <*сы
1/100 6.31 е+0 — 2.15 е+0 1.05 е+0 — 3.05 е-1
1/200 1.30 е+0 2.28 4.79 е-1 2.17 2.28 е-1 2.21 7.88 е-2 1.95
1/400 3.18 е-1 2.03 1.17 е-1 2.03 5.55 е-2 2.04 1.86 е-2 2.08
1/600 1.41 е-1 2.01 5.19 е-2 2.01 2.46 е-2 2.01 8.17 е-3 2.03
1/800 7.90 е-2 2.00 2.92 е-2 2.00 1.38 е-2 2.01 4.63 е-3 1.98
1/1000 5.05 е-2 2.00 1.87 е-2 2.00 8.83 е-3 2.00 2.97 е-3 1.98
В таблице 2 указана зависимость от h норм погрешности (29), (30), а также соответствующих им значений а при т = h.
Тест 3. Рассмотрим при d =2 уравнение Клейна-Гордона (1) на отрезке [-7, 7] при ^(u) = sin u и T = 7, имеющее точное решение
u(x, t) = 4atan(exp(x + y — t)).
На рис. 2 представлен график приближенного решения и его погрешности на грубой сетке узлов. В таблице 3 указана зависимость от h норм погрешности (29), (30), а также соответствующих им значений а при т = h.
Тест 4• Рассмотрим при d = 2 уравнение Клейна-Гордона (1) на отрезке [—10,10] при ^(u) = —u + u3 и T =10, имеющее точное решение [16]
и{х, t) = tgh{{x - ct)/d), с = 0.5, d = л/2(2 - с2).
В таблице 4 указана зависимость от h норм погрешности (29), (30), а также соответствующих им значений а при т = h .
Рис. 2. Тест 3: приближенное решение (слева) и его погрешность на грубой сетке узлов с шагами К = 0.35, т = К/2
Табл. 3
Нормы погрешности и порядки сходимости решения теста 3 при т = К/2
к Цч) ®Ь(и) ад ®Ь(1>) С (и) &С(и) С(у)
1/40 7.58 е-01 — 7.47 е-01 — 1.18 е-01 — 1.85 е-01 —
1/80 1.84 е-01 2.05 1.82 е-01 2.03 2.89 е-02 2.03 3.92 е-02 2.24
1/100 1.17е-01 2.02 1.17е-01 2.01 1.84 е-02 2.02 2.48 е-02 2.05
1/140 5.95 е-02 2.01 5.94 е-02 2.00 9.39 е-03 2.00 1.26 е-02 2.01
Табл. 4
Нормы погрешности и порядки сходимости решения теста 4 при т = К/2
к ад ®Ь(и) ад ®Ь(1>) С (и) &С(и) С(у)
1/40 3.36 е-1 — 1.70 е-1 — 2.41 е-2 — 1.33 е-2 —
1/80 8.34 е-2 2.01 4.32 е-2 1.98 5.84 е-3 2.05 3.41 е-3 1.96
1/100 5.34 е-2 2.00 2.77 е-2 1.99 3.73 е-3 2.01 2.20 е-3 1.97
1/140 2.72 е-2 2.00 1.42 е-2 1.99 1.90 е-3 2.01 1.13е-3 1.97
Из результатов вычислений, представленных в таблицах 1-4, следует, что предложенная в работе двухслойная схема позволяет определить как решение задачи и, так и его производную по времени V = и' с погрешностью порядка О (к2 + т2) в среднеквадратической норме Ь2 (ф).
Заключение
Предложено семейство сеточных схем Петрова-Галеркина на основе метода конечных элементов для решения нелинейного уравнения Клейна-Гордона. Дискретные схемы определены в терминах решения задачи и его производной по переменной времени. Показано, что они обеспечивают сохранение полной энергии на дискретном уровне. Проведено численное исследование простейшей двухслойной схемы подобного типа на основе линейных треугольных конечных элементов. Численным решением ряда тестовых задач, имеющих гладкие решения, показано, что она позволяет определить как решение задачи, так и его производную по времени с погрешностью порядка О (к2 + т2) в среднеквадратической норм.
Литература
1. Додд Р., Эйлбек Дж., Гиббон Дж., Моррис Х.М. Солитоны и нелинейные волновые уравнения. М.: Мир, 1988. 694 с.
2. Браун О.М., Кившарь Ю. С. Модель Френкеля-Конторовой: Концепции, методы, приложения. М.: Физматлит, 2008. 519 с.
3. Bratsos A.G. A modified predictor-corrector scheme for the two-dimensional sine-Gordon equation // Numer. Algorithms. 2006. V. 43, No 4. P. 295-308. https://doi.org/10.1007/s11075-006-9061-3.
4. Bratsos A. G. The solution of the two-dimensional sine-Gordon equation using the method of lines // J. Comput. Appl. Math. 2007. V. 206, No 1. P. 251-277. https://doi.org/10.1016/jxam.2006.07.002.
5. Bao W., Dong X. Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime // Numer. Math. 2012. V. 120, No 2. P. 189229. https://doi.org/10.1007/s00211-011-0411-2.
6. Wang L., Chen W., Wang C. An energy-conserving second order numerical scheme for nonlinear hyperbolic equation with an exponential nonlinear term //J. Comput. Appl. Math. 2015. V. 280. P. 347-366. https://doi.org/10.1016/j.cam.2014.11.043.
7. Kang X., Feng W., Cheng К., Guo C. An efficient finite difference scheme for the 2d sine-Gordon equation // J. Nonlinear Sci. Appl. 2017. V. 10, No 6. P. 2998-3012. http://dx.doi.org/10.22436/jnsa.010.06.14.
8. Argyris J., Haase M., Heinrich J.C. Finite element approximation to two dimensional sine-Gordon solitons // Comput. Methods Appl. Mech. Eng. 1991. V. 86, No 1. P. 1-26. https://doi.org/10.1016/0045-7825(91)90136-T.
9. Wang Q., Cheng D. Numerical solution of damped nonlinear Klein-Gordon equations using variational method and finite element approach // Appl. Math. Comput. 2005. V. 162, No 1. P. 381-401. https://doi.org/10.1016/j.amc.2003.12.102.
10. Asgari Z., Hosseini S.M. Numerical solution of two-dimensional sine-Gordon and MBE models using Fourier spectral and high order explicit time stepping methods // Comput. Phys. Commun. 2013. V. 184, No 3. P. 565-572. https://doi.org/10.1016/j.cpc.2012.10.009.
11. Jiwari R., Pandit S., Mittal R.C. Numerical simulation of two-dimensional sine-Gordon solitons by differential quadrature method // Comput. Phys. Commun. 2012. V. 183, No 3. P. 600-616. https://doi.org/10.1016/j.cpc.2011.12.004.
12. Dehghan M., Shokri A. Numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions // Math. Comput. Simul. 2008. V. 79, No 3. P. 700-715. https://doi.org/10.1016/j.matcom.2008.04.018.
13. Browder F.E. On non-linear wave equations // Math. Z. 1962. V. 80, No 3. P. 249-264.
14. Лионе Ж.Л. Некоторые методы решения нелинейных краевых задач. М: Мир, 1972. 587 с.
15. Dehghan M., Mirzaei D. The boundary integral equation approach for numerical solution of the one-dimensional sine-Gordon equation // Numer. Methods Partial Differ. Equations. 2008. V. 24, No 6. P. 1405-1415. https://doi.org/10.1002/num.20325.
16. Khan K., Akbar M.A. Exact solutions of the (2+1)-dimensional cubic Klein-Gordon equation and the (3+1)-dimensional Zakharov-Kuznetsov equation using the modified simple equation method // J. Assoc. Arab Univ. Basic Appl. Sci. 2014. V. 15, No 1. P. 74-81. https://doi.org/10.1016/j.jaubas.2013.05.001.
Поступила в редакцию 10.07.2023 Принята к публикации 04.09.2023
Даутов Рафаил Замилович, доктор физико-математических наук, профессор
Институт вычислительной математики и информационных технологий, Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: [email protected] Салимзянова Гулина Ринатовна, студентка 2 курса магистратуры
Институт вычислительной математики и информационных технологий, Казанский (Приволжский) федеральный университет
ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: [email protected]
ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)
UCHENYE ZAPISKI KAZANSKOGO UNIVERSITETA. SERIYA FIZIKO-MATEMATICHESKIE NAUKI (Proceedings of Kazan University. Physics and Mathematics Series)
2023, vol. 165, no. 3, pp. 190-207
ORIGINAL ARTICLE
doi: 10.26907/2541-7746.2023.3.190-207
A Conservative Fully Discrete Finite Element Scheme for the Nonlinear Klein—Gordon Equation
R.Z. Dautov*, G.R. Salimzyanova**
Kazan Federal University, Kazan, 420008 Russia E-mail: *[email protected], **[email protected]
Received July 10, 2023; Accepted September 04, 2023 Abstract
This article proposes a family of the Petrov-Galerkin-FEM methods that can be used to solve the nonlinear Klein-Gordon equation. The discrete schemes were formulated based on the solution of the problem and its time derivative. They ensure that the total energy is conserved at a discrete level. The simplest two-layer scheme was studied numerically. Based on the solution of the test problems with smooth solutions, it was shown that the scheme can determine the solution of the problem, as well as its time derivative with an error of the order of O(h2 + t2) in the continuous L2 norm, where t and h characterize the grid steps in time and space, respectively.
Keywords: Petrov-Galerkin method, finite element method, Klein-Gordon equation, implicit scheme
Figure Captions
Fig. 1. Test 1: approximate solution (on the left) and its error on the grid of nodes with steps h = 0.25, t = 0.25.
Fig. 2. Test 3: approximate solution (on the left) and its error on the coarse grid of nodes with steps h = 0.35, t = h/2.
References
1. Dodd R., Eilbeck J., Gibbon J., Morris H.M. Solitony i nelineinye volnovye uravneniya [Solitons and Nonlinear Wave Equations]. Moscow, Mir, 1988. 694 p. (In Russian)
2. Braun O.M., Kivshar Yu.S. Frenkel-Kontorova Model: Concepts, Methods, Applications. Ser.: Theoretical and mathematical physics. Berlin, Heidelberg, Springer. xviii, 472 p. https://doi.org/10.1007/978-3-662-10331-9.
3. Bratsos A.G. A modified predictor-corrector scheme for the two-dimensional sine-Gordon equation. Numer. Algorithms, 2006, vol. 43, no. 4, pp. 295-308. https://doi.org/10.1007/s11075-006-9061-3.
4. Bratsos A.G. The solution of the two-dimensional sine-Gordon equation using the method of lines. J. Comput. Appl. Math., 2007, vol. 206, no. 1, pp. 251-277. https://doi.org/10.1016/jxam.2006.07.002.
5. Bao W., Dong X. Analysis and comparison of numerical methods for the Klein-Gordon equation in the nonrelativistic limit regime. Numer. Math., 2012, vol. 120, no. 2, pp. 189229. https://doi.org/10.1007/s00211-011-0411-2.
6. Wang L., Chen W., Wang C. An energy-conserving second order numerical scheme for nonlinear hyperbolic equation with an exponential nonlinear term. J. Comput. Appl. Math., 2015, vol. 280, pp. 347-366. https://doi.org/10.1016/jxam.2014.11.043.
7. Kang X., Feng W., Cheng K., Guo C. An efficient finite difference scheme for the 2d sine-Gordon equation. J. Nonlinear Sci. Appl., 2017, vol. 10, no. 6, pp. 2998-3012. http://dx.doi.org/10.22436/jnsa.010.06.14.
8. Argyris J., Haase M., Heinrich J.C. Finite element approximation to two dimensional sine-Gordon solitons. Comput. Methods Appl. Mech. Eng., 1991, vol. 86, no. 1, pp. 1-26. https://doi.org/10.1016/0045-7825(91)90136-T.
9. Wang Q., Cheng D. Numerical solution of damped nonlinear Klein-Gordon equations using variational method and finite element approach. Appl. Math. Comput., 2005, vol. 162, no. 1, pp. 381-401. https://doi.org/10.1016/j.amc.2003.12.102.
10. Asgari Z., Hosseini S.M. Numerical solution of two-dimensional sine-Gordon and MBE models using Fourier spectral and high order explicit time stepping methods. Comput. Phys. Commun., 2013, vol. 184, no. 3, pp. 565-572. https://doi.org/10.1016/j.cpc.2012.10.009.
11. Jiwari R., Pandit S., Mittal R.C. Numerical simulation of two-dimensional sine-Gordon solitons by differential quadrature method. Comput. Phys. Commun., 2012, vol. 183, no. 3, pp. 600-616. https://doi.org/10.1016/j.cpc.2011.12.004.
12. Dehghan M., Shokri A. Numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions. Math. Comput. Simul., 2008, vol. 79, no. 3, pp. 700-715. https://doi.org/10.1016/j.matcom.2008.04.018.
13. Browder F.E. On non-linear wave equations. Math. Z., 1962, vol. 80, no. 3, pp. 249-264.
14. Lions J.L. Nekotorye metody resheniya nelineinykh kraevykh zadach [Some Methods for Solving Non-Linear Boundary Value Problems]. Moscow, Mir, 1972. 587 p. (In Russian)
15. Dehghan M., Shokri A. Numerical method for solution of the two-dimensional sine-Gordon equation using the radial basis functions. Math. Comput. Simul., 2008, vol. 79, no. 3, pp. 700-715. https://doi.org/10.1016/j.matcom.2008.04.018.
16. Khan K., Akbar M.A. Exact solutions of the (2+1)-dimensional cubic Klein-Gordon equation and the (3+1)-dimensional Zakharov-Kuznetsov equation using the modified simple equation method. J. Assoc. Arab Univ. Basic Appl. Sci., 2014, vol. 15, no. 1, pp. 74-81. https://doi.org/10.1016/j.jaubas.2013.05.001.
Для цитирования: Даутов Р.З., Салимзянова Г.Р. Консервативная полностью / дискретная схема МКЭ для нелинейного уравнения Клейна-Гордона // Учен. \ зап. Казан. ун-та. Сер. Физ.-матем. науки. 2023. Т. 165, кн. 3. С. 190-207. URL: https//doi.org/10.26907/2541-7746.2023.3.190-207.
For citation: Dautov R.Z., Salimzyanova G.R. A conservative fully discrete finite element / scheme for the nonlinear Klein-Gordon equation. Uchenye Zapiski Kazanskogo Univer-\ siteta. Seriya Fiziko-Matematicheskie Nauki, 2023, vol. 165, no. 3, pp. 190-207. URL: https//doi.org/10.26907/2541-7746.2023.3.190-207. (In Russian)
