Investigation of Explicit Numerical Methods for Solving Parabolic Equations

Full Text

Введение

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

Цель данной работы состоит в исследовании явных численных алгоритмов для решения нелинейных уравнений теплопроводности.

Обзор литературы

В данной области проведено большое количество исследований. В настоящее время активно развиваются новые алгоритмы решения уравнений параболического типа. Так, рассмотрены принципы построения схем численного интегрирования по времени параболических уравнений [1]. Разработан подход к интегрированию по времени системы нестационарных уравнений динамики сжимаемого теплопроводного газа [2]. Создан метод численного решения нелинейного уравнения, описывающего диффузионный перенос энергии излучения [3]. Алгоритмы решения уравнений параболического типа являются явными с отсутствием диффузионного ограничения на шаг интегрирования по времени, что может значительно увеличить скорость получения результатов при решении практических задач [4–6]. Метод, основанный на использовании многочленов П. Л. Чебышева, позволяет уйти от ограничения шага интегрирования по времени, используя итерации внутри одного шага интегрирования [1]. Метод, основанный на гиперболизации исходного уравнения, заключается в введении в исходное уравнение второй производной по времени, умноженной на малый параметр [3]. В силу явности таких методов для них легко могут быть построены параллельные вычислительные алгоритмы.

Материалы и методы

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

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

где функция u(x, t) является искомой; u(x, 0) – начальные условия; u(0, t), u(l, t) – граничные условия.

Рассмотрим пространственно-временную сетку:

Ω_h,τ=Ω_h×Ω_τ ,

где Ω_τ ={t_n nτ, 0 ≤ n ≤ N_τ} – сетка по времени с шагом τ = T/N_τ, Ω_h = {x_i = ih, 0 ≤ i ≤ N_h}; Ω_h – пространственная сетка, зависящая от шага сетки h = L/N_h, характеризующего размер ячеек; T – время, до которого ведется расчет; L – длина расчетной области.

Явная схема для уравнения (1) представлена в следующем виде:

Коэффициенты a(x_i, t) определяются из выражения a(x_i, t) = 0,5 (k (x_i, t) + k (x_{i – 1}, t)).

Отсюда получаем явное выражение для нахождения u_i^{n + 1} на следующем шаге по времени.

Алгоритм схемы локальных итераций ЛИ-М подробно описан в работах [3; 5]. Алгоритм перехода от слоя t ⁿ к слою t ⁿ⁺¹ явно-итерационный. В цикле по  вычисляется u^l:

$u^l=\frac{1}{1+\tau b_l}\left(u^n+\tau b_l{u}^{l-1}+\tau L_h{u}^{l-1}+\tau f^{n+\mathrm{0,5}}\right).$         

Результат (2p – 1)-ой итерации принимается в качестве функции на верхнем слое u^{n + 1} = u^l. Здесь L_h – диффузионный оператор;  – степень чебышевского многочлена; λ_max – значение верхней границы спектров дискретных операторов, отвечающих аппроксимациям процессов диссипативных членов. В соответствии с теоремой Гершгорина о кругах^[2] величину λ_max можно рассчитать, вычисляя для каждого узла сетки сумму модулей коэффициентов разностных аппроксимаций диссипативных членов. b_l – итерационные параметры, составленные особым образом из упорядоченного множества нулей многочлена Чебышева [1; 7].

Метод гиперболизации [3] основан на введении в параболическое уравнение второй производной по времени с малым параметром ω:

$\rho \left(x,t\right)\frac{\partial u}{\partial t}+\omega \frac{{\partial }^{2}u}{\partial t^2}=\frac{\partial }{\partial x}\left(k\left(x,t\right)\frac{\partial u}{\partial x}\right)+f\left(x,t\right).$

Трехслойная разностная схема для этого уравнения строится следующим образом:

$\frac{u_i^{n+1}-u_i^{n-1}}{2\tau }+\omega \frac{u_i^{n+1}-2u_i^n+u_i^{n-1}}{{\tau }^2}=\frac{1}{\rho \left(x_i,t^n\right)h}\left(a\left(x_{i+1},t^n\right)\frac{u_{i+1}^n-u_i^n}{h}-a\left(x_i,t^n\right)\frac{u_i^n-u_{i-1}^n}{h}\right)+f\left(x_i,t^n\right),$

где ω – малый параметр. Отсюда получаем явное выражение для нахождения u_iⁿ⁺¹ на следующем шаге по времени.

Результаты исследования

Проведено исследование трех вышеизложенных методов для решения следующего уравнения теплопроводности с переменными коэффициентами:

$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left((e^x+t)\frac{\partial u}{\partial x}\right)+(e^x-2t{e}^{2x}-tu),0<x<\mathrm{1,}0<t\le T,$ $u\left(x\mathrm{,0}\right)=\mathrm{0,}u\left(\mathrm{0,}t\right)=t,u\left(\mathrm{1,}t\right)=t{e}^1.$           (2)

Данное уравнение (2) имеет аналитическое решение:

$u\left(x,t\right)=t{e}^x.$

Проводились расчеты на сетке с шагом h = 10^–3 м до времени T = 0,2 c.

Для явной схемы шаг по времени τ = 10^–7 с определен из условия устойчивости. Для схемы локальных итераций примем шаг интегрирования τ = 10^–3, для метода гиперболизации – τ = 10^–5 с.

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

Для оценки экспериментального порядка точности [8] по пространству были проведены расчеты на последовательности сгущающихся сеток h = 10^–3 м, h = 10^–3/2 м с сохранением числа Куранта. В результате для всех методов был получен порядок аппроксимации, равный 2, что соответствует теоретическому порядку аппроксимации каждой схемы.

Также были проведены замеры расчетного времени с заданными τ и h для каждого алгоритма (таблица).

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

 

Рисунок. Сравнение численных и аналитических решений: a) явная схема; b) метод локальных итераций; c) метод гиперболизации

Figure. Comparison of numerical and analytical solutions: a) explicit scheme; b) local iteration method; c) hyperbolization method

 

 a)

b)

c)

Источник: составлено авторами.

Sources: compiled by the authors.

 

Таблица. Анализ алгоритмов

Table. Analysis of algorithms

Метод / Method	Явный / Explicit	ЛИ / LI	Гиперболизация / Hyperbolization
Отклонение между точным и приближенным решением / Deviation between the exact and approximate solution
Ошибка / Error	5,42846e-07	0,00155567	3,70791e-06
Порядок аппроксимации по пространству / The order of approximation
k	2,01954	1,90624	2,3244
Время выполнения алгоритмов, h = 10–3 м / The execution time of algorithms, h = 10–3m
τ, с	10–7	10–3+173 итерации	10–5
Время, мс / Time, ms	65428	6654	1336

Источник: таблица составлена авторами на основе данных вычислительных экспериментов.

Sources: the table was compiled by the authors based on the data from computational experiments.

 

Обсуждение и заключение

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

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

Дополнительная информация

• Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.

• Заявленный вклад авторов: Е. Е. Пескова – разработка концепции; научное руководство; валидация результатов; написание рукописи – рецензирование и редактирование. М. С. Мустайкин – проведение исследования; разработка программного обеспечения; визуализация; написание черновика рукописи.

 

^[1] Численное решение многомерных задач газовой динамики / С. К. Годунов [и др.]. М. : Наука, 1976. 400 с.; Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М. : БИНОМ. Лаборатория знаний, 2003. 632 с.

^[2] Гантмахер Ф. Р. Теория матриц. М. : Наука, 1966. 576 с.

About the authors

Elizaveta E. Peskova

National Research Mordovia State University

Email: e.e.peskova@math.mrsu.ru
ORCID iD: 0000-0003-2618-1674
SPIN-code: 4115-3762
Scopus Author ID: 57192978349
ResearcherId: U-7971-2019

Cand.Sci. (Phys.-Math.), Associate Professor of the Chair of Applied Mathematics

Russian Federation, 68 Bolshevistskaya St., Saransk 430005

Maksim S. Mustaykin

National Research Mordovia State University

Author for correspondence.
Email: maksimmustajkin@mail.ru
ORCID iD: 0009-0000-8690-0787

Student in the Mathematics and IT Department

Russian Federation, 68 Bolshevistskaya St., Saransk 430005

References

Supplementary files