CC BY
23
УДК 532.526, 532.529
Д.И. Попов, А.М. Сагалаков
Эффективный алгоритм расчета квадратичных нелинейностей в спектральной аппроксимации для дифференциальных уравнений второго порядка
D.I. Popov, A.M. Sagalakov
The Effective Algorithm of Quadratic Nonlinearity Calculation in Spectral Approximation for the Second Order Differential Equation
Рассмотрена специальная форма спектральной аппроксимации для представления нелинейных членов в виде квадратичной формы с разреженной матрицей. Обсуждаются результаты тестирования.
Ключевые слова: спектральная аппроксимация, уравнение Бюргерса, метод Галеркина, ортогональная система.
Введение. Обсуждается спектральный метод решения начально-краевой задачи для уравнений второго порядка в случае ограниченной области. Основное внимание в работе уделяется дискретизации билинейного слагаемого. В случае неограниченной области обзор сравнительно недавно полученных систематических результатов по существованию и условиям сходимости разложений по ортогональным базисам можно найти, например, в работах [1, 2]. Причем основные результаты справедливы для решений, убывающих экспоненциально или алгебраически. В случае экспоненциального затухания решения традиционно применяется метод усечения области, когда исходная бесконечная или по-лубесконечная область заменяется надлежащей ограниченной областью. Таким образом, для определенного класса задач возможно непосредственное обобщение спектрального метода, построенного для ограниченной области.
Схема аппроксимации. Известно [3], что первая начально-краевая задача уравнений параболического типа разрешима на пространстве W21•0(QT), где QT =Iх(0,Т), I=(-1,1), а через W21•0(QT) обозначено гильбертово пространство со скалярным произведением
(и, v)= | (т + uxvx)dxdt.
Вообще говоря, можно ограничиться рассмотрением таких решений и(хД), которые являются элементами W21(I), т.е. имеют первую обобщенную производную по х, и непрерывно зависят от времени. Далее будем рассматривать базис {фк , где (рк(х)=Тк(х)-Тк+2(х) (Тк(х) - полином Чебышева,
The special form of spectral approximation is concerned to represent nonliner terms as quadratic form with sparse matrix. The testing results are discussed.
Key words: spectral approximation, Burgers equation,
Galerkin method, orthogonal system.
xel), предложенный в работе [4]. Нам понадобятся пространства функций, квадратично суммируемых с весом а>( x)=(1-x2)-1/2, WjI) и L2j(I), соответственно. Например, очевидно, что L2(I) нормально вложено в L2m(I) и т.д., поскольку полиномы {Tk (x)}„ полны и образуют ортогональную систему в L2m(I) [5]. Однако сходимость в L2j( I) не гарантирует сходимость в L2 (I), поэтому условия для пробного решения необходимо усилить. А именно, известна сходимость метода Галеркина для функций {pk}" [4], когда пробное решение является
элементом HJI)<^WlJI), где j>2, HloJI) -замыкание элементов класса С” (I) в норме W2J(I). Представим пробное решение в виде
N
uN(x, t)=Z zk (t' (x).
k=0
Тогда проекция нелинейного слагаемого запишется как
N
bh (Un ,Un ) = (uNDxUN ' ) j= Z ZkZs ' )j .
k, s =0
Примем следующее обозначение:
i
Ghk =(Pkp,Ph)j = \ 9k(x)'S(x)Ph(xJ(x)dx .
-1
Пробные и поверочные функции можно записать в виде
(pn =cos n6-cos(n + 2)6=2sin6 sin(n +1)6 , где x= cos6 , 6e[0, л].
Тогда производная равна
d' =-2sin(n+1)6coS6-2(n+1)cos(n+1)6 .
dx sin 6
Рассмотрим интеграл
Gm = 4 j sm(^+\)0 cos Q sin Q s^+1)^(cos k0-cos(k+2)0)d0 . 0 sin0
Нетрудно установить, что
GT=Mk SS" S'+2+SH1+2),
где M" =(j/2)c, c0 = 2, ct =1, iФ0, а величина i определяется комбинацией h, k, s.
Перейдем к следующему интегралу:
G"k(2) =2(s+1)J cos(s+10(cos k0-cos(k+2)0(cosh0-cos("+2)0d0 .
0
Откуда получим g;;(2)=(s+i)Mhk s+"+1+s;- "+1 -s+"+3-ss-"-1)
+(s+1) M"k (sk+"+3 +sk:'"-1 -sk+"+1 S+1).
Таким образом, нелинейное слагаемое представлено в En+1 семейством из N+1 квадратичных форм, а сложность вычислений h-й проекции, как несложно заметить, является величиной порядка O(N).
Результаты тестирования. Для тестирования аппроксимации использовался одномерный вариант уравнения Бюргерса
д и / д t+ид и / д у-уд2 и / д у2 = 0 с начально-краевыми условиями
и (-да, t)=1, и (да, t)=0; и ( у,0)=1 (-да < у < 0), и(у, 0)=0 (0 < у < да), а сравнение производилось с точным решением и результатами работы [6]. Как известно, эта задача описывает распространение фронта ударной волны. При большой вязкости фронт достаточно быстро сглаживается, а при малой сохраняет свою крутизну, характеризуемую локально резким поведением решения. Однако, несмотря на бесконечный интервал, решение экспоненциально стремится к граничным условиям, поэтому применим метод усечения области. Обозначим через H радиус усечения, тогда преобразованием координат у=Hx можно привести исходную задачу к следующему виду:
д и / д t+(1/H )ид и / д x-(v/H 2)д2 и / д x2 = 0,
R=1/v , и(-1,t)=1, и(1,t)=0; и(x,0)=1 (-1< x < 0), и (x,0)=0 (0< x <1).
Для интегрирования по времени использовалась схема Крэнка-Николсона, неявная в ее линейной части,
(St )-1 (и;-ит-1, pp )+(v/2)K +ит-1, p ]+
+b" (ит-1, ит-\ph)=0, а пробное решение бралось в виде
un (x, t)=(1-x)I2+ Е zk (t)Pk (x).
k=0
Влияние гладкости начальных данных на решение можно проиллюстрировать на примере решения для Я = 10, 1 = 0.5 при N = 16. Именно, начальные данные задавались достаточно грубо с помощью метода коллокаций в узлах Гаусса-Лобатто полинома Чебышева порядка N + 1 и применением метода Га-леркина к начальным данным. Соответствующие значения среднеквадратичного отклонения равняются: ||и1а -ие||« 0.0013 и ||м2а-ие||« 0.0008, причем коэффициент затухания погрешности порядка 1/420.
На рисунке 1 приведено распределение невязки решения при различных N видно, что погрешность распределена по всей области и уменьшается с ростом величины N. Основным источником погрешности служит то обстоятельство, что приближенному решению сложно уследить за изменением точного решения, описывающего профиль ударной волны.
Таблица 1 показывает, что уже при N = 16 метод обеспечивает приемлемую точность. Однако дальнейшее увеличение N вдвое незначительно отражается на уменьшении погрешности. В таблице 2 представлены результаты сравнения спектрального метода и традиционного метода Галеркина (см. [6]) для различных чисел Рейнольдса.
0.04-г
0.02-
-0.02-
-0.04
ие-иа
1 ЧР 'v© / ~ «
>% \
“V
.А../..
-1 -0.5 0 0.5 х 1
Рис. 1. Погрешность решения уравнения Бюргерса при Я = 10
Таблица 1
Решений уравнения Бюргерса спектральным методом Я = 10, 1 = 0.5
X Приближенное решение
N = 5 N = 7 N = 9 N = 1б N = 32
RMS 0.0183 0.0079 0.0040 0.0008 0.0007
Таблица 2
Сравнение решений уравнения Бюргерса с результатами [б]
X R = 10, t = 0.92 R = 100, t = 0.92
N = 9* N = 7 N = 1б N = 9* N = 9 N = 104
RMS 0.0083 0.0058 0.0013 0.1049 0.0974 0.000б
- решение, построенное традиционным методом Галеркина в [б].
*
Анализ среднеквадратичной погрешности на основе данных таблицы 2 показывает, что при одинаковом разрешении спектральный метод дает несколько лучшие результаты, хотя в работе [6] гладкость по времени определялась схемой Рунге-Кутты четвертого порядка. Столбец при N = 104 подтверждает сходимость метода с увеличением разрешения и способность приближенного решения уловить крутой профиль ударной волны.
В таблице 3 представлены результаты сравнения со спектральным методом работы [6], в которой в качестве
базиса рассматривались полиномы Лежандра, для различных чисел Рейнольдса; в таблице 4 - сравнения спектрального метода и метода конечных элементов Галеркина. Отметим, что способ, позволяющий избежать зависимости граничных условий от времени (по сути, превращение начально-краевой задачи в задачу Коши, что неприемлемо), заключается в своевременном увеличении Н. В наших расчетах параметр Н зависел от дискретного времени и в моменты, когда численная невязка превосходила пороговое значение, производилось увеличение радиуса усечения.
Таблица 3
Сравнение решений уравнения Бюргерса с результатами [6]
X Я II і I о 40 Я = 10, 1 = 0.80 6 О II 0, 0 II &
N = 9* N = 7 N = 9* N = 7 N = 9* N = 104
ЯМ8 0.0029 0.0011 0.0086 0.0072 0.1175 0.0005
- решение, построенное спектральным методом Галеркина с использованием полиномов Лежандра [6].
Сравнение решений уравнения Бюргерса с результатами [6]
Таблица 4
X Я = 10, 1 = 0.5 2 т О II 0, 0 II &
* * 21 1 1 г 00 II 2 00 II 21 * N п 6 * * 21 1 1 г N = 40
ЯМ8 0.0011 0.0008 0.0005 0.0002 0.0119 0.0107
- решение, построенное методом конечных элементов Галеркина с использованием квадратичных КЭ [6];
- решение, построенное методом конечных элементов Галеркина с использованием линейных КЭ [6].
На рисунках 2 и 3 проиллюстрирована эволюция ударной волны при числах Рейнольдса Я =10 и Я = 100. Проведенный сравнительный анализ показал,
что предложенная схема аппроксимации является эффективной в вычислительном плане, сходящейся и устойчивой, поэтому может применяться для решения реальных задач.
Рис. 2. Эволюция фронта ударной волны при Я = 10
Рис. 3. Эволюция фронта ударной волны при Я = 100
Библиографический список
1. Shen J. and Wang L. Some recent advances on spectral methods for unbounded domains // Commun. Comput. Phys. - 2009. - Vol. 5, No. 2-4,
2. Boyd J.P. Chebyshev and Fourier spectral methods. -2nd ed. - Dover Publication Inc., Mineola. - NY, 2001.
3. Ладыженская О.А., Солонников В.А., Уральце-ва Н.Н. Линейные и квазилинейные уравнения параболического типа. - М., 1967.
4. Shen J. Efficient spectral-Galerkin methods II. Direct solvers of second and fourth order equations by using Chebyshev polinomials // SIAM J. Sci. Comput. - 1995. -Vol. 16, No. 1.
5. Натансон И.П. Конструктивная теория функций. -М., 1949.
6. Флетчер К. Численные методы на основе метода Галеркина. - М., 1988.
