КОНСЕРВАТИВНАЯ СХЕМА МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ УРАВНЕНИЯ КИРХГОФА Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА. СЕРИЯ ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

2023, Т. 165, кн. 2 С. 115-131

ISSN 2541-7746 (Print) ISSN 2500-2198 (Online)

ОРИГИНАЛЬНАЯ СТАТЬЯ УДК 519.63

doi: 10.26907/2541-7746.2023.2.115-131

КОНСЕРВАТИВНАЯ СХЕМА МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ УРАВНЕНИЯ КИРХГОФА

Казанский (Приволжский) федеральный университет, г. Казань, 420008, Россия

Предложена неявная двухслойная схема метода конечных элементов для решения уравнения Кирхгофа, которое представляет собой нелинейное нелокальное уравнение гиперболического типа и включает интеграл Дирихле. Дискретная схема сформулирована в терминах решения задачи и его производной по переменной времени и обеспечивает сохранение полной энергии на дискретном уровне. Показано, что решение схемы на слое по времени может быть эффективно получено методом Ньютона, несмотря на нелокальность уравнения. На основе решения тестовых задач с гладкими решениями установлено, что схема позволяет определить как решение задачи, так и его производную по времени с погрешностью порядка 0(Н2 + т2) в среднеквадратической норме, где т и !г характеризуют шаги сетки по времени и пространству соответственно.

Ключевые слова: уравнение Кирхгофа, метод конечных элементов, метод Петрова-Галеркина, неявная схема, метод Ньютона

известное как уравнение Кирхгофа, при d = 1 и краевом условии = 0 было получено в 1876 г. Кирхгофом [1] как нелинейная модель поперечных колебаний струны, закрепленной в концевых точках. Уравнение Кирхгофа изучалось Берн-штейном [2] и Похожаевым [3, 4] в специальном классе аналитических функций, см. также статью Лионса [5].

В последние годы большое внимание уделяется изучению нелокальных уравнений различных типов, включающих интеграл Дирихле. Интерес к таким задачам связан как с моделированием разнообразных физических и биологических явлений, так и с рядом интересных математических особенностей этих задач. В связи с этим отметим статью [6], посвященную изучению корректности уравнения (1), обзорную работу [7] и обзор литературы в [8], а также статьи [9-12], в которых рассматривались уравнения, более общие, чем (1).

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

Р. З. Даутов, M. В. Иванова

Аннотация

Введение

1. Следующее интегро-дифференциальное уравнение

(1)

численные методы. Однако существенная особенность уравнения (1) в виде интегрального коэффициента приводит к необходимости его отдельного рассмотрения. Численному решению таких задач посвящено лишь небольшое число работ. Среди известных нам отметим: статьи [13,14], посвященные соответственно обоснованию конформного и смешанного методов конечных элементов для решения эллиптического уравнения, соответствующего (1), и статьи [15-17], посвященные исследованиям неявных дискретных методов решения уравнений параболического типа. Еще меньше работ посвящено решению гиперболического уравнения (1) [18-20]. В статье [18] предложен спектральный метод Галеркина для решения одномерной задачи Кирхгофа (1) и получены априорные оценки погрешности решения модифицированной разностной схемы Кранка-Николсона. В статье [19] также для решения одномерной задачи (1), предварительно сведенной к системе уравнений первого порядка, построена и численно исследована консервативная конечно-разностная схема. Для решения двумерного уравнения (1) с затуханием в [20] предложен и исследован неконформный метод конечных элементов, основанный на четырехугольных квадратичных элементах.

Поскольку уравнение Кирхгофа описывает консервативные системы, следует отметить, что в практических вычислениях консервативные численные схемы его решения являются более предпочтительными, чем неконсервативные: при прочих равных условиях они более точны. Ключевым моментом является то, что такие схемы сохраняют некоторые важные инвариантные свойства дифференциального уравнения и более точно воспроизводят особенности решения при вычислениях на больших отрезках времени. Для решения уравнения Кирхгофа нам известна только консервативная конечно-разностная схема для одномерной задачи [19].

3. Основной целью настоящей работы является численное исследование точности предлагаемого нами двухслойного консервативного неявного метода конечных элементов решения уравнения (1) в произвольной области ^. Первоначально введением новой неизвестной V = ди/дЬ уравнение сводится к системе уравнений, которая далее аппроксимируется на основе комбинации метода конечных элементов (по пространственным переменным) и метода Петрова-Галеркина (по временной переменной). На основе решения ряда тестовых задач с гладкими решениями показано, что такая схема позволяет определить как решение задачи, так и его производную по времени с погрешностью порядка 0(Н2 + т2) в среднеквадрати-ческой норме, где т и Н характеризуют шаги сетки по времени и пространству соответственно.

1. Постановка задачи

1.1. Исходная задача. Пусть Т > 0, ^ - ограниченная область в Кй с липшицевой границей Г, й > 1. Обозначим через Qт = ^ х (0, Т] цилиндр в Ял+1 с боковой поверхностью Гт = Г х (0, Т]. Через и' далее будем обозначать производную по времени Ь функции и(х,Ь), определенной в области Q. Рассмотрим уравнение Кирхгофа

и''(х,Ь) - (1 + ||Уи||2) Ди(х,Ь) = /(х,Ь), (х,Ь) € Qт, (2)

при краевых условиях Дирихле1:

и(х,Ь) = 7(х,Ь), (х,Ь) € Гт. (3)

хДля удобства изложения будем считать, что функция 7 определена в Qт .

В начальный момент времени 4 = 0 решение уравнения и его производная по времени считаются известными:

и(х, 0) = ио(х), и'(х, 0) = «о(х), х € О. (4)

Здесь Уи означает градиент функции и,

|а| = (а • а)1/2 есть длина вектора а, а • Ь = ^4 аф^ - скалярное произведение векторов а и Ь.

Отметим одно важное свойство рассматриваемой задачи (см., например, [5]). Пусть данные задачи 7 и f не зависят от времени, т. е. 7 = 7(х) и f = f (х). Тогда в силу краевого условия (3) имеем и'(х,4) =0 на Гт. Умножим уравнение (2) на и'(х, 4) и проинтегрируем полученное равенство по области О. В результате получим

— — [ (и'(х, ¿))2 йх+ (1 + ||Ум||2) [ \7и(х,1)-\7и'(х, = — ( ¡{х)и{х^)с1х, (5) 2 о ио & и о

поскольку согласно формуле Остроградского-Гаусса

— / Ди(х,4) С(х,4) «х = Уи(х, 4) •УС (х, 4) ¿х, (6)

оо

если ^ обращается в нуль на Г. Определим функцию П(в) = в + в2/2. Теперь легко видеть, что равенству (5) можно придать вид

= г е (о,т], (7)

f (х)и(х,

2 о о

где функционал (полной энергии) Е определяется следующим образом:

1

оо а и(4) означает функцию и в момент времени 4. Из (7) следует равенство

Е(и(£),и'(£)) = Е(и0,«0) для всех 4 € [0,Т],

которое принято называть законом сохранения (полной) энергии.

Для определения дискретного метода решения сформулированной задачи (2)-(4) нам понадобится его обобщенная формулировка. С этой целью введем дополнительные обозначения и определим пространства функций.

1.2. Пространства функций. Используем стандартное обозначение Ь2(О) для пространства функций, измеримых в области О и интегрируемых в ней с квадратом. Через Н 1(О) обозначим пространство Соболева первого порядка, через НО (О) = {и € Н 1(О) : и|г = 0} - его подпространство. Для сокращения записи положим

Н = £2(О), V = Н 1(О), V0 = Н1(О). Через (•, •) обозначим как скалярное произведение в Н, т.е.

(и = /

о

и V ¿х,

о

так и отношение двойственности между пространством V0 и его сопряженным H-1(Q). Пусть также

(Vu, Vv) = / Vu -Vv dx. JQ

Далее для заданной функции u переменных x G Q и t G [0, T] под u(t), t G [0, T] будем понимать функцию x ^ u(x, t), x G Q. Подобные функции, определенные на [0,T] со значениями в некотором банаховом пространстве B функций от x, будем рассматривать как элементы пространства Lp(0, T; B). Напомним, что Lp(0, T; B) есть пространство (классов) измеримых функций [0, T] ^ B таких, что

ifT \ 1/р

IM|lp(0,t;B) = (У0 \\u(t)\\pB dtj < Ж, 1 < p< ж.

При p = ж норма в нем определяется функционалом ess sup te(0,T)l|u(t)||B (подробнее см., например, [13, с. 469]). Отметим, что Lp(0,T; Lp(Q)) = Lp(Q).

1.3. Обобщенное решение задачи. Относительно исходных данных задачи будем предполагать, что выполнены следующие условия:

(H) f G L2(0,T; H), u0 G V, v0 G H, y является следом на Г х [0,T] некоторой функции ц из L^(0, T; V) такой, что ц' G L^(0, T; H).

Обобщенное решение задачи (2)-(4) будем искать в классе функций, гарантирующих конечность функционала энергии E(u(t), u'(t)) почти всюду на (0, T). Это приводит к ограничениям u G L^(0, T; V), u' G L^(0, T; H).

Введем новую неизвестную v = u' и запишем исходную задачу в виде системы

v'(t) - (1 + ||Vu||2) Au(t) = f (t), (8)

u'(t)= v(t), (9)

u - y G L^(0,T; V0), u(0) = щ, v(0) = v0,

где t G (0, T). Для определения обобщенного решения задачи умножим обе части уравнения (8) на произвольную функцию Z G L1(0,T; V0), проинтегрируем полученное равенство по области Qt и учтем равенство (6). Аналогично умножим равенство (9) на произвольную функцию ц G L1(0,T; H). В результате придем к следующим тождествам:

,-т т ,-т

(v'(t),Z(t)) dt +i a(u(t),Z(t)) dt =f (f (t),Z(t)) dt, (10)

00

(u'(t),n(t)) dt = j (v(t),n(t)) dt. (11)

и о и 0

Здесь использовано обозначение

а(и,() = (1 + \\Чи\\2) (Чи, Ч(). Определение 1. Пару функций и и V, удовлетворяющих условиям

и - 7 е ь^(0,Т; V0), v,u' е ь^(0,Т; н), V' е ь^(0,Т; V*),

и(0) = ио, v(0) = vо,

а также тождествам (10) и (11) для всех С е Ь^(0, Т; V0) и ц е Ь1(0, Т; Н), назовем обобщенным решением задачи (2)—(4).

Заметим, что при выполнении условий (Н1), (Н2) данное определение корректно в том смысле, что все слагаемые в (10), (11) являются конечными, а начальные условия имеют смысл (см., например, обсуждение обобщенной формулировки гиперболической задачи в [21, Замечание 1.1, с. 21]).

Далее будем предполагать, что определенное выше обобщенное решение существует (это так по крайней мере для достаточно гладких исходных данных ио, vо, f [5]), и сосредоточимся на его приближенном определении.

2. Определение неявной схемы МКЭ

Укажем способ дискретизации поставленной выше задачи, обеспечивающий сохранение энергии на дискретном уровне. Его можно рассматривать как метод Петрова-Галеркина, основанный на аппроксимациях, принятых в методе конечных элементов. Полностью дискретную схему получим последовательной дискретизацией задачи сначала по пространственным переменным, а затем - по переменной времени. Для облегчения изложения далее будем предполагать, что О является многогранником в .

2.1. Полудискретная схема МКЭ. Пусть Н > 0 означает максимальный размер конечных элементов (симплексов в ), Т/ = {К., К2,..., К^} - множество конечных элементов (отрезков при « = 1, треугольников при « = 2, тетраэдров при « = 3), образующих, как это принято в МКЭ, точное разбиение (триангуляцию) области О. Множество всех вершин симплексов из Т/ обозначим ш/ и назовем сеткой узлов на О. Пусть = {а^}|=1, N = N(Н). Далее для удобства изложения предположим, что внутренние узлы сетки имеют номера от 1 до п, а граничные (лежащие на Г) - от п + 1 до N.

На основе триангуляции Т/ определим пространство конечных элементов

Vh = К € С (О) : vh|к € Р V К € Т/},

где Р1 есть множество алгебраических многочленов первой степени, т. е. функций вида а0 + а1х1 + ... + а^х^. Определим также его подпространство

= К € V, : Vh|г = 0}.

Пространства Vh и Vh° будем рассматривать как аппроксимации пространств Н 1(О) и Нд (О) соответственно. Отметим, что Vh С Н 1(О), У^ С Н(О).

Через {^¿(х)}|=1 обозначим базис Лагранжа в V/,,, т. е. множество таких функций € Vh, что ^¿(а^) = ^, г, ] = 1,...,^ где - символ Кронекера. В этом базисе произвольная функция vh € Vh имеет представление

N

Vh(x) = VI ^¿(х), (12)

¿=1

где v¿ = vh(a¿), г = 1,...,N. Вектор-столбец коэффициентов V = (V!, v2,..., vN)т принято называть вектором узловых параметров функции Vh € Vh. Будем использовать запись V ^ Vh. Функцию и/ из V;,, назовем интерполянтом функции и (Е С(Г2), если и1(щ) = 11(0-1), г = 1,..., N.

Далее через uоh и vоh обозначим функции из Vh, являющиеся некоторыми аппроксимациями функций ио и vо соответственно. Через 7/ обозначим функцию [0,Т] ^ Vh такую, что 7/ (4) есть некоторая аппроксимация 7(4). Например, их можно определить как наилучшие среднеквадратические приближения или как интерполянты соответствующих функций.

Определение 2. Функции ин,ун из Ь^(0,Т; Vн) такие, что

пн — ^н € Ь^(0,Т; VI0), € Ь^(0,Т; Vh),

удовлетворяющие условиям ин(0) = и0н, ун(0) = у0н, а также тождествам /• т Г т ,-Т

(у'к(к),ф)) ¿к +[ а(ин(к),ф)) ¿к =[ (I(к),Ф)) ¿к, Зн Jo

(и'(к),п(к)) ¿к = I (ун(к), п(к)) ¿к

для всех ( € Ьг(0, Т; V'£), п € Ьг(0, Т; V') назовем приближенным решением задачи (2)—(4) по методу конечных элементов.

2.2. Неявная схема. Пусть т > 0. На отрезке [0, Т] определим сетку узлов шт = {ко,3 = 0,1,...,М} с шагами то = ко — , т3 < т. Через Б г(0,Т; V') обозначим множество всех функций, непрерывных на [0, Т] и линейных на каждом отрезке ] со значениями в Vн. Функция инт из Б 1(0,Т; V') имеет пред-

ставление

инт(к) = г4-1 + (и{ - «Г1) г е \tj-uhl (13)

тз

где и' = инт (ко), 3 =0,1,...,М.

Определим также пространство Б(0,Т; Vн) как множество функций £ : [0,Т] ^ V'н, являющихся постоянными на каждом отрезке \bj-i, к3]. Отметим, что обобщенная производная функции из Бг(0,Т; Ун) принадлежит Б(0,Т; Ун).

Искомую дискретную схему для задачи (2)-(4) определим следующим образом:

Определение 3. Найти функции инт,унт из Б 1(0,Т; Уh) такие, что инт — 1нт € бг(0,Т; V0), инт(0) = и0н, Унт(0) = УнН, которые удовлетворяют тождествам

[ (Унт(к), ((к)) ¿к +( а(инт(к),((к)) ¿к = [ (I(к),((к)) ¿к, (14)

0 0 0

I (и^(к), п(к)) ¿к = [ (унт(к), п(к)) ¿к (15)

00

для всех С € б(0, Т; ^), п € Б(0, Т; Уh).

Рассмотрим реализацию этой схемы. Пусть в момент времени ко-1 решение схемы, а именно, пара функций У' 1 = инт(к^-1), у' 1 = унт(к3-1), известно. Это так при 3 = 1, поскольку в силу начальных условий имеем у' = инн, Ун = инн. Выберем в тождествах (14) и (15) функции £ (к) и п(к) равными нулю всюду, кроме отрезка [ко- 1,к3], на котором примем их равными соответственно произвольно зафиксированным € Уh0 и щ € Уh. Тогда придем к следующим равенствам:

ГЬз ГЬз ГЬз

/ (у'кт(к),Сн) ¿к + а(инт(к), Сн) ¿к = (I(к),Сн) ¿к, (16)

/■ Г Ьз

/ «т(к),пн) ¿к = (унт(к), пн) ¿к. (17)

3 — 1 3 — 1

00

Воспользуемся здесь формулой (13) для представления мЛт(4) и аналогичной формулой для (4). Учитывая, что

/ а(м^т(¿), Сл) ^ = тЛ «(^Ч + (1 - а)Ч—1, Сл) ./0

после вычисления остальных интегралов в (16) и (17) придем к следующим уравнениям, справедливым для всех з = 1,...,М:

о — ¿-1 /• 1

Гн ,0) + / а{аи{ + (1 - а)4~\О) <Ь = (/Л О), (18)

4 то у ./0

= (19)

о

Здесь использовано обозначение

Я = / / Ио + (1 - ^Хо-1)

0

В силу произвольности € из (19) следует равенство

о 1 о I 1

ч - ч <+<

то 2

Найдем отсюда «Л и подставим полученное выражение в (18). В результате придем

7

к уравнению для определения искомого неизвестного :

Г* и°к 1 - г{-\о) + Г «(^Ч + (1 - 6.) ^ = (Я, а), (20)

4 то 7 70

1 / <7

где то = то/2, О. - произвольная функция из у0. Уравнение (20) сводится к системе (нелинейных) алгебраических уравнений. После его решения найдем , а потом определим «Л по формуле

о о—1

о 7 — 1

< = % —< ■

' о

2.3. Сохранение энергии на дискретном уровне.

Теорема 1. Пусть / = /(х), 7 = 7(х), мЛт, V Лт - решения схемы (14), (15). Тогда Е(мЛт(¿), (¿)) = Е(м0Л, «0Л) для всех 4 € .

Доказательство. Поскольку 7 не зависит от 4, 7лт совпадает с 7л и также не зависит от 4. Поэтому мЛт (¿) € для всех 4. Это позволяет выбрать в (14) и (15) пробные функции равными

с<«>={млт<()' ;;о, *«={т<«• о,

где з > 1 .В результате придем к следующим равенствам

Г1э Г1э Г1э

/ («Лт (^Чт ^ + / а(мЛт (4),мЛт ^ = (/,Чт^ ^

0 0 0

Г 1з Г 1з

/ (Чт (*), Чт (*)) ^ = («Лт (¿), «Лт (¿)) <Й.

00

1

3

Очевидно, из них вытекает соотношение

Ггз Ггз Ггз

/ (у'Т (к),унт(к)) ¿к + а(инт (к),и'нт(к)) ¿к = / (1,инт(к)) ¿к. (21)

0 0 0

Преобразуем это равенство. Заметим, что

И*)У(*)) = 1/2<*ДЙ(г>(*),г>(*)), <»(«;(*),«/(*)) = | П(||УЦ*)||2).

Поэтому равенство (21) можно записать в следующем виде:

Г 3 ¿

Е(инт (к),унт (к)) ¿к = 0. Следовательно, г'Лт(^)) = Е(ион, и0к), з = \ ,...,М. □

0 ¿к

2.4. Реализация дискретной схемы. Рассмотрим вопросы, возникающие при решении уравнения (20). Для определенности будем считать, что есть ин-терполянт 7. В этом случае условие инт —7нт € Б^(0, Т; Уh0) равносильно условиям <|Г = (кз) для всех з или условиям и3н(аъ) = 7(аъ,к3) для всех п + 1 < г < N, 1 < 3 < М. Поэтому согласно (12) неизвестная функция у' имеет представление

п N

<(х) = £ и Фъ(х) + £ Фъ(x), Ц = 1 (аъ,к3), (22)

1=1 ъ=п+1

где иъ = и"н(аъ). Обозначим через и = (и1,..., ип, 7П+1,..., вектор узловых

параметров и"н, а через и = (и1,... ,и3)т - его неизвестные узловые параметры. Поскольку функции из Уh0 равны нулю на границе Г, то

Сн(х) = £ Съ Фъ(х) VСн € Уh0. 1

Поэтому уравнение (20) эквивалентно системе равенств 3 и3-1 \ г 1

тз 4 тз ' 30

Гк 1'Ъ ~ + / а{<ти{ + (1 - с1а - (Я, = 0, (23)

^ т° / Зп

где 1 < г < п, у' задана формулой (22).

Обозначим через (и) левую часть (23), а через Я (и) - вектор-функцию длины п с компонентами Яъ(и). Тогда решение уравнения (20) равносильно решению системы алгебраических уравнений

Я (и) = 0, и € Еп. (24)

Для реализации итерационных методов решения уравнения (24) необходимы удобные для использования формулы вычисления вектора Я (и) при заданном векторе и. Для их получения определим матрицы М (матрицу масс) и К (матрицу жесткости)

м = {(Фз,Фъ)}&=1, К = {№3, 4<ръ)}%=1.

Эти матрицы могут быть экономно вычислены стандартной процедурой сборки матриц, принятой в МКЭ. Это относится также к вектору (сил) д с компонентами

9ъ = (I3,Фъ) = I3Фъ ¿х, 1 < г < N,

для вычисления которого используются квадратурные формулы. По определению

) = Ми • V, (Уи^, Упь) = Ки • V, (25)

где и ^ и^, П ^ Пь .

1. Вычисление Р(и). Учитывая, что вектор узловых параметров совпадает с ортом г-ой оси в , то согласно первой формуле (25) имеем

1

где и

7_1

7 7 — 1

и! - и

7 — 1 7_1

,7—1

= — [ — М (й - и*-1) - Мг^-1

(26)

1. Под интегралом во втором слагаемом в (23) стоит многочлен третьей степени от а. Поэтому интеграл можно точно вычислить при помощи квадратурной формулы Симпсона. Учтем, что ||Уи^||2 = |Ки|2, (Уи^, = (Ки);, и определим величины

Л^' = 1 + Ки ■ й, Х>+1/'2 = 1 + 2 К{й + и^1) • (й + и*'1). (27)

Теперь нетрудно вычислить

•_ 1 1 • 4

а(аи{ + (1-а)и{ - а(и{, + - а

7 I 7 _ 1 л

Щ + < \ 1 7—1

" " 1 ¥г) + ^ а(иН ,¥><) =

2

А7

-1

А7 А7+1/2

— А'м +-К (и + г^"1)

6 3 1 '

Из формул (23), (26) и (28) следует, что компоненты Р(и) равны

Ки7

-1

. (28)

Р (и) =

— М(й - и*'1) - — Мг^'"Ч т0 т0 т0

А7 А7+1/2 , . + — Ки +---К(й + и1'1)

А7-

6

3

Ки7-1 - ^

где 1 < г < п. Отметим, что Р;(и) зависит от всех компонент и, поскольку величины А7 и зависят от всех компонент и (см. (27)). Следствием этого является полная заполненность матрицы Якоби Р'(и). Это делает непрактичным использование при больших значениях п итерационных методов, требующих вычисления матрицы Якоби (например, метода Ньютона).

2. Переформулировка дискретной задачи. Дадим эквивалентную формулировку дискретной задачи (24). Для этого введем две новые неизвестные А = А7, и расширенный вектор неизвестных V = (и, А, . Расширенная система уравнений для определения V включает следующие п +2 уравнения:

( — М(й - иэ-1) - — Мг^"Ч

6 3 1 ;

А7_

Ки7-1 - # =0, (29)

1 12

(1 + Ки ■ й) - — А = 0,

/I \ 2

(1 + - К (И + и*-1) ■ (И + и*-1)) " | М = 0,

(30)

7

7

7

7

0

6

1

6

1

6

Рис. 1. Приближенное решение тестовой задачи 1 при £ = Т (слева) и его погрешность (справа) на сетке узлов с шагами ^ = 1/40, т = 3Н/8

где А7 1 = 1 + Ки7 • и7 1. Запишем эту систему в виде Ф(и) = 0. Матрица Якоби

этой расширенной системы является симметричной разреженной невырожденной

2

матрицей вида2

Ф' (и) =

А ъ с

ът -1/12 0

ст 0 -2/3

А =

Т7 Т7

1 6 3

К,

Ъ = - Ки, с = - К (и + М-7"1), 6 3

поскольку А является разреженной симметричной положительно-определенной матрицей. Для решения системы (29)-(31) уже можно использовать метод Ньютона. На к -ой итерации он сводится к решению системы линейных алгебраических уравнений вида

р'(и(к))(и(к+1) - и(к)) = -Р(и(к)), к = 0,1,...,

для нахождения приближения и(к+1) к решению и. Начальное приближение можно выбрать (например) с предыдущего слоя по времени.

1

3. Вычислительные эксперименты

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

Функции в начальных условиях дискретной схемы равны интерполянтам ио и г>о. Для решения расширенной системы нелинейных алгебраических уравнений на каждом слое времени использован итерационный метод Ньютона с начальным приближением с предыдущего слоя. Во всех рассмотренных случаях метод сходился квадратично. Итерации завершались следующим образом: после того как равномерная норма разности соседних итераций становилась меньше 10-8, выполнялась еще одна итерация. При таком условии итерационный метод не оказывал влияния на погрешность схемы.

Приближенный метод зависит от двух параметров к и т. Учитывая приближения, использованные при построении схемы, естественно ожидать, что погрешность

2Верхний индекс Т означает операцию транспонирования.

Табл. 1

Нормы погрешности и порядки сходимости решения теста 1 при т = 3Я/8

h Ци) аь(и) L(v) CXL(V) С (и) «С (и) C(v) aG(v)

1/10 1.13 е-1 — 1.14 е-1 — 2.86 е-2 — 5.59 е-2 —

1/20 2.83 е-2 1.99 2.90 е-2 1.98 7.06 е-3 2.02 1.34 е-2 2.06

1/40 7.09 е-3 2.00 7.29 е-3 1.99 1.75 е-3 2.01 3.41 е-3 1.98

1/80 1.77 е-3 2.00 1.83 е-3 2.00 4.40 е-4 1.99 8.51 е-4 2.00

1/160 4.44 е-4 2.00 4.57 е-4 2.00 1.10 е-4 2.00 2.12 е-4 2.01

решения в среднеквадратичной норме имеет второй порядок малости как по h, так и по т, т. е. справедлива оценка

L(u) = ||u - UhT||l2(q) < C (h2 + т2). (32)

Для экспериментальной проверки этой оценки были проведены следующие вычисления: для набора точек сетки n = ni, n2,... (в одном направлении) и соответствующих им шагов сетки h = hi,h2,... при т = ch были найдены приближенные решения задачи с известным точным решением и вычислены величины L(u) и L(v). Одновременно были вычислены погрешности C(u) и C(v), где

C(u) = max max |u(ai,tj) — uhT(aj,tj)|. (33)

Для вычисления интегралов в определениях L(u) и L(v) использована высокоточная составная квадратурная формула. В предположении Eh « C ha, где Eh -любая из указанных выше норм, и что постоянная C слабо зависит от h, были найдены значения aj по формуле

aj = log (Ehi+i/EhJ/ log (hj+i/hj) .

Значения aj « 2 соответствуют ожидаемому значению a = 2.

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

ДЕ = max Eh(uhT(t),VhT(t)) — min Eh(uhT(t),VhT(t)).

Во всех случаях она имела порядок 10-16 ~ 10-12 в зависимости от шага h, что свидетельствует о постоянности функционала энергии.

Тестовая задача 1. Рассмотрим одномерное уравнение Кирхгофа на отрезке [—4,4] при T = 3 с точным решением u(x, t) = sin(x — t). Для решения задачи была использована равномерная сетка узлов с шагами h и т. На рис. 1 представлен график приближенного решения и его погрешности. В таблице 1 указана зависимость от h норм погрешности (32), (33), а также соответствующих им значений a.

Тестовая задача 2. Рассмотрим двумерное уравнение Кирхгофа в квадрате ^ = [0, 2] х [0,2] при T = 2 с точным решением u(x,y,t) = sin(2x + 2y)cos(t). Область ^ разбивалась на треугольные конечные элементы на основе равномерной ортогональной сетки из (n +1) х (n +1) узлов после деления одним и тем же способом каждого элементарного квадрата со стороной h на два треугольника. Равномерная сетка узлов использовалась также по переменной t (с шагом т).

Рис. 2. Приближенное решение тестовой задачи 2 при £ = Т (слева) и его погрешность (справа) на сетке узлов с шагами И = 1/20, т = И

Табл. 2

Нормы погрешности и порядки сходимости решения теста 2 при т = И

к Ь(и) ®Ь(и) т ®Цг>) С (и) ОС С (и) СИ «С!»

1/20 1.33 е-2 — 4.68 е-2 — 1.01 е-2 — 7.59 е-2 —

1/40 3.33 е-3 1.99 1.15 е-2 2.02 2.41 е-3 2.07 2.19 е-2 1.80

1/80 8.20 е-4 2.02 2.85 е-3 2.01 5.02 е-4 2.26 5.28 е-3 2.05

1/100 5.23 е-4 2.01 1.81 е-3 2.01 3.19 е-4 2.03 3.39 е-3 1.99

На рис. 2 представлены график приближенного решения и его погрешности на грубой сетке узлов. В таблице 2 указана зависимость от к норм погрешности (32), (33), а также соответствующих им значений а.

Результаты, аналогичные приведенным выше, были получены при решении других тестовых задач с гладкими решениями. Из них следует, что предложенная схема позволяет определить как решение задачи, так и его производную по времени с погрешностью порядка О (к2 + т2) в среднеквадратической норме.

Заключение

Предложена неявная двухслойная схема МКЭ для решения нелинейного уравнения Кирхгофа. Дискретная схема определяется в терминах решения задачи и его производной по переменной времени. Показано, что она обеспечивает сохранение полной энергии на дискретном уровне. Проведено численное исследование точности предложенной схемы. Численным решением ряда тестовых задач, имеющих гладкие решения, показано, что она позволяет определить как решение задачи, так и его производную по времени с погрешностью порядка О (к2 + т2) в средне-квадратической норме.

Благодарности. Работа выполнена за счет средств Программы стратегического академического лидерства Казанского (Приволжского) федерального университета ("ПРИОРИТЕТ-2030").

Литература

1. Kirchhoff G. Vorlesungen über mathematische Physik. Bd. 1: Mechanik. Leipzig: B.G. Teubner, 1876. 466 S.

2. Бернштейн С.Н. Новые приложения почти независимых величин // Изв. АН СССР. Сер. матем. 1940. Т. 4, № 2. С. 137-150.

3. Pohozaev S.I. On a class of quasilinear hyperbolic equations // Math. USSR-Sb. 1975. V. 25, No 1. P. 145-158. URL: https://doi.org/10.1070/SM1975v025n01ABEH002203.

4. Похожаев С.И. Об одном квазилинейном гиперболическом уравнении Кирхгофа // Дифференц. уравнения. 1985. Т. 21, № 1. С. 101-108.

5. Lions J.L. On some questions in boundary value problems of mathematical physics. Ser.: North-Holland Mathematics Studies. Vol. 30: Contemporary developments in continuum mechanics and partial differential equations. De La Penha G.M., Medeiros L.A.J. (Eds.). Amsterdam: North-Holland Publ. Co., 1978. P. 284-346. URL: https://doi.org/10.1016/S0304-0208(08)70870-3.

6. Arosio A., Panizzi S. On the well-posedness of the Kirchhoff string // Trans. Am. Math. Soc. 1996. V. 348, No 1. P. 305-330.

7. Arosio A. Averaged evolution equations. The Kirchhoff string and its treatment in scales of Banach spaces // Proc. 2nd Workshop on Functional-Analytic Methods in Complex Analysis and Applications to Partial Differential Equations. ICTP, Trieste, Jan. 25-29, 1993. Tutschke W., Mshimba A. (Eds.). River Edge, NJ: World Sci. Publ., 1995. P. 220254. URL: https://doi.org/10.1142/2927.

8. Lin X., Li F. Global existence and decay estimates for nonlinear Kirchhoff-type equation with boundary dissipation // Differ. Equations Appl. 2013. V. 5, No 2. P. 297-317. URL: https://doi.org/10.7153/dea-05-18.

9. Carrier G.F. On the non-linear vibration problem of the elastic string // Q. Appl. Math. 1945. V. 3, No 2. P. 157-165.

10. Cousin A.T., Frota C.L., Lar'kin N.A., Medeiros L.A. On the abstract model of the Kirchhoff-Carrier equation // Commun. Appl. Anal. 1997. V. 1, No 3. P. 389-404.

11. Cordeiro S.M.S., Pereira D.C., Ferreira J., Raposo C.A. Global solutions and exponential decay to a Klein-Gordon equation of Kirchhoff-Carrier type with strong damping and nonlinear logarithmic source term // Partial Differ. Equations Appl. Math. 2021. V. 3. Art. 100018. URL: https://doi.org/10.1016/j.padiff.2020.100018.

12. Зайцев В.В., Никулин А.В., Никулин В.В. Нелинейный резонанс в струнном ЭМР // Вестн. СамГУ. Естественнонаучн. сер. 2005. Т. 39, № 5. C. 125-130.

13. Gudi Т. Finite element method for a nonlocal problem of Kirchhoff type // SIAM J. Numer. Anal. 2002. V. 50, No 2. P. 657-668. URL: https://doi.org/10.1137/110822931.

14. Dond A.K., Pani A.K. A priori and a posteriori estimates of conforming and mixed FEM for a Kirchhoff equation of elliptic type // Comput. Methods Appl. Math. 2017. V. 17, No 2. P. 217-236. URL: https://doi.org/10.1515/cmam-2016-0041.

15. Srivastava V., Chaudhary S., Kumar V. V.K.S., Srinivasan B. Fully discrete finite element scheme for nonlocal parabolic problem involving the Dirichlet energy // J. Appl. Math. Comput. 2017. V. 53, No 1-2. P. 413-443. URL: https://doi.org/10.1007/s12190-015-0975-6.

16. Kundu S., Chaudhary S., Pani A., Khebchareon M. Fully discrete finite element scheme for nonlocal parabolic problem involving the Dirichlet energy // Numer. Funct. Anal. Optim. 2016. V. 22, No 37. P. 719-752.

17. Chaudhary S., Srivastava V., Kumar V.V.K.S. Finite element scheme with Crank-Nicolson method for parabolic nonlocal problems involving the Dirichlet energy // Int. J. Comput. Methods. 2017. V. 14, No 5. Art. 1750053. URL: https://doi.org/10.1142/S0219876217500530.

18. Peradze J. A numerical algorithm for the nonlinear Kirchhoff string equation // Numer. Math. 2005. V. 102, No 2. P. 311-342. URL: https://doi.org/10.1007/s00211-005-0642-1.

19. Bilbao S., Smith J. Energy-conserving finite difference schemes for nonlinear strings // Acta Acust. Acust. 2005. V. 91, No 2. P. 299-311.

20. Shi D., Wu Y. Nonconforming quadrilateral finite element method for nonlinear Kirchhofftype equation with damping // Math. Methods Appl. Sci. 2020. V. 43, No 5. P. 2558-2576. URL: https://doi.org/10.1002/mma.6065.

21. Лионе Ж.Л. Некоторые методы решения нелинейных краевых задач / под ред. О.А. Олейник; пер. с фр. Л.Р. Волевич. М.: Мир, 1972. 587 с.

Поступила в редакцию 10.07.2023 Принята к публикации 15.07.2023

Даутов Рафаил Замилович, доктор физико-математических наук, профессор, Институт вычислительной математики и информационных технологий Казанский (Приволжский) федеральный университет

ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: rafail.dautov@gmail.com Иванова Мария Валерьевна, студентка 2 курса магистратуры, Институт вычислительной математики и информационных технологий

Казанский (Приволжский) федеральный университет

ул. Кремлевская, д. 18, г. Казань, 420008, Россия E-mail: mashylik. 1503@gmail. com

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. 2, pp. 115-131

ORIGINAL ARTICLE

doi: 10.26907/2541-7746.2023.2.115-131

A Conservative Finite Element Scheme for the Kirchhoff Equation R.Z. Dautov*, M.V. Ivanova**

Kazan Federal University, Kazan, 420008 Russia E-mail: *rafail.dautov@gmail.com, **mashylik.1503@gmail.com

Received July 10, 2023; Accepted July, 15, 2023 Abstract

This article presents an implicit two-layer finite element scheme for solving the Kirchhoff equation, a nonlinear nonlocal equation of hyperbolic type with the Dirichlet integral. The discrete scheme was designed considering the solution of the problem and its derivative for the time variable. It ensures total energy conservation at a discrete level. The use of the Newton method was proven to be effective for solving the scheme on the time layer despite the nonlocality of the equation. The test problems with smooth solutions showed that the scheme can define both the solution of the problem and its time derivative with an error of O(h2 + t2) in the root-mean-square norm, where t and h are the grid steps in time and space, respectively.

Keywords: Kirchhoff equation, finite element method, Petrov-Galerkin method, implicit scheme, Newton method

Acknowledgments. This study was supported by the Kazan Federal University Strategic Academic Leadership Program (PRIORITY-2030).

Figure Captions

Fig. 1. Approximate solution of test problem 1 at t = T (on the left) and its error (on the right) on the grid of nodes with step size h = 1/40, t = 3h/8.

Fig. 2. Approximate solution of test problem 2 at t = T (on the left) and its error (on the right) on the grid of nodes with step size h = 1/20, t = h.

References

1. Kirchhoff G. Vorlesungen über mathematische Physik. Bd. 1: Mechanik. Leipzig, B.G. Teubner, 1876. 466 S. (In German)

2. Bernstein S.N. Nouvelles applications des grandeurs aleatoires presqu'independantes. Izv. Akad. Nauk SSSR. Ser. Mat., 1940, vol. 4, no. 2, pp. 137-150. (In Russian and French)

3. PohoZaev S.I. On a class of quasilinear hyperbolic equations. Math. US'SR-Sb., 1975, vol. 25, no. 1, pp. 145-158. URL: https://doi.org/10.1070/SM1975v025n01ABEH002203.

4. Pokhozhaev S.I. The Kirchhoff quasilinear hyperbolic equation. Differ. Equations, 1985, vol. 21, no. 1, pp. 82-87.

5. Lions J.L. On some questions in boundary value problems of mathematical physics. Ser.: North-Holland Mathematics Studies. Vol. 30: Contemporary developments in continuum mechanics and partial differential equations. De La Penha G.M., Medeiros L.A.J. (Eds.). Amsterdam, North-Holland Publ. Co., 1978, pp. 284-346. URL: https://doi.org/10.1016/S0304-0208(08)70870-3.

6. Arosio A., Panizzi S. On the well-posedness of the Kirchhoff string. Trans. Am. Math. Soc., 1996, vol. 348, no. 1, pp. 305-330.

7. Arosio A. Averaged evolution equations. The Kirchhoff string and its treatment in scales of Banach spaces. Proc. 2nd Workshop on Functional-Analytic Methods in Complex Analysis and Applications to Partial Differential Equations. ICTP, Trieste, Jan. 25-29, 1993. Tutschke W., Mshimba A. (Eds.). River Edge, NJ, World Sci. Publ., 1995, pp. 220254. URL: https://doi.org/10.1142/2927.

8. Lin X., Li F. Global existence and decay estimates for nonlinear Kirchhoff-type equation with boundary dissipation. Differ. Equations Appl., 2013, vol. 5, no. 2, pp. 297-317. URL: https://doi.org/10.7153/dea-05-18.

9. Carrier G.F. On the non-linear vibration problem of the elastic string. Q. Appl. Math., 1945, vol. 3, no. 2, pp. 157-165.

10. Cousin A.T., Frota C.L., Lar'kin N.A., Medeiros L.A. On the abstract model of the Kirchhoff-Carrier equation. Commun. Appl. Anal., 1997, vol. 1, no. 3, pp. 389-404.

11. Cordeiro S.M.S., Pereira D.C., Ferreira J., Raposo C.A. Global solutions and exponential decay to a Klein-Gordon equation of Kirchhoff-Carrier type with strong damping and nonlinear logarithmic source term. Partial Differ. Equations Appl. Math., 2021, vol. 3. Art. 100018. URL: https://doi.org/10.1016/j.padiff.2020.100018.

12. Zaitsev V.V., Nikhulin A.V., Nikhulin V.V. Nonlinear resonance in a string resonator. Vestn. SamGU. Estestvennonauchn. Ser., 2005, vol. 39, no. 5, pp. 125-130. (In Russian)

13. Gudi T. Finite element method for a nonlocal problem of Kirchhoff type. SIAM J. Numer. Anal., 2002, vol. 50, no. 2, pp. 657-668. URL: https://doi.org/10.1137/110822931.

14. Dond A.K., Pani A.K. A priori and a posteriori estimates of conforming and mixed FEM for a Kirchhoff equation of elliptic type. Comput. Methods Appl. Math., 2017, vol. 17, no. 2, pp. 217-236. URL: https://doi.org/10.1515/cmam-2016-0041.

15. Srivastava V., Chaudhary S., Kumar V.V.K.S., Srinivasan B. Fully discrete finite element scheme for nonlocal parabolic problem involving the Dirichlet energy. J. Appl. Math. Comput., 2017, vol. 53, nos. 1-2, pp. 413-443. URL: https://doi.org/10.1007/s12190-015-0975-6.

16. Kundu S., Chaudhary S., Pani A., Khebchareon M. Fully discrete finite element scheme for nonlocal parabolic problem involving the Dirichlet energy. Numer. Funct. Anal. Optim., 2016, vol. 22, no. 37, pp. 719-752.

17. Chaudhary S., Srivastava V., Kumar V.V.K.S. Finite element scheme with Crank-Nicolson method for parabolic nonlocal problems involving the Dirichlet energy. Int. J. Comput. Methods, 2017, vol. 14, no. 5, art. 1750053. URL: https://doi.org/10.1142/S0219876217500530.

18. Peradze J. A numerical algorithm for the nonlinear Kirchhoff string equation. Numer. Math., 2005, vol. 102, no. 2, pp. 311-342. URL: https://doi.org/10.1007/s00211-005-0642-1.

19. Bilbao S., Smith J. Energy-conserving finite difference schemes for nonlinear strings. Acta Acust. Acust., 2005, vol. 91, no. 2, pp. 299-311.

20. Shi D., Wu Y. Nonconforming quadrilateral finite element method for nonlinear Kirchhofftype equation with damping. Math. Methods Appl. Sci., 2020, vol. 43, no. 5, pp. 2558-2576. URL: https://doi.org/10.1002/mma.6065.

21. Lions J.L. Nekotorye metody resheniya nelineinykh kraevykh zadach [Some Methods of Solving Nonlinear Boundary Value Problems]. Oleinik O.A. (Ed.), Volevich L.R. (Transl.). Moscow, Mir, 1972. 587 p. (In Russian)

Для цитирования: Даутов Р.З., Иванова М.В. Консервативная схема метода ко/ нечных элементов для уравнения Кирхгофа // Учен. зап. Казан. ун-та. Сер. Физ.-матем. науки. 2023. Т. 165, кн. 2. С. 115-131. URL: https//doi.org/10.26907/2541-7746.2023.2.115-131.

For citation: Dautov R.Z., Ivanova M.V. A conservative finite element scheme for the / Kirchhoff equation. Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matemati-\ cheskie Nauki, 2023, vol. 165, no. 2, pp. 115-131. URL: https//doi.org/10.26907/2541-7746.2023.2.115-131. (In Russian)