ПРИБЛИЖЕННЫЕ И ТОЧНЫЕ РЕШЕНИЯ ВЫРОЖДАЮЩЕГОСЯ НЕЛИНЕЙНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ С ПРОИЗВОЛЬНОЙ НЕЛИНЕЙНОСТЬЮ Текст научной статьи по специальности «Математика»

ИНТЕГРО-ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И ФУНКЦИОНАЛЬНЫЙ АНАЛИЗ

INTEGRO-DIFFERENTIAL EQUATIONS AND FUNCTIONAL ANALYSIS

у. о

Онлайн-доступ к журналу: http: / / mathizv.isu.ru

Серия «Математика»

2020. Т. 34. С. 18-34

УДК 517.958:519.633 ЫБС 35К65

Б01 https://doi.org/10.26516/1997-7670.2020.34.18

Приближенные и точные решения вырождающегося нелинейного уравнения

______и и ^

теплопроводности с произвольной нелинейностью*

А. Л. Казаков1, Л. Ф. Спевак2

1 Институт динамики систем и теории управления им. В. М. Матросова СО РАН, Иркутск, Российская Федерация

2 Институт машиноведения УрО РАН, Екатеринбург, Российская Федерация

Аннотация. В статье рассмотрена задача о движении тепловой волны с заданным фронтом для нелинейного параболического уравнения теплопроводности общего вида. Искомая функция зависит от двух переменных. На фронте тепловой волны коэффициент теплопроводности и функция источника обращаются в нуль, что приводит к вырождению параболического типа уравнения. Это является математической причиной появления исследуемых решений, которые описывают возмущения, распространяющиеся по нулевому фону с конечной скоростью. Подобного рода эффекты, вообще говоря, несвойственны для уравнений параболического типа. Ранее авторами для рассмотренной в настоящей работе задачи была доказана теорема существования и единственности, однако она носит локальный характер и не позволяет исследовать свойства решения за пределами малой окрестности фронта тепловой волны. Для преодоления данной проблемы в статье предложен итерационный метод

* Исследования А. Л. Казакова поддержаны РФФИ (проект № 20-07-00407); Л. Ф. Спевака - Комплексной программой УрО РАН (проект № 18-1-1-5).

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

Ключевые слова: нелинейное уравнение теплопроводности, тепловая волна, метод граничных элементов, точное решение, степенной ряд.

1. Введение

Рассмотрим нелинейное параболическое [9] уравнение теплопроводности [10] с источником (стоком)

Т = ё1у(К (Т)УТ) + Н (Т), (1.1)

где К(Т) — коэффициент теплопроводности, К(0) = 0, Н(Т) — заданная функция источника (стока), Н(0) = 0. Уравнение (1.1) в различных формах и модификациях имеет широкий круг приложений в механике сплошной среды [21], популяционной биологии [19] и т. п. Важной его особенностью является вырождение, когда коэффициент теплопроводности обращается в нуль. Это приводит к появлению у него нетривиальных решений, которые, тем не менее, могут обращаться в нуль на некотором многообразии [13]. Такого рода решения (1.1) в литературе часто именуются "тепловыми волнами" [10]. Известны они давно, большое внимание им уделялось в научной школе А. Ф. Сидорова [11; 12], однако многие вопросы, связанные с существованием, единственностью и качественными свойствами тепловых волн, остаются открытыми. Авторы проводят исследования по данной тематике с 2012 г. [3;6]. При этом был разработан научный аппарат, который включает в себя конструкции в виде специальных [12] рядов и метод граничных элементов (МГЭ) [2], а также точные решения [4; 5]. Полный перечень полученных за это время результатов является обширным, поэтому упомянем лишь наиболее близкие к теме настоящей статьи. Ранее нами были рассмотрены решения типа тепловой волны для уравнения (1.1), когда К (Т) — степенная функция [1; 18]. В настоящей статье исследован случай произвольного вида зависимости коэффициента от температуры. Указанная проблематика частично была затронута в работе [1], однако здесь она изучается значительно подробнее. В частности, предложен алгоритм

построения решений рассматриваемого типа, основанный на методе граничных элементов с использованием метода двойственной взаимности [15; 20; 22], а также исследованы точные решения, построение которых сводится к интегрированию задачи Коши (ЗК) для ОДУ.

Уравнение (1.1) в случае, когда функция K(T) является монотонной и достаточно гладкой, можно заменой u = K(T) привести к виду

ut = uAu + F (u)(Vu)2 + h(u), (1.2)

где F(u) = u0"(u)/0'(u) + 1, h(u) = H(0(u))/0'(u), K(0(u)) = u, иначе говоря, 0(u) — обратная к K(T) функция. В свою очередь, уравнение

(1.2) в одномерном случае плоской симметрии приводится к виду

ut = uuxx + F (u)uX + h(u). (1.3)

Можно видеть, что параболический тип уравнения (1.3) вырождается при u = 0. В случае, когда функция K(T) является степенной, функция F(u) = 1/а = const > 0. Для (1.3) рассматривается задача о движении тепловой волны с заданным достаточно гладкой функцией x = a(t), а(0) = 0, а'(0) > 0 фронтом. Краевое условие здесь имеет вид

u|x=o(t) = 0. (1.4)

Ранее авторами для задачи (1.3), (1.4) была доказана теорема существования и единственности решения в классе аналитических функций [1]. К сожалению, указанная теорема, вообще говоря, не позволяет построить решение за пределами некоторой малой окрестности фронта тепловой волны. Кроме того, требование аналитичности функций F(u), h(u), a(t) ограничивает общность рассмотрения. С целью преодоления указанной проблемы авторы традиционно используют подход, основанный на применении метода граничных элементов [2; 6]. Далее представлена модификация указанного подхода для решения задачи

(1.3), (1.4) в случае, когда функция F(u), в отличие от ранее опубликованных работ авторов, не равна константе.

2. Алгоритм решения МГЭ

Для численного решения задачи (1.3), (1.4) на заданном промежутке времени Ь € [0, ¿*] применим пошаговый алгоритм, обобщающий результаты, приведенные в [6]. На каждом шаге по времени = кН, где Н — размер шага, будем решать краевую задачу для уравнения Пуассона, получающуюся из (1.3), (1.4) при Ь = :

пхх = С (п, щ, пх), х € [0, Ь] , п|х=ь = 0, (2.1)

где С(и,щ,их) = [щ — F(и)(их)2 + Н(и)]/и, Ь = а(1к). Из условия (1.4) следует граничное условие для теплового потока:

я1х=ь = —а'(1к )/F (0). (2.2)

Здесь д(х) = их(х)п(х), п(х) — внешняя нормаль в граничных точках, п(0) = —1,п(Ь) = 1. Решение МГЭ задачи (2.1), (2.2) имеет вид [14]

и(£) = дги*(С, 0) + д2и*(£, Ь) — щд*(£, 0) — ((£), (2.3)

где ((£) = / С (и, иг,их) и*(£,х)йх, дх = д(0), д2 = д(Ь) = —а'^к)/F(0),

0

1х{

и1 = и(0), и*(£,х) - фундаментальное решение, д*(£,х) = и%(£,х)п(х).

Неизвестные граничные значения щ и д\ удовлетворяют граничным интегральным уравнениям, получаемым из (2.3) при £ ^ 0 и £ ^ Ь:

их — дхЬ = —2((0), их + д2Ь = 2((Ь). (2.4)

Для вычисления интегралов ((£) в правых частях уравнений (2.3), (2.4) методом двойственной взаимности используем разложение по системе радиальных базисных функций (РБФ) ¡г(х) [16]

т

С (и, иг, их) = ^ аг/г(х). (2.5)

г=1

Значения ¡г(х) зависят от расстояния между текущей точкой и заданными точками коллокации х1,х2,...,хт, лежащими на отрезке [0, Ь]: ¡г(х) = ¡г (|х — хг|). Для каждой функции ¡г(х) существует такая функция «г, что ¡г = д2иг/дх2. После подстановки разложения и двукратного интегрирования по частям получаем

т

Я(£) = 52 агЯг(£), (2.6)

г=1

где (г(£) = д(0)и* (£, 0) + дг(Ь)и*(£,Ь) — щ(0)д*(£, 0) — иг(Ь)д*(£,Ь) — иг(£), дг(х) = и'г(х)п(х). Соотношения (2.3)-(2.6) позволяют построить итерационную процедуру решения задачи (2.1), (2.2) при начальном приближении и(0) = 0. Текущая п-я итерация решения имеет вид

т

(п)(£) = д{Г]и*(£, 0)+ д2и*(£,Ь) — и(П)д*(£, 0) — £ (г(£),

г=1

(п) (п)

где и1 , д1 определяются из граничных интегральных уравнений

2 N

и™ = 2 £ о^ЯгШ - д2Ь, д^ = ^ £ а^ (&(0) + ЯгЩ) - д2

г=1 г=1

(п) _ п~Т » -

Коэффициенты определяются по предыдущей итерации из решения системы уравнений = 1,..., т)

/¿(х, ) = С(>-1) ,п(га-1),пХга-1))

¿=1

Итерационный процесс заканчивается на п-й итерации, когда — п1п 1)|/|п(1га) | < е, где е - заданное число. Тогда приближенным решением задачи (1.3), (1.4) в момент Ь = является , £) = п(п)(£),£ € (0, а(Ьк)). Таким образом, разработанный алгоритм на основе МГЭ позволяет построить непрерывное по пространственной переменной решение на каждом заданном шаге по времени.

Важной составляющей алгоритма решения является выбор точек коллокации, которые должны располагаться внутри и на границе области решения задачи. Опыт авторов показывает, что наилучшим является их равномерное распределение на отрезке [0, Ь]. В рассмотренных ранее задачах для степенной функции К(Т), т. е. постоянной функции ^(и), точки коллокации выбирались следующим образом: x¿ = (г—1)Ь/т. Точка х = Ь не использовалась в качестве точки коллокации, поскольку п(Ь) = 0, и левая часть уравнения (2.5) не определена. При этом сходимость итерационных процессов и точность решения была достаточной. При рассмотрении некоторых случаев с другими зависимостями коэффициента теплопроводности от температуры возникли проблемы со сходимостью. В связи с этим в данной работе точка х = Ь была включена в число точек коллокации, x¿ = (г — 1)Ь/(т — 1). При этом в системе (14) при ] = т правая часть уравнения была заменена на пХ'Х 1)(хт), в соответствии с уравнением (2.1). Такая модификация алгоритма позволила обеспечить сходимость итерационных процедур для рассмотренных функций ^(и), отличных от постоянной.

3. Точное решение

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

Будем искать такого решения задачи (1.3), (1.4) в виде бегущей волны, т. е. и = г>(,г),,г = х — ^ — п. Ввиду инвариантности (1.3), можно

принять ^ > 0, п = 0. После подстановки в (1.3) получим ОДУ

гг'' + Е (г)(г' )2 + до' + Н(г) = 0. (3.1)

решение которого при условии г(0) = 0 определяет решение задачи

(1.3), (1.4), являющееся тепловой волной с фронтом х = Для уравнения (3.1) имеем условия Коши

г(0) = 0, г' (0) = —^/Е (0). (3.2)

При этом первое из условий (3.2) является следствием (1.4), а второе обеспечивает совместность задачи: из-за вырождения уравнения условие на производную не может быть задано произвольно, задача будет совместна только при г'(0) = — ^/Е(0) и г'(0) = 0. Первое из условий обеспечивает ненулевой тепловой поток на фронте тепловой волны. Второе приводит к решению г = 0. Понизим порядок уравнения (3.1) с помощью замены г' = р. Тогда задача (3.1), (3.2) примет вид

грр'(г) + Е(г)р2(г) + до(г) + Н(г) = 0, р(0) = —^/Е(0). (3.3)

Задача (3.3) является более простым объектом, нежели исходная (1.3),

(1.4), но по-прежнему является не разрешенной относительно производной. Тем не менее, в случае аналитических функций Е(г) и Н(г) можно построить ее решение в виде сходящегося ряда по степеням г.

Рассмотрим теперь случай степенной функции Н(п), который является наиболее распространенным в литературе [21]. Пусть Н(п) = апв, где а, в — положительные константы. Отметим, что при нецелых в задача (1.3), (1.4) при указанном выборе функции источника уже не подпадает под действие общей теоремы из работы [1]. Уравнение (3.3) имеет вид

грр' + Е (г)р2 + до + агв = 0. (3.4)

Сделаем в (3.4) замену ад = гв, в результате (3.3) запишется как

в^рр'(ад) + С(ад)р + до + аад = 0,р|ад=0 = — ^/Е (0), (3.5)

где С (ад) = Е(ад1/в). В случае аналитичности функции С(ад), используя стандартную технику [1], можно построить решение задачи (3.5) в виде ряда по степеням ад, коэффициенты которого определяются по рекуррентной процедуре, и доказать его сходимость методом мажорант. Из проведенных рассуждений следует следующее утверждение. Утверждение 1. Пусть функция С(ад) обладает следующими свойствами: 1) С(0) = 0, 2) является аналитической в некоторой окрестности ад = 0. Тогда задача (3.5) имеет единственное аналитическое 'решение р = р*(ад) в некоторой окрестности точки ад = 0.

Далее решение задачи (3.5) подставим в правую часть уравнения г' = р*(гв), которое совместно с условием г(0) = 0 образует задачу

Коши (зК). На сей раз задача разрешена относительно производной и при в > 1 подпадает под действие классических теорем существования и единственности. Итак, установлена справедливость еще одного математического факта.

Следствие. Пусть = ,в > 1, Р(0) = 0, = Р(и1/в) -аналитическая функция. Тогда задача (3.1), (3.2) имеет единственное классическое (непрерывно-дифференцируемое) 'решение.

Доказанные утверждения не позволяют определить, какова область существования построенного решения и каковы его свойства. Между тем, для того чтобы можно было уверенно использовать решения задачи (3.1), (3.2) при верификации результатов расчетов, данные вопросы являются первостепенными. Рассмотрим их в наиболее содержательном случае, когда функция Р(и) = 1/а > 0 - константа, а Л,(и) = аив, в > 0. В этом случае задача (3.5) имеет вид

1

Визр---1—и2 + ар + а:«; = 0, р(0) = —иа. (3.6)

а-ш а

Выполним в (3.6) замену переменных р = Ар, ш = = в/^, В =

ав/^2, в результате которой задача (3.6) примет вид

ыр^-+-^-р2+р + ъи = 0, р(0) =-я/3. (3.7)

аш ар

Отсюда, в частности, следует, что решение задачи (3.6) при произвольных а, в, а, ^ > 0 может быть получено из решения задачи (3.7) с помощью линейного преобразования, т. е. задача (3.6), фактически, зависит от одного параметра ва, а не от четырех.

Будем строить решение задачи (3.7) в виде ряда по степеням ш

= рп

п!

п=0

(3.8)

ад=0

Из начального условия (3.7) следует, что р0 = —ав < 0. Остальные коэффициенты ряда (3.8) определим, последовательно дифференцируя (3.7) по ш и полагая ш = 0. Коэффициенты с номерами 1 и 2 имеют вид

1 2 / N

= ГТ^2= (1 + ат+ 2аРУ■ (3'9)

И так далее. На п-м шаге получаем следующие соотношения:

= (1 + ^-41 + п./?)' П " (ЗЛ0)

= + п>2, Ь1 = Ь2 = 1. (3.11)

1 + ^

Можно видеть, что коэффициенты (3.10) положительны.

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

Оценим радиус сходимости ряда (3.8) в частном случае а в = 1, когда выкладки упрощаются. Пусть дп = рп/п! Формулы (3.10), (3.11) тогда могут быть переписаны в виде

дп = —т^—г,п>1; 4 = ^ = £¿2 = 1, (3.12)

4 2п~1(п +1) - 'га ¿.^ к+ 1 У '

Вначале построим оценку для элементов последовательности (п. Пусть Лк < М^-1/к. Здесь М1 > 0 — константа, значение которой будет определено далее. Обоснование оценки проведем индукцией по к.

База индукции Л1 = 1 < М0/1, Ь2 = 1 < М|/2 выполняется при М\ > 2. Допустим, что оценка справедлива для всех к = 1,... ,п — 1. Тогда из (3.12) получаем неравенство

п-

йп < 2М?~2 V —-—--, п > 2. (3.13)

1 к= к(к + 1)(п — к) ~ у '

Разложив дробь в (3.13) на два слагаемых, перепишем ее в виде

(п < 2М^-2

1 ( 1 \ 2 П- 1

1--+

п + 1 V п ! п(п + 1) ^ к

к=1

Осталось определить, когда выполняется оценка

п п + 1 п + 1^ к ~ 2 у '

к=1

Можно без труда убедиться, что последовательность ¿п является монотонно убывающей при п > 5. Численные расчеты показывают, что оценка (3.14) справедлива с константой М1 = М* = 2.425744 до п = 27 включительно. Можно видеть, что (^7 = 1.2 < М*/2, и значит, (3.14) также выполняется при п > 27 с константой М1 = М*.

С другой стороны, из того, что при п > 2 справедливо неравенство ¿п > 1 (это несложно проверить), следует, что Лк > 2к-1/к, к = 1, 2,... ,п,... С учетом того, что (*п ^ 1, данную оценку вряд

ли можно существенно улучшить.

Утверждение 2. В случае а в = 1 для радиуса сходимости R ряда (3.8) справедлива двусторонняя оценка 0.824489 < R < 1.

Доказательство проводится на основе полученных ранее неравенств. С помощью формулы Адамара [7] 1/R = lim \f\qn\ получаем, что 2 <

1/R < M*. Отсюда следует искомая оценка. Утверждение доказано.

Используем полученные выше результаты для анализа свойств решения задачи (3.7). Поскольку все коэффициенты ряда (3.8), начиная с pi, являются положительными, то функция p(w) является на интервале [0, R) монотонно возрастающей вместе со своими производными любого порядка. При этом p(0) = -ав < 0, а в точках w > R сумма ряда будет равна Значит, найдется точка w = wmax > 0, в которой

p(wmax) = 0. C другой стороны, решение уравнения (3.7) может обратиться в нуль при w > 0 только в том случае, если dp/dw = Таким образом, можно высказать обоснованное предположение, что wmax = R, причем ряд (3.8) в этой точке сходится и его сумма равна нулю, а ряд из производных расходится; иначе говоря, что справедлива следующая оценка для величины vmax — максимального значения решения v = v(z) задачи (3.1), (3.2) при F (v) = 1/а, h(v) = ave, а,в,а,ц, > 0, ва = 1:

f/o.824489^ < vmax < (3.15)

V а У а

Неравенство (3.15) нельзя считать строго доказанным, однако, как будет показано далее, оно хорошо согласуется с расчетами.

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

В общем случае оценить радиус сходимости аналитического решения задачи (3.1), (3.2) затруднительно. В связи с этим актуальной задачей является построение решений задачи (1.3), (1.4) типа бегущей волны с помощью решения задачи (3.1), (3.2). Сделаем для удобства замену v(z) = v(-z) и разрешим уравнение (3.1) относительно старшей производной. С учетом (3.2) имеем ЗК

= I

дв< _ F (S) (s')2 - h (5) J , 5(0) = 0, ï(0) = (4.1)

Здесь q(z) = v'(z)n(z). Непрерывное решение задачи (4.1) на заданном отрезке z € [0,L], на котором v > 0, может быть получено итерационным методом граничных элементов [18], аналогично построению решения на каждом шаге в разделе 2. Отметим, однако, необходимость применения различных подходов к решению в зависимости от наличия источника в уравнении (1.3). В случае, когда h(u) = 0, решение задачи (4.1) является монотонным. Для построения решения задачи (1.3), (1.4)

на интервале £ € [0, £*] нужно решить задачу (4.1) на отрезке г € [0, В рассмотренном в разделе 3 случае, когда Л,(и) = аив, гладкое неотрицательное решение (4.1) существует на конечном интервале г € [0, г*], при этом г(г*) = 0 и гг'(г) ^ то при г ^ г*. Неограниченность производной в граничной точке г* не позволяет применить стандартный алгоритм решения на отрезке г € [0,г*]. Будем строить решение в два этапа. На первом задача (4.1) решается на отрезке г € [0, Ь], где Ь < г* подобрано так, что г'(Ь) < 0. Решение этой задачи может быть получено итерационным МГЭ [18]. Поскольку решение на отрезке г € [Ь, г*] в постановке (4.1) получить невозможно из-за неограниченности производной, на втором этапе поменяем местами независимую переменную и искомую функцию. Получим ЗК для новой искомой функции г(гг):

г" = ^ [Р (г;) - до' + к (у) {г')2] , г(Ь*) = Ь, д2(Ь*) = 1/дЩ, (4.2)

где V € [0, Ь*], (г) = г'(гг)п(гг), Ь* = гг(Ь); г(Ь) и ;(Ь) - значения, полученные из решения на первом этапе. Задача (4.2) также может быть решена итерационным МГЭ. В том числе будет определено значение г(0) = г*. Непрерывный вид найденной функции г(гг) позволяет найти гг(г) для г € [Ь, г*] и полное решение задачи (3.1), (3.2).

Пример 1. Сравнение решений МГЭ с точными решениями. Рассмотрим задачу (1.3), (1.4) без источника (Л,(и) = 0) при Р(и) = (Л + и)/(А — и), что соответствует К(Т) = АТ/(Т + С), А, С > 0. В этом случае при а(£) = точное решение (4.1) имеет вид

г(г) = Л [1 — ехр (—до/А)]. (4.3)

Соответствующее решение задачи (1.3), (1.4) имеет вид

и(£, ж) = Л [1 — ехр (^(ж — /Л)]. (4.4)

Решение (4.3) было использовано для оценки точности решения задачи (3.1), (3.2) итерационным МГЭ. В табл. 1 приведены максимальные погрешности численного решения на отрезке г € [0,1] в зависимости от числа точек коллокации т в (2.5). В качестве РБФ были приняты

мультиквадратичные функции ^ = <^1 + 5г2. Наилучшие по точности решения результаты были получены, когда параметр формы 5 находился в интервале [0.351, 0.751 ], где 51 = 1/(0.815й), d = (^™1 а) /т, а — расстояние от х^ до ближайшей соседней точки коллокации [17].

Сравнение с точным решением (4.4) использовалось для верификации пошагового алгоритма решения (1.3), (1.4), описанного в разделе 2. В табл. 2 представлены погрешности решения МГЭ в отдельных точках в различные моменты времени при разных шагах по времени. Здесь А = 3, ^ = 0.5. Приведенные результаты расчетов демонстрируют

Таблица 1

Погрешности решения задачи (4.1) итерационным МГЭ

Л М гп = 26 гп = 51 т = 101 т = 201

3 0.5 7.5Е-05 2.4Е-05 6.7Е-06 1.6Е-06

3 1 2.7Е-04 6.8Е-05 2.2Е-05 6.8Е-06

2 0.5 9.7Е-05 2.4Е-05 8.7Е-06 2.6Е-06

2 1 3.5Е-04 1.1Е-04 3.2Е-05 9.5Е-06

Таблица 2

Погрешности пошагового решения задачи (1.3), (1.4)

Н £ ж = 0.3 а(£) ж = 0.5 а(£) ж = 0.7 а(£) ж = 0.9 а(£)

0.1 0.5 1.9Е-04 1.1Е-04 4.80Е-05 6.5Е-06

0.1 1 4.8Е-04 ЗЛЕ-04 1.50Е-04 2.6Е-05

0.01 0.5 3.5Е-05 2.5Е-05 1.4Е-05 4.4Е-06

0.01 1 7.6Е-05 5.5Е-05 3.2Е-05 8.9Е-06

численную сходимость разработанных алгоритмов МГЭ относительно как числа точек коллокации (табл. 1), так и шага по времени (табл. 2).

Пример 2. Сравнение пошагового решения МГЭ с решением типа бегущей волны. Проведенные расчеты показали достаточную точность решений ЗК (3.1), (3.2) итерационным МГЭ. Это позволяет сделать вывод, что найденные с их помощью решения задачи (1.3), (1.4) могут быть использованы для тестирования пошагового алгоритма решения МГЭ. В качестве иллюстрации приведем сравнение пошагового решения с построенным с помощью МГЭ решением типа бегущей волны при К(Т) = А (ехр(аТ) — 1). В этом случае Р(и) = А/(А+и). Сравнение двух решений при Л,(и) = и3, А = 3,^ = 1, приведенное на рис. 1а, показывает хорошее совпадение результатов расчетов.

Пример 3. Сравнение решения МГЭ с решением в виде ряда. Полученная в утверждении 2 оценка радиуса сходимости ряда (3.5) позволяет использовать его для верификации решения МГЭ, по крайней мере, до нижней оценки. На рис. 1б в фазовой плоскости представлено сравнение решения МГЭ с отрезками ряда различных степеней. Расчеты выполнены при в = а = а = ^ = 1, в этом случае ш = V. Отметим, что решения двумя методами совпадают, пока значение V не превосходит нижнюю оценку радиуса сходимости, которая в рассматриваемом случае приближенно равна 0.82. Это свидетельствует о высокой точности решения МГЭ. Для демонстрации различия в решениях зависимость V (V') представлена на рис. 1б при V > 0.45. По графикам видно, что решение в виде отрезка ряда приближается к решению МГЭ с увеличением числа членов. Анализ приведенных ре-

зультатов позволяет предположить, что радиус весьма близок к нижней оценке. С другой стороны, мы наблюдаем еще одно подтверждение корректности применяемой методики решения с помощью МГЭ. Дополнительным показателем этого является то, что представленный на рис. 1б график имеет вертикальную касательную в точке, где V = 0, и это соответствует качественному анализу решаемой ЗК.

аб

Рис. 1. Сравнение решений задачи (1.3), (1.4) (а) и задачи Коши (б)

Пример 4. Численная оценка параметра г>тах. В заключение проверим соответствие значений г>тах, найденных при решении задачи (3.1), (3.2) методом граничных элементов, оценкам (3.15). Результаты представлены в табл. 3. Здесь а = ц = 1, в = 1/а, гчтах — расчетные значения параметра, V-ах и ^+ах — нижняя и верхняя оценки из неравенства (3.15). Представленные данные показывают хорошее соответствие результатов решения полученным оценкам. При этом с ростом а значение Утах отдаляется от нижней оценки и приближается к верхней.

Таблица 3

Значения параметра ьшах

СГ ^шах ^шах v+ 17 max СГ ^шах ^тах v+ 17 max

0.5 0.642063 0.654997 0.707107 1 0.824489 0.858849 1

0.6 0.655544 0.671285 0.736022 2 2.719128 3.034428 4

0.75 0.697324 0.718924 0.805927 3 15.132768 21.201591 27

0.8 0.716836 0.740548 0.836512 3.5 40.821145 70.037634 80.211780

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

А. Л. КАЗАКОВ, Л. Ф. СПЕВАК Заключение

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

Дальнейшие исследования в данном направлении могут быть связаны с увеличением количества искомых функций и/или независимых пространственных переменных. Также интересной задачей является получение новых классов точных решений рассмотренной в статье задачи.

Авторы признательны к.ф.-м.н. П.А. Кузнецову, при участии которого была получена формула (3.10).

Список литературы

1. Казаков А. Л. О точных решениях краевой задачи о движении тепловой волны для уравнения нелинейной теплопроводности // Сибирские электронные математические известия. 2019. Т. 16. С. 1057-1068. https://doi.org/10.33048/semi.2019.16.073

2. Казаков А. Л., Кузнецов П. А., Спевак Л. Ф. Об одной краевой задаче с вырождением для нелинейного уравнения теплопроводности в сферических координатах // Труды ИММ УрО РАН. 2014. Т. 20, № 1. С. 119-129.

3. Казаков А. Л., Лемперт А. А. Аналитическое и численное исследование одной краевой задачи нелинейной фильтрации с вырождением // Вычислительные технологии. 2012. Т. 17, № 1. С. 57-68.

4. Казаков А. Л., Орлов С. С. О некоторых точных решениях нелинейного уравнения теплопроводности // Труды ИММ УрО РАН. 2016. Т. 22, № 1. С. 112-123.

5. Казаков А. Л., Орлов С. С., Орлов С. С. Построение и исследование некоторых точных решений нелинейного уравнения теплопроводности // Сибирский математический журнал. 2018. Т. 59, № 3. С. 544-560. https://doi.org/10.17377/smzh.2018.59.306

6. Казаков А. Л., Спевак Л. Ф. Методы граничных элементов и степенных рядов в одномерных задачах нелинейной фильтрации // Известия Иркутского государственного университета. Серия Математика. 2012. Т. 5, № 2. С. 2-17.

7. Корн Г., Корн Т. Справочник по математике. М. : Наука, 1974. 832 с.

8. Косов А. А., Семенов Э. И. O точных решениях уравнения нелинейной диффузии // Сибирский математический журнал. 2019. Т. 60, № 1. С. 123-140. https://doi.org/10.33048/smzh.2019.60.111

9. Ладыженская О. А., Солонников В. А., Уральцева Н. Н. Линейные и квазилинейные уравнения параболического типа. М. : Наука, 1967. 738 с.

10. Самарский А. А., Галактионов В. А., Курдюмов С. П., Михайлов А. П. Режимы с обострением в задачах для квазилинейных параболических уравнений. М. : Наука, 1987. 480 с.

11. Сидоров А. Ф. Избранные труды. Математика. Механика. М. : Физматлит, 2001. 576 с.

12. Филимонов М. Ю. Использование метода специальных рядов для представления решений начально-краевых задач для нелинейных уравнений с частными производными // Дифференциальные уравнения. 2003. Т. 39, № 8. С. 1100-1107.

13. Antontsev S. N., Shmarev S. I. Evolution PDEs with nonstandard growth conditions. Existence, uniqueness, localization, blow-up. Paris, Atlantis Press, 2015. 409 p. https://doi.org/10.2991/978-94-6239-112-3

14. Banerjee P. K., Butterfield R. Boundary element methods in engineering science. London, McGraw-Hill Book Company, 1981. 494 p.

15. Divo E., Kassab A. J. Transient non-linear heat conduction solution by a dual reciprocity boundary element method with an effective posteriori error estimator // Comput. Mater. Contin. 2005. Vol. 2, N 4. P. 277-288. http://doi:10.3970/cmc.2005.002.277

16. Golberg M. A., Chen C. S., Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM // Eng. Anal. Boundary Elem. 1999. Vol. 23. P. 285-296.

17. Hardy R. L. Multiquadric equations of topography and other irregular surfaces // Journal of Geophys. Res. 1971. Vol. 76. P. 1905-1915.

18. Kazakov A. L., Kuznetsov P. A., Spevak L. F. Analytical and numerical construction of heat wave type solutions to the nonlinear heat equation with a source. Journal of Math. Sciences. 2019. Vol. 239, N 1. P. 111-122. https://doi.org/10.1007/s10958-019-04294-x

19. Murray J. Mathematical Biology: I. An Introduction. N. Y. : Springer, 2002. 551 p.

20. Spevak L. F., Nefedova O. A. Solving a two-dimensional nonlinear heat conduction equation with degeneration by the boundary element method with the application of the dual reciprocity method // AIP Conf. Proc. 2016. Vol. 1785. P. 040077. https://doi.org/10.1063/1.4967134

21. Vazquez J. L. The Porous Medium Equation: Mathematical Theory. Oxford : Clarendon Press, 2007. 648 p.

22. Wrobel L. C., Brebbia C. A., Nardini D. The dual reciprocity boundary element formulation for transient heat conduction. Finite elements in water resources VI. Berlin : Springer-Verlag, 1986. P. 801-811.

Александр Леонидович Казаков, доктор физико-математических наук, главный научный сотрудник, Институт динамики систем и теории управления им. В. М. Матросова СО РАН, Российская Федерация, 664033, г. Иркутск, ул. Лермонтова, 134, тел.: +7(3952)453033, email: [email protected], ORCID iD https://orcid.org/0000-0002-3047-1650.

Лев Фридрихович Спевак, кандидат технических наук, заведующий лабораторией прикладной механики, Институт машиноведения УрО РАН, Российская Федерация, 620049, г. Екатеринбург, ул. Комсомольская, 34, тел.: +7(343)3753592, email: [email protected], ORCID iD https://orcid.org/0000-0003-2957-6962.

Поступила в 'редакцию 07.08.2020

Approximate and Exact Solutions to the Singular Nonlinear Heat Equation with a Common Type of Nonlinearity

A.L.Kazakov L.F.Spevak2

1 Matrosov Institute for System Dynamics and Control Theory SB RAS, Irkutsk, Russian Federation

2 Institute of Engineering Science UB RAS, Ekaterinburg, Russian Federation

Abstract. The paper deals with the problem of the motion of a heat wave with a specified front for a general nonlinear parabolic heat equation. An unknown function depends on two variables. Along the heat wave front, the coefficient of thermal conductivity and the source function vanish, which leads to a degeneration of the parabolic type of the equation. This circumstance is the mathematical reason for the appearance of the considered solutions, which describe perturbations propagating along the zero background with a finite velocity. Such effects are generally atypical for parabolic equations. Previously, we proved the existence and uniqueness theorem for the problem considered in this paper. Still, it is local and does not allow us to study the properties of the solution beyond the small neighborhood of the heat wave front. To overcome this problem, the article proposes an iterative method for constructing an approximate solution for a given time interval, based on the boundary element approach. Since it is usually not possible to prove strict convergence theorems of approximate methods for nonlinear equations of mathematical physics with a singularity, verification of the calculation results is relevant. One of the traditional ways is to compare them with exact solutions. In this article, we obtain and study an exact solution of the required type, the construction of which is reduced to integrating the Cauchy problem for an ODE. We obtained some qualitative properties, including an interval estimation of the wave amplitude in one particular case. The performed calculations show the effectiveness of the developed computational algorithm, as well as the compliance of the results of calculations with qualitative analysis.

Keywords: nonlinear heat equation, heat wave, boundary element method, exact solution, power series.

1.

2.

3.

4.

5.

6.

7.

8.

9.

10.

11

12

13

14

15

16

17

18

References

Kazakov A.L. Solutions to a heat wave propagation boundary-value problem for a nonlinear heat equation. Sib. Electronic Math. Reports, 2019, vol. 16, pp. 10571068. (in Russian) https://doi.org/10.33048/semi.2019.16.073 Kazakov A.L., Kuznetsov P.A., Spevak L.F. On a Degenerate Boundary Value Problem for the Porous Medium Equation in Spherical Coordinates. Trudy Instituta matematiki i mehaniki UrO RAN, 2014, vol. 20, no. 1, pp. 119-129. (in Russian)

Kazakov A.L., Lempert A.A. Analytical and Numerical Investigation of a Nonlinear Filtration Boundary-Value Problem with Degeneration. Vychislitel'nye tehnologii, 2012, vol. 17, no. 1, pp. 57-68. (in Russian)

Kazakov A.L., Orlov S.S. On some exact solutions of the nonlinear heat equation. Trudy Instituta matematiki i mehaniki UrO RAN, vol. 22, no. 1, pp. 112123. (in Russian)

Kazakov A.L., Orlov S.S., Orlov S.S. Construction and Study of Exact Solutions to a Nonlinear Heat Equation. Sib. Math. J., 2018, vol. 59, no. 3, pp. 427-441. https://doi.org/10.1134%2FS0037446618030060

Kazakov A.L., Spevak L.F. Boundary Elements Method and Power Series Method for One-dimensional Nonlinear Filtration Problems. The Bulletin of Irkutsk state University. Series Mathematics, 2012, vol. 5, no. 2, pp. 2-17. (in Russian) Korn G., Korn T. Mathematical Handbook. NY, McGraw-Hill Book Company, 1968, 1130 p.

Kosov A.A., Semenov E.I. Exact solutions of the nonlinear diffusion equation. Sib. Math. J., 2019, vol. 60, no. 1, pp. 93-107. https://doi.org/10.1134/S0037446619010117

Ladyzenskaja O.A., Solonnikov V.A., Ural'ceva N.N.Linear and quasi-linear equations of parabolic type. Providence, American Mathematical Society, 1968, 648 p.

Samarskii A.A., Galaktionov V.A., Kurdyumov S.P., Mikhailov A.P. Blow-up in Quasilinear Parabolic Equations. NY, Walter de Gruyte Berlin, 1995, 534 p. Sidorov A.F. Izbrannye trudy: Matematika. Mekhanika [Selected Works: Mathematics. Mechanics]. Moscow, Fizmatlit Publ., 2001, 576 p. (in Russian) Filimonov M.Yu. Representation of solutions of initial-boundary value problems for nonlinear partial differential equations by the method of special series. Differential Equations, 2003, vol. 39, no. 8, pp. 1159-1166.

Antontsev S.N., Shmarev S.I. Evolution PDEs with nonstandard growth conditions. Existence, uniqueness, localization, blow-up. Paris, Atlantis Press, 2015, 409 p. https://doi.org/10.2991/978-94-6239-112-3

Banerjee P.K., Butterfield R. Boundary element methods in engineering science. London, McGraw-Hill Book Company, 1981, 494 p.

Divo E., Kassab A.J. Transient non-linear heat conduction solution by a dual reciprocity boundary element method with an effective posteriori error estimator. Comput. Mater. Contin., 2005, vol. 2, no. 4, pp. 277-288. http:// doi:10.3970/cmc.2005.002.277

Golberg M.A., Chen C.S., Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Eng. Anal. Boundary Elem., 1999, vol. 23, pp. 285-296.

Hardy R.L. Multiquadric equations of topography and other irregular surfaces. Journal of Geophys. Res., 1971, vol. 76, pp. 1905-1915.

Kazakov A.L., Kuznetsov P.A., Spevak L.F. Analytical and numerical construction of heat wave type solutions to the nonlinear heat equation with a source. Journal of

Math. Sciences, 2019, vol. 239, no. 1, pp. 111-122. https://doi.org/10.1007/s10958-019-04294-x

19. Murray J. Mathematical Biology: I. An Introduction. NY, Springer, 2002, 551 p.

20. Spevak L.F., Nefedova O.A. Solving a two-dimensional nonlinear heat conduction equation with degeneration by the boundary element method with the application of the dual reciprocity method. AIP Conf. Proc., 2016, vol. 1785, p. 040077. https://doi.org/10.1063/L4967134

21. Vazquez J.L. The Porous Medium Equation: Mathematical Theory. Oxford, Clarendon Press, 2007, 648 p.

22. Wrobel L.C., Brebbia C.A., Nardini D. The dual reciprocity boundary element formulation for transient heat conduction. Finite elements in water resources VI. Berlin, Springer-Verlag, 1986, pp. 801-811.

Alexander Kazakov, Doctor of Sciences (Physics and Mathematics), professor, chief scientist, Matrosov Institute for System Dynamics and Control Theory SB RAS, 134, Lermontov st., Irkutsk, 664033, Russian Federation, tel.: +7(3952)453033, email: [email protected], ORCID iD https://orcid.org/0000-0002-3047-1650

Lev Spevak, Candidate of Sciences (Engineering Sciences), Head of Laboratory of Applied Mechanics, Institute of Engineering Science UB RAS, 34, Komsomolskaya st., Ekaterinburg, 620049, Russian Federation, tel.: +7(343)3753592, email: [email protected], ORCID iD https://orcid.org/0000-0003-2957-6962.

Received 07.08.2020