CC BY
9
УДК 519.632.4+519.651 БСТ: 10.14529/mmp220401
ПРИМЕНЕНИЕ ДРОБНО-РАЦИОНАЛЬНЫХ ИНТЕРПОЛЯЦИЙ ДЛЯ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ С ОСОБЕННОСТЯМИ
Б.В. Семисалов, Новосибирский государственный университет, г. Новосибирск, Российская Федерация; Институт математики им. С.Л. Соболева СО РАН, г. Новосибирск, Российская Федерация
Статья посвящена разработке, реализации и тестированию нового метода решения сингулярно-возмущенных краевых задач для нелинейных уравнений с частными производными второго порядка в прямоугольной области. Для приближения решения в методе использованы прямые (тензорные) произведения дробно-рациональных функций, полученных из интерполяционных полиномов с узлами Чебышева, записанных в барицентрической форме, с помощью специальной замены переменной. Замена делается с целью адаптации положения узлов интерполяции к особенностям искомой функции и приводит к их сгущению в окрестности больших градиентов решения. Для аппроксимации нелинейных уравнений используется сочетание итерационного метода установления и метода коллокаций, что позволяет свести задачу на каждой итерации к решению матричного уравнения Сильвестра. Такой подход приводит к существенному снижению времени вычислений. Высокая эффективность метода продемонстрирована на примере тестовой краевой задачи в квадрате, решение которой имеет пик в центре области, обусловленный наличием у неизвестной функции полюса в комплексной плоскости.
Ключевые слова: еингулярно-возмущенная краевая задача; дробно-рациональная интерполяция; метод коллокаций; быстрая сходимость.
Введение и постановка задачи
Численное моделирование стационарных процессов, описываемых сингулярно-возмущенными уравнениями с частными производными (УЧП), представляет фундаментальную проблему современной науки. Учет априорной информации о гладкости разыскиваемого решения, а также о расположении и типах его особых точек, связан с целенаправленным применением нелокальных полиномиальных и дробно-рациональных приближений при построении вычислительных алгоритмов. В рамках этой работы под особой понимается точка, в которой производная решения некоторого порядка терпит разрыв, либо, если решение есть аналитическая функция, -неустранимая особая точка решения, как функции одной из переменных, лежащая в комплексной плоскости на малом расстоянии от области задачи. Соответственно, краевые задачи с особенностями - это те задачи, решения которых имеют особые точки. Классические теоремы о погрешности полиномиального приближения (см., например, [1,2]) и недавний прогресс в области анализа дробно-рациональных приближений функций с особыми точками [3, 4] подтверждают значительную перспективу такого подхода. Его преимущества над методами, основанными на локальных аппроксимациях решения типа конечных разностей или конечных элементов, состоят прежде всего в экспоненциальном характере убывания погрешности нелокальных приближений с ростом числа степеней свободы.
Тем не менее, для эффективного применения нелокальных приближений при решении краевых задач для нелинейных УЧП необходимо разработать математический и вычислительный аппарат для построения и быстрого решения систем линейных
уравнений с заполненными матрицами, возникающих при аппроксимации исходной задачи. Возможности существенного повышения точности и скорости работы соответствующих алгоритмов связаны с применением тензорных произведений одномерных приближений в областях прямоугольных форм и новых способов адаптации позиций узлов коллокации к особенностям искомых функций.
В настоящей работе предложен, реализован и протестирован алгоритм решения краевых задач для УЧП второго порядка с нелинейной правой частью. Для адаптации алгоритма к особенностям решения использованы дробно-рациональные интерполяционные формулы, основанные на применении специального отображения, «уносящего> особые точки решения, лежащие в комплексной плоскости, на достаточно большое расстояние от области задачи. Такое отображение приводит к сгущению узлов интерполяции в окрестности особой точки и к существенному повышению скорости сходимости метода. Другая идея, позволившая снизить время работы алгоритма, связана с применением специального подхода к записи задачи линейной алгебры, соответствующей исходной дифференциальной постановке, в виде матричного уравнения Сильвестра, см. [5,6]. В рамках данной работы этот подход использован в комбинации с дробно-рациональными интерполяционными формулами, что обеспечило выигрыш как по точности приближения, так и по скорости решения задачи линейной алгебры.
При разработке метода будем пользоваться модельной краевой задачей в квадрате для уравнения с оператором Лапласа в левой части
д2 и д2и
А« =9-2 + д-2= /■:•'••//•":■• ■:•'••//:■ е » = [-1,1]2, (1)
где и(х, у) - достаточно гладкая неизвестная функция, f (х, у, и) - нелинейная правая часть. Предполагается, что на каждой стороне квадрата О заданы значения функции и, либо ее производной по нормали к границе, либо комбинации значений функции и производной. Отметим сразу, что описанный далее метод можно применять также для решения квазилинейного уравнения, при этом к левой и правой частям такого уравнения нужно добавить оператор Лапласа, а нелинейные члены с коэффициентами при производных решения перенести в правую часть. Это позволит записать уравнение в форме (1), где функция f будет зависеть также от производных и. Будем считать, что решение (1) имеет особую точку, что приводит к возникновению большого градиента внутри области О. Конкретный пример такой задачи рассмотрен в последнем разделе статьи. Тестовые расчеты и сравнения результатов, полученных с адаптацией сетки к особенности и без адаптации, показывают высокую эффективность метода.
1. Барицентрическая интерполяционная формула с узлами Чебышева и ее модификации
Идея о представлении интерполяционного полинома в барицентрической форме восходит к работам Лагранжа и Якоби. Некоторые важные результаты в этом направлении были опубликованы в первой половине 20-го века, см. [7-9] и обзор в книге [2, гл. 5]. Однако широкое распространение в вычислительной практике такие приближения получили после работы Сэлзера (Бакег) [10]. Здесь было впервые получено барицентрическое представление для интерполяционного полинома с узлами Чебышева Xе-1 = со8(^|^-7г), ] = которое оказалось простым и чрезвычайно
эффективным с вычислительной точки зрения:
N ш- и{ХсН) / ^
Х х А / Х Х А
1 = 1 1 ' 1 = 1 1
En (и) <
где и - приближаемая функция, ujj = т, = (—1)J 1 - веса интерполя-
ционной формулы, Tn (x) - полиномы Чебышева степени N.
Для интерполяционного полинома с узлами Чебышева, как и для формулы (2), имеет место неулучшаемая оценка погрешности, которую также именуют неравенством Лебега, см. [11, гл. 3]:
к 2
En{u) = ||'u(x) - pN(u, х) || < (1 + An) ||м - P%f(u) ||, AN < 1 + - log N, (3)
n
где Л^ - константа Лебега (ее оценка дана, например, в [2, гл. 15]), PN(u) - наилучшее полиномиальное приближение функции u(x) степени меньше N. Ниже сформулированы два важных следствия (3) (см. [12, гл. 4, 6]).
• Приближение (2) сходится тогда и только тогда, когда ||u — P^(u)|| = o(1/log(N)) при N то. В силу неравенства Джексона, обобщенного на случай непериодических функций, последнее условие эквивалентно тому, что приближаемая функция лежит в классе Дини - Липшица.
• Для функций, имеющих конечный порядок гладкости u Е Cr-1[— 1,1], у которых производная u(r) лежит в классе Гельдера (или Липшица) порядка а Е (0,1], при достаточно больших N имеет место соотношение < Cr( 1 + ^ log N)N~r~a, где Cr - некоторая величина, зависящая только от r (см. также [11, гл. 3, § 1]).
Отметим также важный результат, имеющий место для бесконечно гладких функций u Е C[—1,1]. Если u(x) допускает аналитическое продолжение в некоторую область комплексной плоскости C, содержащую эллипс Бернштейна Ep с фокусами в точках ±1 и суммой полуосей р > 1, то верно соотношение
^ 4Mu
■^N{U) < -- ,
Р — 1
где Mu - максимальное значение модуля аналитического продолжения u на контуре Ep, см. [2, гл. 8].
Основное преимущество барицентрического представления для интерполяционного полинома состоит в высокой скорости вычислений: для расчета приближенного значения функции u(x) в любой точке x Е [—1,1], x = xch требуется порядка O(N) операций. Кроме того, такие формулы обладают повышенной устойчивостью к ошибкам округления, см. [13], что позволяет строить на ЭВМ полиномиальные приближения степени 106 с точностью близкой к погрешности машинной арифметики, см. [2, гл. 5]. Алгоритмы, использующие барицентрические формулы, и их преимущества подробно обсуждаются в [4,14,15].
В общем случае, когда узлы и веса не связаны специальным образом, нетрудно видеть, что (2) представляет дробно-рациональное приближение. Положим сначала, что приближаемая функция является аналитической, и ее продолжение в комплексную плоскость имеет неустранимую особую точку, лежащую на малом расстоянии от отрезка [—1,1]. Следуя [16], зафиксируем значения весов Wj и будем адаптировать положения узлов xch к положению особой точки функции. Для этого будем использовать отображение д : [—1,1] ^ [—1,1], обладающее рядом свойств, описанных в следующей теореме из [16].
Теорема 1. Пусть D1 и D2 являются подобластями C, содержащими отрезок I = [—1,1], д - конформное отображение D1 ^ V2, такое что g(I) = I, и функция u : D2 ^ C такова, что композиция u о д : D1 ^ C является аналитической внутри эллипса Бернштейна £р, где р - сумма полуосей эллипса. Рассмотрим интерполяцию
1 N N
rN(u,x) = ^ Y1 Мх) = Х_3Х.' ®N(X) = Y1 (4)
N ( ) j=1 j j=1
где х- = д(хск). Тогда для любого х Е I при N ^ <х имеет место следующее асимптотическое представление погрешности:
|и(х) — Тм(и,х)| = 0(р -м).
Идея построения отображения д состоит в том, чтобы <унести> особую точку приближаемой функции в комплексной плоскости на достаточно большое расстояние от отрезка [—1,1], что приводит к существенному увеличению значения р и, как следствие, скорости сходимости приближения.
В этой работе для решения краевых задач для УЧП используются приближения (2) и (4) и метод коллокаций. Однако для учета типа краевых условий необходимо предварительно модифицировать эти приближения. Смысл таких модификаций состоит в том, чтобы после аппроксимации УЧП получить задачу линейной алгебры с низкой обусловленностью (см. замечание 1). Отметим, что непосредственная подстановка (2) или (4) в дифференциальное уравнение с последующим применением метода коллокаций ведет к системам с вырожденными матрицами, спектр которых состоит из одного нулевого значения, которому отвечают лишь два собственных вектора, что делает такой подход абсолютно непригодным для вычислений. Обоснование этого факта для интерполяционных полиномов с узлами Чебышева с полезными комментариями содержится в [11, стр. 277, 278].
Предлагаемая модификация выглядит следующим образом:
1 N
rN(u,x) = V\;i.r;-/?(.rh/(.r?; + (-i(x)u(—l) + <i(a:)w(l), (5)
dN (x) • i .7=1
где функции С(х), С±1(х) задаются так, чтобы приближение гN(и,х) автоматически
**! ~>3
удовлетворяло однородным граничным условиям того же типа, что и условия краевой задачи.
1. Для условий Дирихле u(±1) = 0
1 — x
2
1 — x2
2. Для условий Неймана и'(±1) = 0
=1
3. Для смешанных граничных условий аи(1)+ ви'(1) = 0, 7и(—1) + ви'(—1) = 0, в = 0, в = 0 функция £ *(х) определяется точно также, как в предыдущем пункте,
(x ± 1)2 ^ х7[к±1(х ^ 1) ± x — 2] — 2x ± 3 + к±1(1 ^ x)
UJ-i
С±1И ^N(x)^ (1TX3)2
a y
где Kl = -, = -.
4. Для случая, когда на одной границе заданы условия Неймана, а на другой - условия Дирихле (u(±1) = 0, u'(^1) = 0), возможны 2 варианта: когда реализуется верхний из знаков <±> и и когда реализуется нижний. Для этих вариантов соответственно
имеем:
Заметим, что, в зависимости от свойств искомого решения, существует множество подходов к построению отображения g(x), см., например, [15-17]. Однако, в этой работе применяется формула
г • тГЛ 1-1 1 - 8 -1-1 1 + 8 \ x - 1 . п -1 1 - 8 1
спх) = о + е smn < smn--h smn ---Ь smn -> , о
IV - - J 2 e j
где выражение <sinh-1> обозначает обратную функцию к гиперболическому синусу. Такой выбор функции g обусловлен возможностью явного указания координат (8, е) Е C особой точки решения. Действительно, при е ^ 1 несложно видеть, что g по построению отображает точку (8, е) в точку, лежащую значительно дальше от отрезка [-1,1]. Это приводит к существенному росту значения р в теореме 1. Кроме того, формула (6) допускает естественное обобщение на случай произвольного количества особых точек, см. [18]. Больше деталей о построении g(x) можно найти в [15]. Прежде чем приступать к разработке алгоритма решения краевых задач, отметим, что при наличии разрыва у некоторой производной приближаемой функции в точке Xd Е [-1,1], в формуле (6) следует задать 8 = Xd и установить значение е порядка 10-10. Тогда, как показывают численные эксперименты, погрешность дробно-рационального приближения (5) с ростом N демонстрирует экспоненциальный характер убывания.
2. Приближение решения и его производных
Для аппроксимации решения краевой задачи в прямоугольной области будем использовать тензорные произведения формул (5). После подстановки таких приближений в дифференциальное уравнение краевой задачи и использования метода колло-каций с узлами Xj приходим к матричным аппроксимациям операторов дифференцирования. Для выражения элементов соответствующих матриц необходимо вычислить первые и вторые производные от f(u,x) и перейти в полученных выражениях к пределу x ^ Xi .В итоге, используя правило Лопиталя и формулы из [14], приходим к выражениям:
N +1 N +1
г /(u,xi) = ^^ D(1) (xi,xj-)uj, г "(u,xi) = ^^ D(2) (xi,xj-)uj, (7)
j=o j=0
где использованы обозначения x0 = - 1, xN+1 = 1 ,
D(1)(xi ,xj) = < при i, j = 1,..., N;
Wi (xi xj )
N W
{с:У{Хг)~ T, Ui{xi-xky „ k=1,k=i
dw(x0,Xj) = DW(XN+1,Xj) = (с;у(1)^т при j = 1.....iv;
DN (-1) DN (1)
d(1) ( xi,x0) = lim Z-1(x), D(1) (xi,xN+1 )= lim (1(x) при i = 0, ...,N.
D(2) (xi ,xj) = <
Wi(xi xj)
NN Wk
tT^-Wi(xi - xk)
k=1,k=i k=1,k=i
j
• • 1 AT 2Шк ( £ Ш1 при i,j = 1,..., N, где Yífc =----- ^ —--- +
Wi (xi - Xk ) \ l=1,l=i W (xi - Xi) Xi - Xk
N (
j V ,Qn (-1) ' j
1 лт л ± WJ Ukxk
при j = 1,..., N, где Л, = J^
D<2>( xi,x0) = lim Z-i(x), D(2)(xi,xN+1) = lim Z1'(x) при i = 0,...,N. Запишем матрицы, аппроксимирующие операторы первых и вторых производных по переменной x, в виде , . ,1Ь ,
Aa = (D(1)(xi,x¿)), Afl = (D(2)(xi,xj)).
Нижний индекс <g> в обозначениях матриц подчеркивает зависимость значений элементов матриц от отображения g. Отметим также, что в случае краевых условий Дирихле строки и столбцы этих матриц с индексами 0 и N +1 исключаются, и матрицы имеют размер N х N .В случае, когда на одной границе задается краевое условие Неймана, на другой - условие Дирихле, возможны 2 варианта. Если значение функции задается в точке x = -1, а значение производной - в точке x = 1, то из матриц исключаются столбец и строка с индексами 0. В противном случае исключаются строка и столбец с индексами N +1. В итоге для любого из указанных вариантов матрицы имеют размер (N + 1) х (N + 1). В случае краевых условий Неймана или смешанных краевых условий, строки и столбцы из матриц не исключаются, и матрицы имеют размер (N + 2) х (N + 2).
В силу свойств приближения (5) матрицы Ag и Ag аппроксимируют операторы дифференцирования только для функций, удовлетворяющих однородным краевым условиям рассмотренных типов. Однако, как показано ниже, неоднородные краевые условия могут быть реализованы за счет применения добавочных функций простого вида. В дальнейших описаниях для краткости ограничимся случаем условий Дирихле и рассмотрим в качестве примера применения предложенных аппроксимаций одномерную и двумерную краевые задачи.
Одномерный случай. Пусть функция u(x) £ C2([-1,1]) является решением краевой задачи
Uxx = f (x,u), u(-1) = a, u(1) = b. (8)
Введем выражения для векторов значений функции u(x) и ее производных в узлах коллокации: U = (u(x1),..., u(xN ))T, Ux = (ux(x1),..., ux(xN ))T и Uxx = (uxx(x1), ...,uxx(xN))T. Для реализации неоднородных краевых условий положим
u(x) ~ r(u - v, x) + v(x),
где v(x) - линейная функция, удовлетворяющая условиям v(-1) = a, v(1) = b.
Теперь, используя линейные комбинации (7), запишем приближенное представление производных в задаче (8):
Ux - Afl(U - V) + Vx, Uxx - Afl(U - V) + Vxx, (9)
где V, Vx - вектора, составленные из значений добавочной функции v(x) и ее производной в точках x = xj, j = 1,..., N соответственно, Vxx - нулевой вектор.
Для реализации неоднородных краевых условий других типов в качестве добавочной функции нужно задавать полиномы второй или третьей степени. В таком случае вектор Vxx, составленный из значений вторых производных v (x), уже не будет нулевым. Кроме того, в зависимости от конкретного типа условий, к узлам коллокации
1
должны быть добавлены точки —1 и/или 1, а вектора и, Их, Ихх, V,, V,, должны быть дополнены соответствующими компонентами.
Ниже нам потребуется спектральное разложение матрицы :
Л = Да^АД-1, (10)
где Да - матрица собственных векторов, а Да - диагональная матрица собственных значений . Обозначим эти собственные значения ¿А, ] = 1,..., N.
Замечание 1. Модификация (5) барицентрического приближения (4) обеспечивает важные с вычислительной точки зрения свойства матрицы . Мы имеем в виду отсутствие нулевых, совпадающих, или комплексных собственных чисел. Строгое математическое обоснование этого факта требует достаточно громоздких выкладок, поэтому на данный момент он был проверен только в вычислениях для нескольких видов отображения д(ж). Более строгий анализ свойств матрицы с привлечением интервальных методов для случая интерполяционных полиномов с узлами Чебыше-ва проведен в [6]. Важным обстоятельством является также медленный рост числа обусловленности матрицы Да с ростом N, что обеспечивает устойчивость алгоритма, описанного в следующем разделе, к погрешностям округления действительных чисел в памяти ЭВМ, см. приложение статьи [5].
Двумерный случай. Положим Д = [—1,1]2, и и(ж,у) € С2(Д) - решение краевой задачи Дирихле для уравнения (1) в Д. Введем в области Д сетку с узлами ■!•'•/• .'//.К где Ху = д(со8^=±тг), ук = Л, (сое Щ^тт), ] = 1,...,^, к = 1,...,А', и положим, что д, к - отображения, удовлетворяющие условиям теоремы 1. Положим , = и (ж, ), (и^, = (ж ,у&), где ^ € {ж, у, жж, уу, жу} обозначает производную. Пусть и = (м^), = ((иД^) - N х К-матрицы.
Для реализации граничных условий введем функцию и(ж, у) = и (ж, у) — ^(ж, у), где ^(ж,у), задающая неоднородные граничные условия, рассчитывается специальным образом (см. ниже). Запишем сначала тензорное произведение интерполяций (5) для приближения функции и, удовлетворяющей однородным граничным условиям:
1 км
ГмхкЫ,х,у) = —---г Х^СЛ^ССу)Мх)МуЫъ,Ук), (11)
хк(ж,у) к=1 .=1
где хК (ж, у) = (ж)®к (у). Затем используем приближение и(ж,у) ~
гм хк (и, ж, у) + а(ж,у).
Функцию ^(ж, у) выразим в виде суммы V = Vх + vy, где Vх(ж, у) = ах(у)ж + в,(у), vy (ж, у) = ау (ж)у + ву (ж),
vx|y=±l = vy |х=±1 = 0. (12)
Множества значений {ах,вх}, {ау,ву} вычисляются исходя из условий на границах ж = ±1, у = ±1 соответственно, поэтому, в силу (12), v(ж, у) удовлетворяет заданным условиям на всей границе области Д.
Подставляя в (11) сначала у = у^ с последующим дифференцированием по ж и переходом к пределу ж ^ ж,, а затем подставляя ж = ж, с последующим дифференцированием по у и переходом к пределу у ^ ук, = 1,...^, к = 1,...,К, проводя те же выкладки, что и в одномерном случае (см. также [5]), получаем матрицы для аппроксимации производных и (ж, у) по переменным ж и у соответственно. Таким образом, для аппроксимации производных по ж имеем N х N-матрицы и ; а для приближения производных по у - К х К-матрицы В, и В:
и, « Л (и — V) + V,, иу « (и — V )ВТ + к, « А (и — V), иуу ъ{и- У)В1, иху « ла([/ - + Уху.
(13)
Здесь матрицы V и V^ содержат значения v(x,y) и (x, y) в узлах сетки (xj ,yk), матрицы VM вычисляются следующим образом:
Vx « VX + Ag Vy, Vy « VхBT + Vyy, Vxy « VxXBT + Ag Vyy, (14)
где VX,y - матрицы, содержащие значения (xj ,yk).
Заметим, что такой подход к аппроксимации производных подобен предложенному в работе [19], где он назван differentiation matrix technique. Однако использование модифицированных дробно-рациональных интерполяций (5) определяет новые важные свойства и перспективы применения этого подхода. Среди прочих отметим возможности естественного обобщения алгоритма этой работы для решения краевых задач произвольной размерности в областях канонических форм с различными типами краевых условий (по поводу задач Дирихле произвольной размерности см. [6]).
3. Алгоритм решения краевой задачи
Решение краевой задачи для нелинейного уравнения (1) будем искать методом установления. Для реализации метода введем новую переменную t, играющую роль времени, и нестационарный оператор Bt, который далее будем называть регуляризацией. Метод установления состоит в том, чтобы от уравнения (1) перейти к уравнению для функции u(t,x,y):
Btu = △u — f (x,y,u) (15)
с теми же краевыми условиями и некоторыми начальными данными u(0,x,y) = u0(x,y), а затем искать решение (1) как предел решений (15) при t ^ œ. Сходимость такого процесса неразрывно связана с устойчивостью по Ляпунову решения (1). В большинстве нелинейных прикладных задач строгое доказательство устойчивости провести затруднительно. Один подход к анализу устойчивости на основе вывода априорных оценок на норму решения нестационарной задачи обсуждаются в [20].
Рассмотрим простую регуляризацию Bt = -щ и введем сетку по временной переменной t с постоянным шагом т и узлами ts = ts, s = 1, 2,... . Пусть u и u - решения, полученные на итерациях с номерами s и s — 1 соответственно. Аппроксимируя производную по времени разностным отношением Btu ~ ^т-, приходим к уравнению для вычисления значений и по значениям u:
u — т Ди = u — Tf (x,y,u). (16)
Критерий остановки итерационного процесса можно записать в виде
||Btu|| = [[△u — f(x,y,u)||< £r, (17)
где £R > 0 - малое число, которое далее будем называть невязкой установления. Таким образом, стартуя с начальных данных u(0, x, y) = u0(x, y) и используя формулу (16), необходимо переходить с предыдущей на следующую итерацию по времени до тех пор, пока решение не установится, то есть пока не выполнится неравенство (17). В результате решение краевой задачи для нелинейного уравнения (1) будет найдено с невязкой, значения которой не превосходят £r.
Для реализации этой схемы будем использовать приближения, описанные в предыдущих разделах статьи, и метод коллокаций с узлами (xj,), j = 1,...,N, k = 1,...,K. При этом для аппроксимации производных в (16) применим формулы (13), (14). Перенося производные от добавочной функции v(x, y) в правую часть, а затем для краткости опуская их, приходим к системе линейных уравнений на текущей итерации метода установления:
(U — т (Ag U + иВТ )) = U — tF (U ) = H (U ), (18)
где и, и - матрицы, содержащие значения и(ж,у) в узлах коллокации на текущей и предыдущей итерациях; матрица Д(и) содержит значения функции f (ж, у, и) в тех же узлах.
Умножая (18) на матрицу спектрального разложения Я-1 слева и на матрицу Д-1 справа, получаем
и - трАи + ) = С(и), (19)
где и = Д-1иД-1, С(и) = Д-1Н(и)Д-1. Решение (19) выражается просто: элементы матрицы и имеют вид
г?, =-^--, (20)
' 1-Г№ + 4)
где д4к - элементы матрицы С(и), ¿Д, ¿В - собственные значения матриц , В^, = 1,..., N, к = 1,..., К. При этом должны быть выполнены следующие условия:
Т =1/(4 + ¿В), Ч?,к.
После вычисления элементов и значения приближенного решения (16) в узлах коллокации могут быть выражены по формуле
и = Да и Дв . (21)
Таким образом, алгоритм решения краевой задачи для уравнения (1) можно сформулировать в следующем виде:
Блок 1. Инициализация
1. Задаем шаг метода установления т и значения невязки
2. Задаем количества узлов коллокации вдоль осей ж и у (Ж и К узлов соответственно); используя имеющуюся информацию об особенностях решения, определяем отображения д(ж) и к(у); рассчитываем узлы (ж4-,), ] = 1,...,Ж, к = 1,...,К как образы узлов Чебышева под действием д и к соответственно; составляем матрицы, аппроксимирующие производные, и рассчитываем их спектральные разложения.
3. Задаем начальные значения неизвестной функции (массив ио = (мо(ж^, у к ))м хк) и всех ее производных, присутствующих в правой части уравнения. В примерах, рассмотренных ниже, начальные значения заданы нулевыми.
Блок 2. Решение задачи линейной алгебры на каждой итерации метода установления
1. Используя значения решения с предыдущей итерации (или начальные данные), рассчитываем элементы матрицы Д(и).
2. В соответствии с формулами (18), (19) рассчитываем элементы матриц Н(и) и с(и).
3. Находим элементы и и и в соответствии с формулами (20), (21) и рассчитываем приближенные значения производных с применением (13); переходим к Блоку 3.
Блок 3. Итерации метода установления
Если
М < гя, (22)
т
то алгоритм останавливается, и в качестве выходных данных представляются приближенные значения решения в узлах коллокации и. Здесь норма ||и|| обозначает максимальный по модулю элемент матрицы и.
Иначе задаем и = и и возвращаемся к блоку 2.
Теорема 2. Пусть задана сетка из N х К узлов коллокации. Предложенный алгоритм требует выполнения ^ + К)] + N + К3) операций. Для его реализации необходимо О^К + К2 + N) байт памяти, где символ «0> обозначает, величину, определенную с точностью до постоянного множителя, ^ - число итераций метода установления, необходимое для реализации неравенства (22).
Доказательство. Ниже указаны шаги алгоритма, требующие наибольшего количества операций.
- Спектральные разложения матриц, аппроксимирующих вторые производные. Применение ^Д-алгоритма требует порядка О^3 + К3) операций.
- Вычисление матрицы С(и) по формуле С(и) = Д-1^(и)Д—1 и вычисление решения и по формуле (21) на каждом шаге метода установления сводится к произведению N х N - и N х К-матриц, а также N х К - и К х К-матриц. Следовательно, в сумме требуется порядка О^К^ + К)) операций.
- Вычисление производных решения, стоящих в правой части (если такие имеются), а также производных добавочной функции v(ж,y) (см. (13), (14)). Как и в предыдущем пункте, эти расчеты сводятся к произведению матриц, что приводит к аналогичному вкладу в общее число операций.
Суммируя количество операций шагов, описанных выше, приходим к числу, указанному в условии теоремы.
Что касается объема памяти, для реализации алгоритма требуется хранить следующую информацию: 1) матрицы аппроксимирующие производные и матрицы их спектральных разложений размеров N х N и К х К; 2) N х К-матрицы, содержащие значения решения, добавочных функций и их производных; 3) значения координат узлов, параметров установления и параметров задачи, что дает пренебрежимо малый вклад по сравнению с предыдущими пунктами.
В сумме получаем, что объемы памяти, необходимые для реализации алгоритма, соответствуют заявленным в теореме. □
4. Численные эксперименты
Рассмотрим в качестве примера краевую задачу Дирихле для уравнения вида (1) в квадрате Д = [—1,1]2:
А" = 7^{(а'2 + х2)2(3у2 - 02) + (02 + у2)2(3х2 - а2)}, (ав )4
г,(±1,у)=(1 + а-2)-1(1+ф2)-1, (23)
и(а:,±1)= (1+(-)2)"1(1 + Г2)"1,
а
где а и в - положительные величины. Точное решение этой задачи имеет вид
Если а или ß достаточно малы, то функция uex(x,y) быстро меняется в окрестности точки (0, 0) (см. рис. 1, а). Отметим, что исследование существования и единственности решений задач с кубической нелинейностью в правой части вида (23) представляет отдельную задачу, которая выходит за рамки этой работы. Существование в данном случае очевидно, чего нельзя утверждать про единственность. Следовательно, вопрос, к какому решению будет сходится метод установления, является открытым. Выше было сказано, что ключевую роль в этом вопросе играет устойчивость
по Ляпунову соответствующего решения. Исследование устойчивости решения задачи (23) по линейному приближению показывает, что существенным является знак выражения при и3 в правой части. Положительный знак обеспечивает устойчивость, наличие отрицательного знака может приводить к неустойчивости. Отметим также, что устойчивости лишь по линейному приближению не всегда достаточно для однозначного вывода о сходимости. В связи с этим решающую роль играет вычислительный эксперимент. В эксперименте наблюдалась сходимость метода установления с невязкой £д порядка 10-11. При этом, как показано ниже, при достаточном числе узлов N максимум отклонений предельного решения от функции иех(х,у) на узлах сетки составляет также 10-11. Этот факт служит существенным аргументом в пользу устойчивости иех(х,у).
Рассмотрим случай а = в = 0, 05 и используем разработанный алгоритм для поиска приближенного решения иар(х, у) задачи (23). В частности используем два варианта алгоритма: вариант, в котором координаты узлов коллокации задаются корнями многочленов Чебышева, и вариант со сгущением таких узлов. Для сгущения применим вдоль координат х и у отображения вида (6) с параметрами 5 = 0, £ = 0, 05. При выборе параметров учтено, что аналитическое продолжение в комплексную плоскость функции иех(х, 0) как функции переменной х имеет полюс в точке 5+%£ £ С. То же самое можно сказать про функцию иех(0,у).
Пусть число узлов вдоль координаты х совпадает с числом узлов вдоль у ^ = К), обозначим £N = тах |иех(х,-, ) — иар(х~, )|. Тогда, увеличивая число узлов
N = К = 5,10,15,..., 200, мы будем наблюдать сходимость приближенного решения к точному, см. рис. 1, б. Здесь верхние индексы и <г» обозначают, что приближенное решение получено с применением полиномов с узлами Чебышева и дробно-рациональных приближений соответственно. Из графика видно, что при использовании адаптированных узлов скорость сходимости существенно увеличивается. В результате для достижения высокой точности требуется намного меньше вычислительных ресурсов.
Рис. 1. (а) Решение задачи (23) при а = в = 0, 05; (б) зависимость log10 от числа узлов коллокации
На рис. 2 показаны распределения узлов коллокации для двух указанных вариантов алгоритма.
В заключении исследуем быстродействие алгоритма, проводя вычисления на одном ядре персональной ЭВМ AMD Ryzen 9 5950X, 3.40 Ггц, DRAM 32 Гб, 2133 МГц. Каждая строка таблицы соответствует фиксированной величине относительной погрешности численного решения задачи (23), указанной в первом столбце. В остальных столбцах показано, сколько узлов и сколько времени в миллисекундах требуется алгоритму, чтобы получить решение с относительной погрешностью не более указанной
Рис. 2. (а) Сетка из 15 х 15 узлов коллокации с координатами в нулях многочленов Чебышева; (б) сетка со сгущением узлов на основе (6) с 8 = 0 и е = 0, 05
в первом столбце. Столбцы «Узлы Чебышева> содержат результаты, полученные при использовании узлов с координатами в нулях многочленов Чебышева, столбцы «Сгущенные узлы> содержат результаты, полученные при использовании сгущения узлов и дробно-рациональных приближений.
Таблица
Затраты времени при решении задачи (23)
Относительная Узлы Чебышева Сгущенные узлы
погрешность N время (мс) N время (мс)
10"1 100 83 17 4
1СГ3 160 297 26 7
1СГ5 220 798 34 12
10"7 304 2271 46 18
1СГ9 390 4835 56 21
10"11 478 8758 68 33
Заключение
В итоге работы сделаем несколько комментариев о применимости разработанного подхода для решения прикладных задач. Во-первых, необходимо отметить работы [15,18], в которых дробно-рациональные барицентрические приближения успешно использовались при решении одномерных нелинейных задач с особенностями для уравнения Бюргерса и уравнения теплового взрыва в рамках теории Франка-Каменецкого. Была показана особая эффективность этих приближений для анализа эволюции особых точек решений в комплексной плоскости. В частности, удалось связать разрушение гладких решений и эффект blow-up с выходом особой точки в область решения задачи.
Применение дробно-рациональных приближений в 2D случае открывает значительные перспективы при решении сингулярно-возмущенных краевых задач для УЧП. Автор планирует использовать разработанный алгоритм для конструирования метода спектральных элементов и расчета течений несжимаемой вязкоупругой полимерной жидкости в экструдере, по поводу этой задачи см. [21, 22]. Известно,
что эти течения демонстрируют различные нелинейные эффекты, включая потерю устойчивости и разрушение стационарных режимов, образование вихрей и переход к турбулентным режимам, и т.п. Исследование таких эффектов представляет значительный интерес для специалистов в области гидродинамики и разработчиков новых технологий экструзии.
Работа выполнена за счет гранта Российского научного фонда (Соглашение № 20-71-00071).
Литература
1. Achieser, N.I. Theory of Approximation / N.I. Achieser. - N.Y.: Frederick Ungar, 1956.
2. Trefethen, L.N. Approximation Theory and Approximation Practice / L.N. Trefethen. - N.Y.: SIAM, 2013.
3. Stahl, H.R. Best Uniform Rational Approximation of xa on [0,1] / H.R. Stahl // Acta Mathematica. - 2003. - V. 190. - P. 241-306.
4. Nakatsukasa, Y. The AAA Algorithm for Rational Approximation / Y. Nakatsukasa, O. Sete, L.N. Trefethen // Journal on Scientific Computing. - 2018. - V. 40, № 3. - P. 1494-1522.
5. Семисалов, Б.В. Быстрый нелокальный алгоритм решения краевых задач Неймана-Дирихле с контролем погрешности /Б.В. Семисалов // Вычислительные методы и программирование. - 2016. - Т. 17, № 4. - С. 500-522.
6. Семисалов, Б.В. Об одном подходе к численному решению задач Дирихле произвольной размерности / Б.В. Семисалов // Сибирский журнал вычислительной математики. -2022. - Т. 25, № 1. - С. 77-95.
7. Riesz, M. Uber einen Satz des Herrn Serge Bernstein / M. Riesz // Acta Mathematica. -1916. - V. 40. - P. 43-47.
8. Taylor, W.J. Method of Lagrangian Curvilinear Interpolation / W.J. Taylor // Journal of Research of the National Bureau of Standards. - 1945. - V. 35. - P. 151-155.
9. Dupuy, M. Le calcul numerique des fonctions par l'interpolation barycentrique / M. Dupuy // Comptes rendus de l'Academie des Sciences. - 1948. - V. 226. - P. 158-159.
10. Salzer, H.E. Lagrangian Interpolation at the Chebyshev Points xn,v = cos(vn/n), v = O(1)n; Some Unnoted Advantages / H.E. Salzer // Computer Journal. - 1972. - V. 15. - P. 156-159.
11. Бабенко, К.И. Основы численного анализа / К.И. Бабенко. - М.: Наука, 1986.
12. Дзядык, В.К. Введение в теорию равномерного приближения функций полиномами / В.К. Дзядык. - М.: Наука, 1977.
13. Higham, N.J. The Numerical Stability of Barycentric Lagrange Interpolation / N.J. Higham // Journal of Numerical Analysis. - 2004. - V. 24, № 4. - P. 547-556.
14. Schneider, C. Some New Aspects of Rational Interpolation / C. Schneider, W. Werner // Mathematics of Computation. - 1986 - V. 47. - P. 285-299.
15. Tee, T.W. A Rational Spectral Collocation Method with Adaptively Transformed Chebyshev Grid Points / T.W. Tee, L.N. Trefethen // SIAM Journal on Scientific Computing. - 2006. -V. 28. - P. 1798-1811.
16. Baltensperger, R. Exponential Convergence of a Linear Rational Interpolant Between Transformed Chebyshev Points / R. Baltensperger, J.-P. Berrut, B. Noel // Mathematics of Computation. - 1999. - V. 68. - P. 1109-1120.
17. Семисалов, Б.В. К вопросу о приближении гладких функций с погранслойными составляющими / Б.В. Семисалов, Г.А. Кузьмин // Труды Института математики и механики УрО РАН. - 2021. - Т. 27, № 4. - С. 111-124.
18. Идимешев, С.В. Дробно-рациональная аппроксимация в начально-краевых задачах с фронтами / С.В. Идимешев // Вычислительные технологии. - 2020. - Т. 25, № 2. -С. 63-79.
19. Gottlieb, D. Theory and Application of Spectral Methods / D. Gottlieb, M.Y. Hussaini, S.A. Orszag // Spectral Methods for Partial Differential Equations. - Philadelphia: SIAM, 1984. - P. 1-54.
20. Blokhin, A.M. Numerical Method for 2D Simulation of a Silicon MESFET with a Hydrodynamical Model Based on the Maximum Entropy Principle / A.M. Blokhin, A.S. Ibragimova // SIAM Journal on Scientific Computing. - 2009. - V. 31. - P. 2015-2046.
21. Schmidt, M. Setup and Test of a Laser Doppler Velocimeter for Investigations of Flow Behaviour of Polymer Melts / M. Schmidt, E. Wassner, H. Münstedt // Mechanics of Time-Dependent Materials. - 1999. - V. 3. - P. 371-393.
22. Kosheleva, K.B. Modeling of the Three-Dimensional Flow of Polymer Melt in a Convergent Channel of Rectangular Cross-Section / K.B. Kosheleva, G.V. Pyshnograib, M.Yu. Tolstykhc // Fluid Dynamics. - 2015. - V. 50. - P. 315-321.
Борис Владимирович Семисалов, кандидат физико-математических наук, доцент, кафедра <Дифференциальные уравнения:», Новосибирский государственный университет (г. Новосибирск, Российская Федерация); старший научный сотрудник, Институт математики им. С.Л. Соболева СО РАН (г. Новосибирск, Российская Федерация), [email protected].
Поступила в редакцию 25 января 2022 г.
MSC 65N35, 41A10, 41A20 DOI: 10.14529/mmp220401
APPLICATION OF RATIONAL INTERPOLATIONS FOR SOLVING BOUNDARY VALUE PROBLEMS WITH SINGULARITIES
B. V. Semisalov, Novosibirsk State University, Novosibirsk, Russian Federation; Sobolev
Institute of Mathematics, Novosibirsk, Russian Federation, [email protected]
In this paper, we develop, implement and test a new method for solving singularly perturbed boundary value problems for non-linear partial differential equations of the second order, which are posed in a rectangular domain. In order to approximate a solution, the tensor products of rational functions are used in the method. These functions are obtained from barycentric interpolation polynomials with Chebyshev nodes by means of a special change of variable. The purpose of this change of variable is to adapt the locations of interpolation nodes to singularities of the desired function, which leads to concentration of them in the neighborhood of large gradients of the solution. To approximate the non-linear equations, a combination of iterative and collocation methods is used. This allows to pass to the matrix Sylvester equation at each iteration and to reduce considerably the run time of algorithm. High computational performance of the method is demonstrated on the example of test boundary value problem in square domain with a known solution, which has a peak in the centre of the domain. Such singular behaviour is related with the presence of pole of the desired function in complex plain.
Keywords: singularly perturbed boundary value problem; rational interpolation; collocation method; fast convergence.
References
1. Achieser N.I. Theory of Approximation. N.Y., Frederick Ungar, 1956.
2. Trefethen L.N. Approximation Theory and Approximation Practice. N.Y., SIAM, 2013.
3. Stahl H.R. Best Uniform Rational Approximation of xa on [0,1]. Acta Mathematica, 2003, vol. 190, pp. 241-306.
4. Nakatsukasa Y., Sete O., Trefethen L.N. The AAA Algorithm for Rational Approximation. Journal on Scientific Computing, 2018, vol. 40, no. 3, pp. 1494-1522.
5. Semisalov B.V. A Fast Nonlocal Algorithm for Solving Neumann-Dirichlet Boundary Value Problems with Error Control. Numerical Methods and Programming, 2016, vol. 17, no. 4, pp. 500-522. (in Russian)
6. Semisalov B.V. On an Approach to Numerical Solutions of the Dirichlet Problem of an Arbitrary Dimension. Siberian Mathematical Journal, 2022, vol. 25, no. 1, pp. 77-95.
7. Riesz M. Über einen Satz des Herrn Serge Bernstein. Acta Mathematica, 1916, vol. 40, pp. 43-47. (in German)
8. Taylor W.J. Method of Lagrangian Curvilinear Interpolation. Journal of Research of the National Bureau of Standards, 1945, vol. 35, pp. 151-155.
9. Dupuy M. Le Calcul numerique des Fonctions par l'Interpolation Barycentrique. Comptes rendus de l'Academie des Sciences, 1948, vol. 226, pp. 158-159. (in French)
10. Salzer H.E. Lagrangian Interpolation at the Chebyshev Points xn,v = cos(vn/n), v = O(1)n; Some Unnoted Advantages. Computer Journal, 1972, vol. 15, pp. 156-159.
11. Babenko K.I. Osnovy chislennogo analiza [Fundamentals of Numerical Analysis]. Moscow, Nauka, 1986. (in Russian)
12. Dzydyk V.K. Vvedenie v teoriyu ravnomernogo priblizheniya funkcij polinomami [Introduction to the Theory of Uniform Approximation of Functions by Means of Polynomials]. Moscow, Nauka, 1977. (in Russian)
13. Higham N.J. The Numerical Stability of Barycentric Lagrange Interpolation. IMA Journal of Numerical Analysis, 2004, vol. 24, no. 4, pp. 547-556.
14. Schneider C., Werner W. Some New Aspects of Rational Interpolation. Mathematics of Computation, 1986, vol. 47, pp. 285-299.
15. Tee T.W., Trefethen L.N. A Rational Spectral Collocation Method with Adaptively Transformed Chebyshev Grid Points. SIAM Journal on Scientific Computing, 2006, vol. 28, pp. 1798-1811.
16. Baltensperger R., Berrut J.-P., Noel B. Exponential Convergence of a Linear Rational Interpolant between Transformed Chebyshev Points. Mathematics of Computation, 1999, vol. 68, pp. 1109-1120.
17. Semisalov B.V., Kuzmin G.A. On the Question of Approximation of Smooth Functions with Boundary Layer Components. Proceedings of Krasovskii Institute of Mathematics and Mechanics UB RAS, 2021, vol. 27, no. 4, pp. 111-124. (in Russian)
18. Idimeshev S.V. Rational Approximation in Initial Boundary Value Problems with Fronts. Computational Technologies, 2020, vol. 25, no. 2, pp. 63-79. (in Russian)
19. Gottlieb D., Hussaini M.Y., Orszag S.A. Theory and Application of Spectral Methods. Spectral Methods for Partial Differential Equations, Philadelphia, SIAM, 1984, pp. 1-54
20. Blokhin A.M., Ibragimova A.S. Numerical Method for 2D Simulation of a Silicon MESFET with a Hydrodynamical Model Based on the Maximum Entropy Principle. SIAM Journal on Scientific Computing, 2009, vol. 31, pp. 2015-2046.
21. Schmidt M., Wassner E., Münstedt H. Setup and Test of a Laser Doppler Velocimeter for Investigations of Flow Behaviour of Polymer Melts. Mechanics of Time-Dependent Materials, 1999, vol. 3, pp. 371-393.
22. Kosheleva K.B., Pyshnograib G.V., Tolstykhc M.Yu. Modeling of the Three-Dimensional Flow of Polymer Melt in a Convergent Channel of Rectangular Cross-Section. Fluid Dynamics, 2015, vol. 50, pp. 315-321.
Received January 25, 2022
