Разработка и анализ быстрого псевдоспектрального метода решения нелинейных задач Дирихле Текст научной статьи по специальности «Математика»

УДК 519.632.4

БО!: 10.14529/ттр 180210

РАЗРАБОТКА И АНАЛИЗ БЫСТРОГО ПСЕВДОСПЕКТРАЛЬНОГО МЕТОДА РЕШЕНИЯ НЕЛИНЕЙНЫХ ЗАДАЧ ДИРИХЛЕ

Б.В. Семисалов, Институт вычислительных технологии СО РАН, г. Новосибирск, Российская Федерация; Новосибирский государственный университет, г. Новосибирск, Российская Федерация

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

Ключевые слова: чебышевекие приближения; краевая задача; нелокальный метод без насыщения; уравнение Сильвестра; метод установления.

Введение

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

Решение указанной проблемы видится во всестороннем учете свойств дифференциальной задачи при построении вычислительной модели. Автоматический учет важнейшего из таких свойств - свойства гладкости искомого решения - обеспечивается при использовании методов без насыщения (см. [1, 2]). Известно, что применение рядов Фурье и интерполяционных полиномов с узлами Чебышева д ,ля приближения функции, удовлетворяющей условию Дини - Липшица, обеспечивает асимптотику погрешности наилучших полиномиальных приближений с точностью до медленнорастущих множителей (констант Лебега). Чем выше гладкость функции, тем выше скорость сходимости приближения.

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

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

Р эздел ы 2, 3 посвящены формированию и решению задачи линейной алгебры, со-ответствующби исходной дифференциальной постановке. Здесь на основе чебытттев~ ских приближений, метода коллокаций, итерационного метода установления, набора нестационарной регуляризации и методов, описанных в [3], получены системы линеиных алгебраических уравнений (СЛАУ) в форме уравнений Сильвестра и их тензорных обобщений. Далее предложен аппарат для работы с этими уравнениями, позволивший на несколько порядков снизить количество операций по сравнению с классической схемой метода коллокаций (КСМК). Под последней здесь и далее понимается алгоритм, работающий со СЛАУ вида Ax = b, имеющей единую разреженную матрицу A большого размера.

С точки зрения развития методов без насыщения основной интерес представляет вопрос о соответствии точности численных решений порядку гладкости искомых функций. В [1] для некоторых линейных задач предложен алгоритм, основанный на методе коллокаций и чебышевских приближениях, и дано строгое обоснование отсутствия его насыщения (см. гл. 9, §5, с. 563-572). Идея аппроксимации из [1] взята за основу метода, предложенного в настоящей работе. Однако здесь речь идет о решении нелинейных задач с использованием итерационного процесса. Общих подходов к обоснованию сходимости и ненасыщаемости подобных схем до сих пор предложено не было. В разделе 4 при решении тестовых задач показано повышение порядка аппроксимации разработанного алгоритма при увеличении порядка гладкости искомых функций в строгом соответствии с оценками точности наилучших приближений, что является косвенным подтверждением его ненасыщаемости. Важной особенностью работы является применение единого (нелокального) приближения неизвестной функции полиномом (1D случай) либо прямым произведением полиномов (2D и 3D случаи). В связи с указанными обстоятельствами разработанный метод будем именовать далее нелокальным методом без насыщения (далее - НМБН).

1. О приближении гладких функций

Ниже приведены теоремы, характеризующие асимптотики погрешности приближений, использованных в методе. Результаты раздела 4 демонстрируют те же асимптотики при использовании метода для решения краевых задач. Пусть D = [a, b] Е R, f Е C(D),

b i

\\f II =max \f (x)\, \\f \\p =(/ \f (t)\d^ 1 < P <*>,

a

при p = ж полагаем \\ • \|p = \| • \\, Wrp(M, D) = {f Е Cr (D) : \\f (r)\\p < M} - множества гладких функций, Lp-нормы r-x производных которых ограничены константой M. Далее будем работать с классами функций Wp (M,D), являющимися замыканиями множеств Wrp(M, D) в норме \\ • В рамках теории приближения функций f Е Wp (M,D) традиционно вводятся следующие понятия (см., например, [1]).

1. Аппроксимирующий конечномерный компакт Kn, элементы которого используются для приближения f (x) и, как правило, представляют ряды пли полиномы, содержащие n слагаемых или мономов.

2. Оператор приближения (или просто приближение) функции f (x) непрерывное отображение Pn, осуществляющее проекцию функционального пространства на компакт (Pn : C(D) ^ Kn или Pn : Wp (M, D) ^ Kn).

3. Наилучшее приближение функции f E C(D) - элемент en(f, Kn) E Kn, на котором достигается нижняя грань ein(f, Kn) = ein(f) = inf \\f — g\\.

g&Kn

4. Константа Лебега Лп = sup \\Pn(f)\\/\\f \\ - норма оператора Pn, харак-

f&C[a,b], f=0

теризующая точность и устойчивость метода приближения:

\\f — Pn(f)\\< (1 + A»)en(f),

\\f — Pn(f + 8)\\< (1+AnK(f)+Лп6.

Здесь (f + 8) - функция, полученная возмущением f(x) в каждой точке x E D не более чем на 8\\f \\, 8 - малое число.

5. Метод без насыщения (нестрого) - это метод приближения, обладающий порядком сходимости наилучшего полиномиального приближения для любой гладкости (регулярности) аппроксимируемой функции. Строгое определение, основанное на анализе асимптотики поперечников Александрова, дано в [2]. Методы без насыщения нашли развитие и применение в работах С.Д. Алгазина, А.Л. Афендикова, В.Н. Белых, A.M. Блохина.

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

I. Для функций с ограниченной гладкостью f E W£ (M,D) в соответствии с теоремой Джексона имеет место степенная сходимость наилучшего приближения.

Теорема 1. Для любой 2п-периодической г-гладкой функции f (t) и числа n E N, n > г существует тригонометрический полином Tn(t) такой, что

\f (t) — Tn(t)\< 12Г+1 - Jf M; Г) , (1)

nr \ n J

где w(f; q) = sup \f (x+h)—f (x) \ - модуль непрерывности f первого порядка.

x&[0,2n-h], h&[0,q]

Конкретный вид Tn(t) и доказательство теоремы см. в [5, с. 204-208]. Теорема 1 была обобщена на случай приближения функций f E Wp(M,D) [1, с. 158-159].

Теорема 2. Для наилучших приближений имеет место оценка

sup e*n(f,Kn) < MCrn-r, (2)

f ewr(M,D)

Cr г

Важно отметить, что в если u(f(r); q) < Mqa, 0 < а < 1, то для дан ной f существует более точная оценка [5, с. 244]:

e*n(f,Kn) < MCln-r-a, (3)

где C} зависит только от r.

II. Для бесконечно дифференцируемых функций, имеющих особенность (полюс, устранимую особую, либо существенно-особую точку) в комплексной плоскости C, Бернштейном доказана возможность построения приближений, имеющих экспоненциальную скорость сходимости.

Теорема 3. [4, 5, с. 465-467] Пусть Ш - ограниченное множество с односвязным дополнением в C. Тогда если f (z) (z Е Ш) - аналитическая в области Ш1 Е Ш функция, имеющая особенность на границе Г = дШ, то существуют число q (0 < q < 1) и последовательность полиномов Pn(z), такие что

\f(z) — Pn(z)\ < AfCqn, Vz Е Шъ C = const, (4)

где Af = max \f (z) \. Выбирая Ш и Mi такими, что D С Mi С Ш, получаем

экспоненциальную сходимость в D со знаменателем q, величина которого зависит от выбора Ш и Ш1.

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

Теорема 4. Пусть f (x) - целая функция на отрезке [a,b], M(n) - последовательность чисел, ограничивающих значения производных f: \\f (n)\\ < M(n). Тогда имеет м,есто оценка

< MЩ - a)"21~2" = M(n)qn ш MM = о. (5)

n! n! п^ж n!

Последнее предельное равенство для целых функций выполнено в силу формулы Коши - Адамара. Отметим, что (5) соответствует приближению интерполяционным полиномом с узлами Чебышева. В случае M(n) < an, a Е R, имеет место фактори-альная сходимость.

Для реализации свойств наилучших приближений

в данной работе будем пользоваться разложением функций в ряд Фурье и в ряд с базисом, состоящим из полиномов Чебышева Tn (x)

= cos(narccos(x)). В [1] обосновано отсутствие насыщения таких приближений. Кроме этого^ их константы Лебега растут медленно с ростом n

Теорема 5. Константа Лебега приближений частичными суммами ряда Фурье

4

Фп логарифмически зависит от n, Лп = — Inn + Rn, где Rn < 3.

п

Теорема 6. Константа Лебега интерполяционного полином,а с у злам,и в нулях

2

многочлена Чебышева Tn логарифмически завис ит от n, Лп = — In n +1 — Qn, 0 <

п

On < 4.

2. Аппроксимация функций и их производных

В данном разделе описан нелокальный метод решения краевых задач Дирихле размерностей й = 1, 2, 3 для уравнений эллиптического типа второго порядка вида

й =1: пхх = /(и,х), й =2: пхх + Щу = /(и,х,у),

й =3: Пхх + Пуу + Щг = / (п,х,у, г),

где и Е С (Б) - неизвестная фун кция, / - известная правая часть (возможно, разрывная). Разрабатываемый метод несложно модифицировать для решения уравнений с переменными коэффициентами и квазилинейных уравнений (см. ниже).

Предлагаемый подход восходит к работам [3, 6] и представляет псевдоспектраль-ныи метод. В нем для приближения решения использованы интерполяционные полиномы на отрезке [ — 1,1] с узлами в нулях многочлена Чебышева х3 = со$,((2] — 1)п/2М), ] = 1,..., М, представленные в форме Лагранжа:

Р1В(п,х) = ^ Xх п(х)> С(х,х) = Т-'Ь. (7)

(х х 3)Т И(х3) 1 х 3

Функция £(х,хз) обеспечивает выполнение поставленных краевых условии. Неизвестную функцию будем приближать суммой п(х) & Рю(п,х) + у(х), где Рщ(и, ±1) = 0, ь(х) - линейная функция, удовлетворяющая неоднородным условиям. Рассмотрим далее метод аппроксимации первых и вторых производных п(х) в предположении, что имеют место условия Дирихле: п(—1) = ф п(1) = ф. Определим ь(х) = ((ф — ф)х + (ф + ф))/2. Тогда

п'(х) & Р'п(п,х) + (ф — ф)/2; п"(х) & Р"в(п,х). Положим в = 1 — х2, 83 = — х2р из = и(хз). Дифференцируя Рщ по х, получаем

Р'ш м=N е — 1У-1{ х+и' ■ (8)

3 = 1 V (х хз )2вз вз(х хз)

р" ( ^ из-1! 2в2 — М2(х — хз, л (х2 — 3ххз + 2) ' , ^

Р'0(и,х) = — у(—1)3 ---Тм(х)----— Тм(х)> из. 9)

N 3=1 ^ (х — хзГвз вз(х — хз)2 )

* Г2в2 — N 2(х — хз )2

(гх _ хх , ) 3 8 . ^ ^ в ■ ( Т _ Т ■ ) 2

х х 3) в3 в3(х х 3)

Переходя в (8), (9) к пределу при х ^ хг и используя правило Лопиталя, находим

3хг

Р'п (их)= £ (—1)г+3 8 (Хгг_ х.) П3 — ^ иг, (10)

3=1,3=7 г 3 г

Р';в(и,хг) = Е ( ~1Г'-1 2в\ + 3хг<-:т х3и,~ +38 +3хг иг. (11)

. -. . 838г(хг х3) 38г

3=1,3=г г

Введем обозначения:

аг3 = {—1УУ\,аг3 = (—1)г+3 28 + (х ^г ~*3), ^ = 1,-.,М,г= (12) 83 (хг х3) 83 8г(хг х3)

П = —3хг/(2в2г), иг = — ((N2 + 5)s~ + 3xi)/(3s4), г = 1,...,N

(13)

U

( u1 \

U2

\Un J

Uxx

Ux

( (Ux)i \ (Ux)2

\ (Ux)n

( Ux(xi) \ Ux(x2)

V Ux(xN)

( (Uxx)i \ ( Uxx(xi) \ (uxx)2 Uxx (x2)

У (Uxx)n J \ Uxx (Xn) /

и сформируем N x N матрицы

(ni ai2 ...

A

a2i n2 ••• a2N

A

( vi ai2 a2i V2

arn\

a2N

\aNi aN2 ... nN J \aNi aN2 ... vN J

Для аппроксимации производных в уравнениях (6), получим формулы

Ux ~ AU + Vx, Ux

AU,

(14)

где ^ - вектор с компонентами (ух)г = -(ф+ф)/2, = 1,..., N. Далее для построения алгоритмов воспользуемся спектральным разложением матрицы Л, аппроксимирующей вторые производные!

Л = ЯлПаЯ-1, (15)

где Я а - матрица собственных вектров Л; Б а - диагональная матрица собственных значений Л - ¿А, ] = 1,..., N.

Замечание 1. При вычислении спектра матрицы Л для N = 2,..., 200 с помощью фЯ-алгоритма, реализованного в МАТЬАВ, было установлено, что собственные чис-Л

между числами есть 6,88). Поэтому разложение (15) корректно. Удачным обстоятельством является медленный рост обусловленностей матриц Яа с ростом N, обеспечивающий устойчивость разрабатываемого алгоритма к погрешностям округления.

Замечание 2. При численной реализации алгоритма матрицы А, Л, ЯА, БА, Я—1 для разных N целесообразно занести в базу дан н ых • чтобы не тратить время на их расчет.

В двумерном случае (d = 2, D = [—1,1]2, U(x,y) Е C(D) - искомая функция) введем в области D сетку с узлами (xj ,yk), где xj = cos 22— П Ук = cos Щг-1 П j = 1,...,N, k = 1,K. Положим Ujk = U(xj,Ук), (Up)jk = Up(xj,yk), где ^ обозначает одну из производных ц Е {x, y, xx, yy, xy}. Пусть U = (Ujk), Uр = ((Uj) ~ N x K матрицы.

U(x, y)

(7): n к

P ( ) = z(x,xj)TN(x) z(y,yk)tk(y)

2D ,y) h k (x — xj )T N (xj) (y — yk)T'k (yk)

-U(xj ,yk).

(16)

Тогда u(x,y) — P2D(u,x,y) + v(x,y), P2D удовлетворяет однородным, a v(x,y) -неоднородным граничным условиям: v = vx + vy, где vx(x,y) = ax(y)x + /3x(y), v (x,y) = ay (x)y + I3y(x),

vx\y=±i = vy |x=±i = 0. (17)

Значения (ax,fix), (ay, /3y) подбираются таким образом, чтобы удовлетворить условиям на границах x = ±1 y = ±1 соответственно. При этом в силу (17) после суммирования v = vx + vy функция v удовлетворяет необходимым краевым условиям на границе D.

Проводя выкладки, аналогичные (8)—(14) (см. также [3]), получаем матрицы, аппроксимирующие производные u(x,y) A, B, A, B:

Ux — AU + Fx, Uy — UBT + Vy, Uxx — AU, Uyy — UBT, Uxy — AUBT + Vxy. (18)

Здесь матрицы значений v^(x,y) на введенной сетке V и задаются формулами

Vx = Vxx + AVy, Vy = VxBT + Vy, Vxy = VxBT + AVyy, (19)

где V/x'y - массивы, содержащие vU'y(xj ,yk)•

В трехмерном случае (d = 3, D = [—1,1]3, u(x,y,z) E C(D) - искомая функция) введем в области D сетку с узлами (xj,yk,zm), где xj = cos n yk = cos П zm = cos 2m—1 П j = 1, ■■■N, к = 1, ...,K, m = 1,..., M. Положим umk = u(xj,yk, zm), U)m = u»(xj ,yk, zm), Ц E {x, y, z, xx, yy, zz, xy, xz, yz}. Пусть U = j), Uu = ((uj)

- трехмерные массивы размера N x K x M. Для приближения u(x,y,z) используем прямое произведение:

р ( ) = N — M Z (x,xj )TN(x) Z (y,yk)TK (y) Z (z,zm)TM(z) ( )

P3D (u,x,y,z) = U m=i (x — xj )TN (x3 ) (y — yk )T— (yk ) (z — zm)TM (zm) ^ ^ ^

(20)

Тогда u(x, y, z) — P3D(u,x,y,z)+ v(x,y,z), P3D удовлетворяет однородным, a v(x,y,z)

- неоднородным граничным условиям и строится аналогично 2D случаю.

Определим слоистые трехмерные массивы A, B, C A, B, C, аппроксимирующие первые и вторые производные функции u(x,y,z) по x, y, z. Пусть k и a j k элементы массивов A и A; j и ajk - элементы матриц A, A. Слоистые масс ивы A и A содержат в слоях матрицы A и A, расположенные вдоль оси z, т. е. Ут j = ajk, Ojk = ajk- Аналогично мае сивы B, B и C, C содержат в слоях матрицы B, B и C, C,

xy

Определение 1. Боковым произведением слоистого массива A со слоями A =

_ s

(ajS)NxN, лежащими вдоль z, и массива U назовем массив а = A x U с элемен-

N

тами a^k = Y1 ajSumk-

S=1

Определение 2. Фронтальным, произведением массива B со слоями B = (bks)KxK,

— f

лежащими вдоль x, и массива U назовем, массив в = B x U с элементами j =

K

bksujms

S=1

Определение 3. Верхним произведением слоистого массива С со слоями С =

_ и

(стз]мхм, лежащими вдоль у, и массива и назовем массив 7 = В х и с элементами м

1¿к стви ^]к-

в=1

Такие произведения, основанные на размещении слоистого аппроксимирующего

и

ных слоях, позволяют записать операции дифференцирования в виде

- 5 - / - и - 5 - / - и

их и А х и, иу и В х и, их и С х и, ихх и А х и, иуу и В х и, ихх и С х и,

- s - f - s - u - f - u

UXy w A x B X U, Uxz w A x C x U, UyZ w B x C x U.

(21)

Здесь для краткости 3D массивы V^ с производными vß(xj,yk,zm) опущены. Основываясь на введенных операциях, можно записать спектральные разложения слоистых массивов, например, A = Ra x Da x R—\ оде Ra, Da, R-1 - слоистые массивы со слоями, расположенными вдоль оси z и содержащими матрицы Ra, Da, R-1- Подробное описание и обоснование корректности введенных операции дано в о .

Описанные приближения неизвестной функции и ее производных можно обобщить на случай произвольного порядка и размерности дифференциального уравне-

3. Реализация алгоритма

После аппроксимации неизвестной функции и ее производных в уравнениях (6) приходим к нелинейным соотношениям, записанным для векторов (d = 1), матриц (d = 2) или 3-мерных массиbob (d = 3). Для решения таких соотношений будем строить итерационные процессы на основе метода установления. Пусть имеется краевая задача для уравнения эллиптического типа

L(x, u)u = f (x, u), u\dD = g(x), (22)

где L = L(x, u) оператор главных коэффициентов, D = [— 1,1]d, x E D С Rd. Метод установления состоит во введении новой временной переменной t и нестационарного оператора - регуляризации Bt с последующей заменой исходной краевой задачи (22) на начально-краевую задачу с неизвестной функцией u(t, x):

Btu = L(x,u)u - f (x,u), u(t, x)\x£dD = g(x), u(0, x) = uo(x). (23)

Решение задачи (22) ищется как предел решений (23) при t ^ ж>.

Рассмотрим простую регуляризацию Bt = —, введем сетку по временной пере-

dt ^ менной ^постоянным шагом т и узлами ts = ts, s = 1, 2,... . Пусть и, u - решения на s-m и (s — шагах то времени, I - единичный оператор. После дис-

u — u

кретизации Btu w - и «выделения» в уравнении (23) оператора Лапласа А

т

(L(x, u) = L(x, u) — А + А) запишем схему

u = u + т (Au — (f + Au — Lu)), ИЛИ u = (I — т A)-1[u — т (f + (А — L)u)]. (24)

Сходимость итерационного процесса для конкретного вида Ь может быть исследована на основе принципа сжимающих отображений. Критерий остановки процесса

\\ВМ\/\\п\\ = \\Ь(х,и)и — /(х,и)\\/\\и\\< ез

(25)

дает приближенное решение задачи (22) с невязкой установления £s■

Записывая простую регуляризацию уравнений (6), где Ь = А, в соответствии с (24). получим неявные схемы метода установления:

с1 = 1: 1 = 2: 1 = 3:

(/ Т д х2)

и = и — т/(и, х),

I -

I -

д

д

дх2 ду2

и = и — т/(и, х),

(26)

д д д , .. . дХ2 + ду-2 + д^2] )и = и — Т/(и,х)■

Аппроксимируя в (26) вторые производные по формулам вида (14). (18), (21) и перенося производные от добавочных функций V в правую часть (с целью сокращения выкладок они далее опущены) получаем

(Е — тА)и = и — тГ (и) = Н (и),

(и — т (Аи + иВт)) = и — тГ (и) = Н (и),

— /

(и — т (А х и + В х и + С х и )) = и — тГ (и) = Н (и),

(27)

(28) (29)

где Е - единичная матрица, и, и - вектора, матрицы или трехмерные массивы, содержащие значения неизвестных функций в узлах интерполяции на текущем и предыдущем временных слоях, Г содержит значения правой части в узлах сетки (см. рис. 1 для случая 1 = 3).

Рис. 1. Представление уравнения (29) в виде операций с трехмерными массивами

Умножая (27) на матрицу Я—1 слева; (28) на матрицу Я—1 слева и Я—1 справа; или записывая боковое произведение (29) на слоистый массив Я—1, фронтальное произведение на Я—11 и верхнее произведение на Я—1 с учетом ассоциативности операций произведения, получаем системы с диагональными матрицами:

(е—тба)\и = с(и), и — т (бау + увв) = с(и),

и — т(лА х у + бвх у + бС х у) = с(и),

(30)

где в случае d = 1: U = R-lU, G(U) = R-lH(U),

в случае d = 2: у = R-lUR-l, G(U) = R-l H (U)R(31)

в случае & = 3: У = Я-1 х Я-1 х Я-1 х и, СЦи) = Я-1 х Я-1 х Я-1 х Н (и).

Решение (30) дается элементарными формулами:

^ о ■ для случая & = 1 компоненты вектора у определяются равенствами г)■ =-—

при d = 2 элементы матр ицы U определяются равенс твами Vjk = ^

1 - rdjA

1 - Т(djA + dkBУ ^ qm

при d = 3 элементы трехмерного массива у - v"i =-:—j—--.

^ jk 1 - Т(dA + dBB + d%)

Здесь gj, gj^ gjk компоненты элементы матрицы пли массива G (в

зависимости от размерности), &a-> dB>, d^ - собственные числа матриц Л, Б, C, j = 1,...,N, k = 1,...,K, m = 1,..., M. №1числения y решения регуляризованных

уравнений на текущей итерации находятся по формулам

d =1: U = RaV , d =2: U = RAV RB, d = 3 : U = RA x RB x RC x V. (32)

При решении необходимо следить, чтобы выполнялись условия

T = 1/dA, T=1/(dA + dkB), T=1/(dA + dkB + dm) Vj,k,m. (33)

Алгоритм решения задачи Дирихле для уравнений (6)

U

ющих в задаче. В случае нелинейных задач, имеющих множество решений, начальные значения нужно выбирать обоснованно.

NKM

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

H(U)

G(U).

4. По приведенным формулам вычисляем значения элементов массивов V и U (см. (32)) и считаем приближенные значения производных решения по формулам вида (14), (18), (21).

5. Если \\U — U\\/\\Uт|| < eS, то U - решение, иначе задаем U = U и вернемся к шагу 3. Здесь \\U\\ - максимальное значение элементов U.

Отметим, что доказательство сходимости такого процесса - отдельная задача для

Т

не только добиться сходимости при решении нелинейных задач с малыми параметрами, но и уменьшить число итераций в десятки и даже сотни раз. Спектральный портрет аппроксимирующих матриц при этом играет ключевую роль [7]. Интересные соображения по этому поводу читатель может найти также в [8].

N=K=M

бует, порядка p((2d + a(d + 1))Nd+l + O(Nd)) операций и s(3N2 + N + [a(d + 1) + 5}Nd + O(Nd~l)) байт оперативной памяти, где d - размерность задачи, а - коли-

чество различных производных от неизвестной функции, присутствующих в уравнении (кратные производные считаются за одну), p - количество итераций метода установления, s - объем памяти в байтах, занимаемый одним действительным, числом.

Доказательство. Поскольку арифметические операции, стоящие в правой части f (u, x), в рамках данного метода «превращаются> в соответствующие one рации с Nd элементами массива решения и его производных, то таких операций насчитывается O(Nd). Основные затраты времени при вычислении H(U) состоят в операциях по расчету производных, каждая из которых эквивалентна умножению матрицы на вектор (d = 1, N2 операций), матрицы на матрицу (d = 2, N3 операций) и слоистого массива на 3D массив (d = 3, N4 операций). Общее количество таких операций - aNd+1. Прочие затраты включают вычисление элементов G(U) по формулам (31) (dNd+1 операций), вычисление элементов у (O(Nd) операций), элементов U по формулам (32) (dNd+l операций), добавочных массивов и их произведений на аппроксимирующие матрицы по формулам вида (21) (dNd+1 + O(Nd) операций для каждой производной в уравнении). Суммируя все затраты, получаем (2d + a(d + 1))Nd+1 + O(Nd) операций на каждой итерации.

В ходе расчетов в оперативной памяти должны храниться

1) аппроксимирующая матрица A, матрицы ее спектрального разложения Ra(,b,c) и обратные к ним (в сумме 3 матрицы размера N х N с учетом того обстоятельства, что при N = M = K имеем A = B = C);

2) диагональные матрицы Da(,b,c) (суммарным объемом N чисел);

3) массивы значений производных в выражениях правой части (всего таких производных а — d), массивы решения на текущем и предыдущем временных слоях, массив правой части и массивы у, G, см. (31) (суммарный объем s(a — d + 5)Nd);

4) добавочные массивы для неизвестной функции и ее производных, реализующие граничные условия (в сумме d(a + 1) массивов);

5) прочие данные малого объема (параметры метода, коэффициенты уравнения).

Суммируя все затраты памяти, получаем выражение из условия теоремы. □

Отметим существенные преимущества разработанного метода перед КСМК, которую традиционно используют при реализации спектральных методов [9]. Такая схема приводит к системе линейных уравнений Ax = I) с разряженной матрицей, где x, b - векторы значений решения u(x) и значений u(x) — rf (u, x) в узлах интерполяции размера NKM при d = 3. Ее решение требует O(N3d) операций при N = K = M. В соответствии с теоремой 7 предложенная схема построения и решения СЛАУ типа уравнений Сильвестра и их тензорных обобщений (см. рис. 1) позволяет снизить количество операций на 1, 3 и 5 порядков для d = 1, 2, 3 соответственно.

4. Тестовые численные эксперименты

Разработанный метод был реализован на языке Java с использованием современных концепций объектно-ориентированного анализа проектирования и программирования. Все расчеты проводились на ЭВМ Intel Core ¡0-3330 CPU 3.00 GHz, DIMM DDE3 1600 MHz 8 Gb. Рассмотрим тестовую задачу

d2u/dx2 = u2e-x, u(—1) = e-1, u(1) = e. (34)

Ее точное решение - целая функция иех(х) = ех. С помощью описанного метода найдены приближенные значения решения и^ в узлах Интерполяции х■ ] = Пусть

и

ex,j

ue

Xх j )'

■■ spec

M

-l

max \uexj — и. \, M = max u,

j=1,..,N j=l,...,N

ex,j

(35)

На рис. 2 а изображена зависимость к^10 £зрес от количества узлов коллокации N. Из графика видно, что порядок аппроксимации решения является переменным, он стремительно растет с ростом N. Пунктирной линией показана зависимость логарифма правой части оценки (5) от N. Отметим, что сплошная и пунктирная линии имеют одинаковый угол наклона и выпуклость, что свидетельствует о том, что асимптотика погрешности построенного метода с высокой точностью совпадает с асимптотикой полиномиального приближения функции иех, т.е. свойство ненасыщаемости исполь-зовсшного приближения переносится на разработанный метод. Более подробно эта тема обсуждается ниже.

Для сравнения решение задачи (34) найдено с помощью метода конечных разностей (МКР) с использованием трехточечной аппроксимации второй производной и метода прогонки. Пусть по условию задачи требуется найти решение с точностью до шести значащих цифр (не менее). Для достижения точности £ = 10-6 (см. рис. 2 б) МКР требуется построить сетку из N = 197 узлов, в то время как разработанному методу достаточно N = 6 узлов. Данное обстоятельство, наблюдаемое также в 2Б и ЗБ задачах, определяет существенные преимущества разработанного алгоритма над МКР в случае задач с решениями высокой гладкости.

О

-10

-15

Y\logioe \ %

\ \ £ £ \ \ spec\ \ \ %

\ N

10

20

30

б)

Рис. 2. Зависимость погрешности численных решений (34) от N в логарифмической шкале: а) погрешность £.3рес', б) сравнение погрешностей МКР £й%ц и НМБН £зрес

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

и(±1) = ±1,

iv" = 6^|U|, Г

|v(±1) = 1, |w(±1) = ±1

и =2Хг, р = 6^' Ы' = 12хШ (36)

с соответствующими решениями ,2

х, х е (0,1]

и(х)

{

- х, х е [-1,0],

2 „ ^ п ь(х)

{

х

х е [-1,0],

х3, х е (0,1],

'ш(х)

{

х

х е [-1,0],

х4, х е (0,1].

Отметим, что и,ь,т являются 1, 2 и 3 раза непрерывно дифференцируемыми функциями соответственно. Кроме того, и' = ь" = т" ' = |х|, модуль непрерывности ш(|х|; д) = д. Поэтому в оценке (1) ш(и'; П) = ш(ь"; П) = ; П) = п и в Для

этих функций а = 1. Значит, порядки аппроксимации их наилучших полиномиальных приближений есть 3 и 4 соответственно.

В дополнение к одномерным задачам (36) рассмотрим еще двумерную и трехмерную нелинейные задачи

[ ихх + Цуу = 12ху(х2 + у2) Щ, I Щ(±1,у) = ±8§п (у)у\ у и(х, ±1) = ± ^п (х)хА,

2у

Ухх + Ууу + у,, = 12г 8ёп(ху)(х2 + у2) ^\У\ + —,

у(±1,у,г) = ±8^ (уг)у4г2,

У (х, ±1,г) = ± ^п (хг)х4г2, у(х, у, ±1) = ± (ху)х4у4

(37)

с решениями

и(х, у)

!

х

х е [-1,0]

х

х е (0,1]

X

{

у(х,у,г) = и(х,у) х

{

— г

- у4, у е [-1,0] у4, у е (0,1],

г е [-1, 0]

г2, г е (0,1].

иу являются 3 и 1 раз непрерывно дифференцируемыми. Обобще-

иу

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

Задачи (36), (37) решены численно с помощью разработанного метода, найдены относительные погрешности решений по формулам вида (35) и порядки аппроксимации. Для этого сделаны расчеты на последовательности сеток с двукратным увеличением числа узлов N = (К = М) = 4, 8,16,..., 128. Пуст ь е^рес(и) - относительная погрешность численного решения и на N узлах. Порядок аппроксимации алгоритма на N узлах найдем по формуле

Ям(и) = 1<Щ2(е%С / е^ес), N = 8,16,..., 128.

В табл. 1 приведены значения погрешностей е^,ес, в табл. 2 - значения по рядков К^ для пяти рассмотренных задач. Из приведенных данных видно строгое соответствие полученных порядков оценкам точности наилучших полиномиальных приближений. Этот факт свидетельствует об отсутствии насыщения предложенного алгоритма на основе метода коллокаций в рассмотренных нелинейных задачах.

Таблица 1

Погрешности численных решений

N £ spec(uU) £spec (v) £spec (w%) £spec(U ) £spec(V )

4 6,686E-03 l,731E-02 4,745E-03 2,945E-03 6Д00Е-04

8 2.173K-03 l,901E-03 3,33E-04 l,63E-04 2.89Ö [{-О!

16 5.372K-0I 2.323 K-0! l,893E-05 9,403E-06 1.29 !!•:-() !

32 1.282 K-0-1 2,89E-05 1Д79Е-06 7,803E-07 6,297E-05

64 ЗД65Е-05 3,608E-06 7,365E-08 5,754E-08 2,559E-05

128 7,888E-06 4,509E-07 I.602K-09 3,993E-09 _

Таблица 2

Порядки аппроксимации метода

N Rn (u) Rn (v) Rn (w) Rn (U) Rn (V)

8 1,4349 3,1867 3,8329 4,1753 1,0753

16 2,2026 3,0323 4,1369 4,1158 1,1617

32 2,0674 3,0071 4,0044 3,591 1,0391

64 2,0178 3,0017 4,0012 3,7614 1,299

128 2,0045 3,0004 4,0003 3,8491 _

Среднее 1,9454 3,0456 3,9951 3,8985 1,1438

Заключение

В заключение отметим, что каждая нелинейная проблема уникальна и требует отдельного анализа количества ее решений, их свойств и сходимости метода установления к тому или иному решению. Помимо решений с различным порядком гладкости, могут встречаться решения с существенными особенностями в комплексной плоскости (см. теорему 3). Этот случай здесь подробно не обсуждался в силу ограничений по объему статьи. Исследование устойчивости стационарного решения регуляризованно-го уравнения, сходимости метода установления и, конечно же, насыщаемости метода должно проводиться отдельно для каждой задачи. Тем не менее косвенное обоснование ненасыщаемости, приведенное выше, говорит о принципиальной возможности реализовать свойства наилучших полиномиальных приближений при решении сложных прикладных задач предложенным методом. Реализация таких свойств в рамках разработанного алгоритма, в свою очередь, обеспечивает минимизацию вычислительных затрат при сохранении необходимой точности решения и существенные преимущества над традиционными схемами МКР и КСМК.

Работа выполнена за счет гранта Российского научного фонда (проект № 17-7110135).

Литература

1. Бабенко, К.И. Основы численного анализа / К.И. Бабенко. - М.: Наука. Гл. ред. физ.-мат. лит., 1986.

2. Бабенко, К.И. О явлении насыщения в численном анализе / К.И. Бабенко // Доклады АН СССР. - 1978. - Т. 241, № 3. - С. 505-508.

3. Семисалов, Б.В. Нелокальный алгоритм поиска решений уравнения Пуассона и его приложения / Б.В. Семисалов // Журнал вычислительной математики и математической физики. - 2014. - Т. 54, № 7. - С. 1110-1135.

4. Trefethen, L.N. Approximation Theory and Approximation Practice / L.N. Trefethen. -Philadelphia: SIAM, 2013.

5. Дзядык, В.К. Введение в теорию равномерного приближения функций полиномами /

B.К. Дзядык. - М.: Наука, 1977.

6. 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 Scientific Computing. - 2009. - V. 31, № 3. - P. 2015-2046.

7. Белов, А.А. Эволюционная факторизация и сверхбыстрый счет на установление / А.А. Белов, Н.Н. Калиткин // Математическое моделирование. - 2014. - Т. 26, № 9. -

C. 47-64.

8. Белых, В.Н. Особенности реализации ненасыщаемого численного метода для внешней осесимметричной задачи Неймана / В.Н. Белых // Сибирский математический журнал. - 2013. - Т. 54, № 6. - С. 1237-1249.

9. Boyd, J. Chebyshev and Fourier Spectral Methods / J. Boyd. - Mineola; N.Y.: University of Michigan, 2000.

Борис Владимирович Семисалов, кандидат физико-математических наук, Институт вычислительных технологии СО РАН; Новосибирский государственный университет (г. Новосибирск, Российская Федерация), [email protected].

Поступила в редакцию 15 ноября 2017 г.

MSC 65N35, 35J60 DOI: 10.14529/mmpl80210

DEVELOPMENT AND ANALYSIS OF THE FAST PSEUDO SPECTRAL METHOD FOR SOLVING NONLINEAR DIRICHLET PROBLEMS

В. V. Semisalov, Institute of Computational Technologies SB EAS, Novosibirsk; Russian Federation, Novosibirsk State University, Novosibirsk, Russian Federation, [email protected]

Numerical method for solving one-, two- and three-dimensional Dirichlet problems for the nonlinear elliptic equations has been designed. The method is based on the application of Chebyshev approximations without saturation and on a new way of forming and solving the systems of linear equations after discretization of the original differential problem. Wherein the differential operators are approximated by means of matrices and the equation itself is approximated by the Sylvester equation (2D case) or by its tensor generalization (3D case). While solving test problems with the solutions of different regularity we have shown a rigid correspondence between the rate of convergence of the proposed method and the order of smoothness (or regularity) of the sought-for function. The observed rates of convergence strictly correspond to the error estimates of the best polynomial approximations and show the absence of saturation of the designed algorithm. This results in the essential reduction of memory costs and number of operations for cases of the problems with solutions of a high order of smoothness.

Keywords: Chebyshev approximation; boundary-value problem; nonlocal method without saturation; stabilization method.

Acknowledgements. This research was done under the financial support of Russian Science Foundation (project No. 17-71-10135).

B.B. CeMHcajiOB

References

1. Babenko K.I. Osnovy chislennogo analiza [Fundamentals of Numerical Analysis]. Moscow, Nauka, 1986. 714 p.

2. Babenko K.I. [On the Saturation Phenomenon in Numerical Analysis]. Dokl. Akad. Nauk SSSR [Doklady Mathematics], 1978, vol. 241, no. 3, pp. 505-508. (in Russian)

3. Semisalov B.V. Non-Local Algorithm of Finding Solution to the Poisson Equation and Its Applications. Computational Mathematics and Mathematical Physics, 2014, vol. 54, no. 7, pp. 1110-1135. (in Russian)

4. Trefethen L.N. Approximation Theory and Approximation Practice. Philadelphia, SIAM, 2013.

5. Dzjadyk V.K. Vvedenie v teoriju ravnomernogo priblizhenija funkcij polinomami [Introduction to the Theory of Uniform Approximation by Polynomials]. Moscow, Nauka, 1977.

6. 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 Scientific Computing, 2009, vol. 31, no. 3, pp. 2015-2046. DOI: 10.1137/070706537

7. Belov A.A., Kalitkin N.N. Evolutionary Factorization and Superfast Relaxation Count. Mathematical Models and Computer Simulations, 2014, vol. 26, no. 9, pp. 47-64. (in Russian)

8. Belykh V.N. Particular Features of Implementation of an Unsaturated Numerical Method for the Exterior Axisymmetric Neumann Problem. Siberian Mathematical Journal, 2013, vol. 54, no. 6, pp. 984-993. DOI: 10.1134/S0037446613060037

9. Boyd J. Chebyshev and Fourier Spectral Methods. Mineola, N.Y., University of Michigan, 2000.

Received November 15, 2017