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

Вычислительные технологии

Том 10, № 2, 2005

ПРИМЕНЕНИЕ МЕТОДА НЬЮТОНА В РАСЧЕТНЫХ СХЕМАХ ВЫСОКОГО ПОРЯДКА ДЛЯ СТАЦИОНАРНЫХ ЗАДАЧ*

Н.Б. Петровская Институт прикладной математики им. М. В. Келдыша РАН,

Москва, Россия e-mail: PetrovskayaN@yandex.ru

High order discontinuous Galerkin discretization schemes are considered for steady state problems. We discuss the issue of oscillations arising when Newton's method is applied to obtain a steady state solution. It is shown that flux approximation near flux extrema may produce spurious oscillations propagating over the computational domain. The control over the numerical flux in the problem allows us to obtain non-oscillating convergent solutions.

Введение

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

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

*Работа выполнена при финансовой поддержке компании "Боинг" (контракт 104AE) и Российского фонда фундаментальных исследований (грант № 03-01-00063).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

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

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

1. Формулировка задачи

Мы рассматриваем краевую задачу для функции и (ж) на отрезке [0,1]

= 0, Ви(ж) = 0, ж € П=[0,1], (1)

где ^(и(ж)) — функция потока, а обозначение В используется для оператора граничного условия. Для численного решения краевой задачи введем расчетную сетку О следующим образом:

N

О =и = [Жг,Жг+1], 1 < г < N

г=1

где жг — координата узла сетки.

В качестве примера схемы высокого порядка будет рассмотрен разрывный метод Га-леркина [2]. Метод определяет приближенное решение и^(ж) в каждой ячейке сетки как

к

иь(ж) = ^икфк(ж), к = 0,1,...,К, х € бг, (2)

к=0

где базисные функции выбраны полиномами соответствующей степени: фк (ж) = ((ж - Хг )/(Хг+1 - Хг))к, к = 0, 1,...,К, X € бг.

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

^(и^(жг+1 ))фк (жг+1) - ^(ин(жг))фк (ж*)- ( ^(ж,иЛ(ж))¿ж = 0, к = 0,1,...,К. (3) I аж

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

^ига)(ига+1 - ип) = -Щи"), (4)

где невязка метода И.(и), которая задается уравнениями (3) в каждой ячейке сетки, и якобиан ^и) = [дИ,/ди] зависят от решения и" на п-й итерации.

2. Численный поток для схем высокого порядка

Определение численного потока в уравнениях (3) требует знания точек экстремума функции Г (п). Вычисление производной ¿Г/ёп, в свою очередь, зависит от того, как определена вариация решения. В отличие от кусочно-постоянной аппроксимации решения, для схемы высокого порядка (К > 0) приближенное решение меняется в области [£¿,£¿+1]. Вариация решения 5пн = п^(хг+1 — 0) — п^(хг + 0) может привести к появлению экстремума потока во внутренней точке ячейки сетки, в то время как поток остается монотонной функцией в узлах ячейки. Таким образом, аппроксимация потока в схеме высокого порядка не будет правильной до тех пор, пока не определено правильное положение экстремума. В качестве примера рассмотрим численный поток Энквиста — Ошера [3]

иг щ

РЕ°(щ,Пг) = У тт(Г'(в), 0)ёв + у'тах(Г/(з), 0)ёв + Г(0), (5)

0 0

где и = пь(хг — 0) и пг = п^(хг + 0) — левое и правое значения решения в узле сетки.

Пример аппроксимации потока в точке экстремума показан на рис. 1. Пусть функция Г (п) имеет экстремум в точке и = и* ив узлах ячейки выполнены условия п^(хг ± 0) > и*,

Щ(х)

1. "г-1 и(х) = и*

-

- \dFiuydu = 0

-......... ......... ..........

х1 +1

Щ -

и)/du = 0

.......... .........

щ{х)

и(х) = и* и1 + 1(х)

.......

X,- + 1 X

Рис. 1. Аппроксимация потока вблизи точки экстремума для численного потока (5): а — кусочно-постоянное решение помещает экстремум потока на интерфейс сетки. Экстремум учтен при аппроксимации потока; б — реконструкция решения полиномами высокого порядка помещает экстремум во внутреннюю точку ячейки сетки. Определение численного потока не учитывает присутствие экстремума в ячейке.

± 0) < и*. Тогда в соответствии с определением (5) значения решения в узлах сетки, необходимые для аппроксимации потока, схематически помечены черными точками. Для кусочно-постоянного восполнения решения, которое показано на рис. 1, а, экстремум потока на интерфейсе сетки учтен при аппроксимации потока. Аппроксимация решения более высокого порядка (например, кусочно-линейное восполнение решения, представленное на рис. 1, б) не принимает во внимание экстремум потока во внутренней точке ячейки сетки. Таким образом, решение в ячейке e^ не участвует в определении численного потока.

3. Метод Ньютона для схем высокого порядка

Покажем, что "фантомное" решение в ячейке сетки может привести к сингулярному якобиану системы (4) в случае, когда для решения задачи используется метод Ньютона. Определим локальный якобиан в ячейке e^ как jk1k2 = dRk1 /duk2, где локальные индексы = 0,1,...,K используются, чтобы обозначить компоненты вектора невязки и решения в данной ячейке.

Каждое из уравнений (3), рассмотренное для показателя степени полинома k > 0, содержит интегральное слагаемое, чья линеаризация гарантирует ненулевые коэффициенты в соответствующей строке локального якобиана. Однако первое из уравнений (3) (дискретный закон сохранения) требует только равенства потоков на интерфейсах ячейки. Таким образом, если определение численного потока касается только решения в соседних ячейках сетки, в матрице jk1k2 появится нулевая строка. Рассмотрим теперь блок якобиана Jim

l е L = 1,2,... ,N(K + 1) + 1, m G M = M1,M1 + 1,... ,M1 + K,

где N — число ячеек сетки, а индекс M1 определен для ячейки e^: M1 = (i — 1)(K + 1) + 1. Выберем произвольное m1 из множества M. Поскольку ранг локального якобиана есть rank(jk1 k2) = K < dim(jk1k2) = K + 1, перестановкой строк и столбцов матрицы Jim мы получаем

Jim1 =0 Vl G L : M1 < l < M1 + K.

Если при этом

Jlm1 = 0 Vl G L : l < M1 или l > M1 + K

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

4. Метод установления для решения стационарных задач

Рассмотрим скалярный закон сохранения

ди

— + FX(x,u(x)) = 0, x G П. (6)

Пусть un и un+1 — решения дискретной задачи на данной сетке на соответствующих временных слоях tn и tn+1 = tn + At. Неявная схема для гиперболического уравнения (6) может быть записана как

M(un+1 — un) = — AtRDG(un+1),

где — невязка стационарной задачи (3), а матрица М определена как

Мы = / фк(ж)фг(ж)^г, к, I = 0,1,..., К,

в каждой ячейке сетки.

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

^ига+1 - и") = -Щи"). Якобианом этой системы является

' = М +

Дг

где JDG — якобиан стационарной задачи, = [5И^/ди].

Таким образом, присутствие в дискретизации положительно определенной диагональной матрицы М обеспечивает наличие диагональных членов в матрице линеаризованной задачи, даже если якобиан стационарной задачи (3) является сингулярным.

e

г

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

В качестве численного примера была рассмотрена следующая задача, известная как задача о расходе массы в канале с переменным сечением [4].

Пусть A(x) — площадь сечения канала, A(x) = 1/2 + 2(x — 1/2)2; u(x) — скорость.

Расход массы в канале описывается уравнением

(x,u(x)) = d(A(x)m(u)) = 0 _

dx " dx ' Х ^ 1J' (7)

где поток массы определяется как m(u) = ц(1 — u2). Дополнительным условием является

1

/u(x)dx = B' (8)

о

где постоянная B определяет тип решения. Выбор B = —1/4 позволяет получить разрывное решение, показанное на рис. 3, а,

и() = f —V(1 — 1/2A(x)), 0 < x < xs, orxsh + 0 < x < 1, (9)

(x) 1 V(1 — 1/2A(x)), xs < x < xsh — 0, ( )

где xsh — положение разрыва, а точка xs = 1/2 определяет экстремум потока: u(xs) = us = 0.

Для численного решения задачи (7), (8) использовалась схема (3), (4). Решение ищется на последовательности однородных сеток, каждая из которых получается путем уменьшения шага предыдущей сетки в два раза. Для получения приближенного решения на каждой сетке решается методом Ньютона система нелинейных уравнений (3). При этом начальное приближение для метода Ньютона на первой сетке выбрано u0(x) = const = —1.0.

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

Проведенные расчеты показали, что метод Ньютона расходится при аппроксимации решения полиномами порядка K > 0. Только кусочно-постоянное восполнение решения позволило получить сходящиеся решения на последовательности сеток. В то же время использование метода Ньютона для аппроксимаций высокого порядка приводит к осцил-ляциям численного решения уже на первой сетке, что делает невозможными вычисления на более мелких сетках.

Развитие осцилляций на начальной сетке показано на рис. 2, где римскими цифрами обозначено промежуточное решение на соответствующей итерации метода Ньютона. В обозначениях рисунка промежуточное решение u1 есть начальное приближение к методу Ньютона, u1 = const = -1.0. Видно, что если на первых ньютоновских итерациях поток остается монотонной функцией в области определения решения, то промежуточное решение u111 уже требует аппроксимации потока в точке экстремума. Аппроксимация потока в ячейке, содержащей экстремум u = 0, приводит к сингулярному якобиану и, как следствие, к нефизическим осцилляциям решения uIV на следующей итерации метода Ньютона. Заметим, что осцилляции в приближенном решении возникают еще до формирования разрыва. Появление осцилляций в промежуточном решении uIV связано с необходимостью аппроксимировать поток в точке экстремума u = 0 для решения uIII на предыдущей итерации метода Ньютона.

Полученные численные результаты подтверждают, что при использовании метода Ньютона в схеме высокого порядка требуется контроль аппроксимации потока. Необходимо заметить, что при использовании метода Ньютона в промежуточном решении могут возникнуть нефизические экстремумы, которые должны быть устранены на каждой итерации, чтобы не допустить неправильной аппроксимации потока. Предложенный нами в [5] алгоритм проверяет наличие экстремума потока во внутренней точке каждой ячейки сетки. При обнаружении "фантомного" решения в ячейке сетки в этой ячейке на следующей итерации метода Ньютона используется кусочно-постоянная аппроксимация решения.

Рис. 2. Численное решение задачи (7), (8) на начальной сетке. Возникновение осцилляций на очередной итерации метода Ньютона при аппроксимации решения полиномами степени К > 0.

Рис. 3. Численное решение задачи (7), (8): а — точное и приближенное (степень полинома К = 3) решение на грубой (число ячеек сетки N = 16) и мелкой (число ячеек сетки N = 128) сетках. Точное решение показано сплошной линией, приближенное решение — пунктирной; б — сходимость приближенного решения (степень полинома К > 0) на последовательности равномерных сеток.

Алгоритм позволяет получить сходящиеся решения для схем высокого порядка (К > 0). Пример приближенного решения показан на рис. 3, а на грубой и мелкой сетках для аппроксимации кубическими полиномами. Из рисунка видно, что предложенный алгоритм понижает порядок аппроксимации в двух смежных ячейках сетки вблизи разрыва, для того чтобы избежать появления на нем "фантомного" решения. Однако в областях, где решение гладкое, порядок аппроксимации сохранен.

Сходимость приближенного решения на последовательности равномерных сеток показана на рис. 3, б. Число узлов начальной сетки выбрано N = 8. Норма ошибки

1

||егг||^1 = J |и(х) -

о

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

Заключение

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

должно быть уделено пристальное внимание, так как плохо обусловленная матрица системы (4) может появиться на последних итерациях, где шаг по времени выбирается обычно

достаточно большим для того, чтобы ускорить сходимость.

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

[1] Куликовский А.Г., ПогорЕлов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с.

[2] Cookburn B., Karniadakis G.E., Shu C.-W. The development of discontinuous Galerkin methods / / Discontinuous Galerkin methods. Theory, Computation and Applications / B. Cockburn, G.E. Karniadakis, C.-W. Shu (Eds). Lecture Notes in Comput. Sci. Engrg. N.Y.: Springer-Verl., 2000. Vol. 11. P. 3-50.

[3] Engquist B., Osher S. One-sided difference equations for nonlinear conservation laws // Math. Comp. 1981. Vol. 36. P. 321-352.

[4] Anderson J.D. (Jr.) Fundamentals of Aerodynamics. N.Y.: McGraw-Hill, 1991.

[5] Петровская Н.Б. Метод Ньютона для схем высокого порядка: почему он не работает? Москва, 2004. (Препр. ИПМ РАН).

Поступила в редакцию 25 мая 2004 г.