CC BY
30
Вычислительные технологии
Том 11, часть 2, Специальный выпуск, 2006
РЕШЕНИЕ ЭЛЛИПТИЧЕСКИХ ЗАДАЧ С ОСОБЕННОСТЯМИ ПО СХЕМАМ ВЫСОКОГО ПОРЯДКА АППРОКСИМАЦИИ*
В. П. ШАПЕЕВ
Институт теоретической и прикладной механики СО РАН
Новосибирск, Россия e-mail: shapeev@itam.nsc.ru А. В. ШАПЕЕВ Национальный университет Сингапура e-mail: shapeev@gorodok.net
Numerical solutions of boundary-value problems with singularities are examined. We consider finite difference schemes with high order of approximation for two specific different elliptic equations, namely Poisson equation and an equation with a small parameter at the highest derivative. For the boundary-value problem for the Poisson equation the singularity arises due to a discontinuity of the second-order derivative at the corner point of the domain boundary. For the equation with a small parameter at the highest derivative an inner thin boundary layer is a manifestation of the singularity. Behavior of the solution and convergence of the numerical solution on a sequence of grids has been analyzed using exact solutions. In both cases the convergence order of the high-order methods was higher compared to the low-order methods. A derivation of the finite difference schemes and analysis of the solutions were carried out with the help of computer algebra software package Mathematica.
Введение
Рассматриваются две краевые задачи с особенностями для двух различных эллиптических уравнений. На основе численных экспериментов анализируются порядок сходимости на последовательности сеток и поведение погрешности численного решения, получаемого с помощью разностных схем различного порядка аппроксимации. Известно, что в случае наличия особенностей в виде сильных разрывов в области решения краевых задач для гиперболических уравнений не удается добиться высокого порядка сходимости численного решения по схемам высокого порядка аппроксимации. Проведенные в данной работе численные эксперименты при решении краевых задач для эллиптических уравнений показывают, что при наличии особенностей в решении эллиптических задач по схемам высокого
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00080-а).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
порядка на последовательности сеток может наблюдаться высокий порядок сходимости численного решения. При этом, естественно, порядок сходимости может не совпадать с порядком аппроксимации и зависит от вида особенности и порядка разностной схемы. Интересно, что и характер поведения погрешности различный для схем различного порядка.
Для построения схем высокого порядка здесь использован метод неопределенных коэффициентов, реализованный в системе компьютерной алгебры МаШеша^са, как это описано в [1]. В развитие этого способа в [2] предложено в случае высокого порядка разностной схемы на конкретной сетке строить ее численно по программе на ЭВМ, затем полученные формулы передавать в другую часть программы (или в другую программу) и решать, как обычно, численно рассматриваемую краевую задачу. Таким способом для конкретных эллиптических уравнений на регулярных и нерегулярных шаблонах удалось построить схемы до восемнадцатого порядка включительно и решить по ним соответствующие краевые задачи. Численные эксперименты по решению различных краевых задач без особенностей для эллиптических уравнений с переменными коэффициентами в области с криволинейной границей на последовательности сеток показывали зачастую соответствие порядка сходимости численного решения порядку аппроксимации [2]. Здесь рассмотрим решение задач с особенностями.
1. Краевая задача с разрывом второй производной на границе для уравнения Пуассона
Рассматривается краевая задача с условиями Дирихле для уравнения Пуассона
Ди(х, у) = 1, 0 < х < 1, 0 < у < 1, (1)
и|г = 0 (1)
в единичном квадрате с однородными краевыми условиями. Здесь Г — граница области. Эта задача имеет особенность: в угловых точках квадрата вторые производные терпят разрыв.
Точное решение этой задачи можно выписать в виде равномерно и абсолютно сходящегося ряда
( ) = -— ^^ ап((2г - 1) пх) 8т((2; - 1) пу)
п4 ¿1 (2г - 1) (2; - 1) ((2г - 1)2 + (2; - 1)2)' ( )
График точного решения в расчетной области представлен на рис. 1. Для оценки погрешности разностного решения в качестве точного решения взят начальный отрезок ряда так, чтобы соответствующий остаток ряда по модулю гарантированно не превосходил 5.96 ■ 10-10. Модуль остатка ряда оценивался с помощью мажоранты
1 „ оо оо 1
16 1
п4 ¿1 2=1 (2г - 1) (2; - 1) ((2г - 1)2 + (2; - 1)2).
Для задачи (1) построены схемы 2-го, 6-го и 10-го порядков аппроксимации [2]. На рис. 2 приведены графики погрешности разностного решения для схемы второго порядка
40
Рис. 1. Точное решение уравнения Пуассона.
Рис. 2. Погрешность схемы второго порядка.
слева и ее диагонального профиля справа. На рис. 3 показана погрешность разностного решения для схемы десятого порядка. Числами 10, 20, 30,40 обозначены номера узлов сетки вдоль осей х, у.
Из расчетов видно, что максимум погрешности схемы второго порядка находится в центре расчетной области и примерно в 2.5 раза больше погрешности решения в угловой внутренней точке, где она минимальна. Погрешность схемы десятого порядка ведет себя по-другому: максимум погрешности достигается в угловых узлах сетки внутри области и гораздо больше (в О(N¡2 ) раз) погрешности в центре расчетной области (здесь — количество узлов сетки).
В табл. 1 приведена погрешность разностного решения Я(Ь.) в равномерной и среднеквадратичной нормах для схем 2-го, 6-го, 10-го порядков на сетках 10 х 10, 20 х 20 и 40 х 40. Из таблицы видно, что все рассмотренные схемы имеют второй порядок сходимости в равномерной норме, но величина погрешности у схемы шестого порядка в 40 раз меньше, чем у схемы второго порядка, а у схемы десятого порядка в свою очередь в 20 раз меньше, чем у схемы шестого порядка. Разностное решение в среднеквадратичной норме для схем
Рис. 3. Погрешность схемы десятого порядка.
Таблица 1. Погрешность разностного решения Д(Н) на последовательности сеток
Порядок схемы к Шаг сетки Н Р1и Р1|2
2 1/10 5.73 ■ 10-4 3.10 ■ 10-4
2 1/20 1.45 ■ 10-4 8.22 ■ 10-5
2 1/40 3.63 ■ 10-5 2.11 ■ 10-5
6 1/10 1.50 ■ 10-5 3.28 ■ 10-6
6 1/20 3.73 ■ 10-6 4.16 ■ 10-7
6 1/40 9.33 ■ 10-7 5.28 ■ 10-8
10 1/10 6.99 ■ 10-7 1.72 ■ 10-7
10 1/20 1.74 ■ 10-7 2.15 ■ 10-8
10 1/40 4.34 ■ 10-8 2.71 ■ 10-9
шестого и десятого порядков аппроксимации сходится с третьим порядком, а поточечно (т. е. в каждой фикисированной точке) — с четвертым порядком.
Здесь наличие локальной особенности в решении эллиптической задачи по схеме повышенного порядка аппроксимации понизило во всей области порядок сходимости разностного решения по отношению к порядку аппроксимации. Но он остался выше порядка сходимости решения, полученного по схеме низкого порядка аппроксимации. Из теории эллиптических уравнений известно, что решение в каждой точке области зависит от всех значений в условиях на границе области. Поэтому отсутствие гладкости в нескольких точках в условиях на границе сказывается на характере поведения решения задачи во всей области. Однако в данном случае схема повышенного порядка аппроксимации точнее воспроизводит решение, чем схема низкого порядка.
2. Краевая задача для уравнения с малым параметром
Далее рассматривается краевая задача в единичном квадрате с границей Г с условиями Дирихле для уравнения типа диффузии-конвекции с малым параметром в при старших
производных:
^ ,'д2u д2u
Я — +
дх2 dy2 u|r = v(x,y)-
du ,лди
+ (x - a) dx + (y - b) dy
0 < x < 1, 0 < y < 1,
(3)
(4)
Функция v(x, y) в краевых условиях модельной задачи взята из точного решения. Уравнение (3) имеет следующее семейство точных решений:
u(x, y) = G((x — a) cosф — (y — b) sin ф) G((x — а) sin ф + (y — b) cosф),
G(£) = 2 +1 / e-t2dt. 0
Решения (4) имеют ярко выраженный внутренний погранслой с большими градиентами решения. Толщина погранслоя порядка \[в.
За недостатком места в качестве примера здесь приведем схему только для уравнения с постоянными коэффициентами
д2 u д2 u du du
Т"2 + + + d^- =0.
dx2 dy2 дх dy
Рассмотрим аппроксимацию на сетке Gh U rh = {(ih, jh), 0 < i, j < N}. Обозначим
. = w¿+1,j — Wi-ij . = Wi,j+i —
Axwi,j = 0, , Ay wi,j =
2h
2h
л = Wi+i,j — 2w i,j + Wi-i,j . = Wi,j+i — 2wi,j + Wi,j-i AxWi,j = h2 ' луWi,j = h2
Тогда схема будет иметь вид Axu + Ay u + с Axu + d Ay u + c d h2
1
+
1
12 + c2 h2 12 + d2 h2
c2 h2
Ax Ayu + —¡^ Ax u+
12 + c2 h2
+
d2 h2 dh2 (144 + 6(7c2 + d2) h2 + c2 (c2 + 2 d2) h4) л 4
Ay u +-^-„ Л . 2,2' ^ . ---- Ay Axu+
12 + d2 h2'^ ' 6 (12 + c2 h2) (12 + d2 h2)
c h2 (144 + 6 (c2 + 7 d2) h2 + d2 (2 c2 + d2) h4)
+
6 (12 + c2 h2) (12 + d2 h2) (72 h2 + 9 (c2 + d2) h4 + c2 d2 h6)
Ax Ayu+
3(12 + c2 h2) (12 + d2 h2) Ax Ayu 0'
Далее для параметров в (3), (4) брались следующие значения: а = 0.6, Ь = 0.3, ф = п/6. На рис. 4 показаны слева график точного решения для в = 10-4 и справа профили точного решения этой задачи для в = 10-2, 10-3,10-4 в сечениях, перпендикулярных к фронту погранслоя. Графики профилей, соответствующие указанным значениям в по порядку их перечисления, помечены индексами 1, 2 и 3.
Проведены численные эксперименты по решению модельной задачи (3) по схемам второго и четвертого порядков [2]. При в = 10-2, в = 10-3 на последовательности равномерных квадратных сеток к = 1/20, 1/40, 1/80, 1/160 наблюдается сходимость разностного
0
+
решения с порядком, соответствующим порядку аппроксимации. В табл. 2 и 3 приведены результаты этих численных экспериментов.
Рис. 4. Решение задачи с малым параметром. Таблица 2. Погрешность разностного решения Я(к) на последовательности сеток при в = 10-2
Порядок схемы к Шаг сетки Н ця|и 1|Я»2
2 1/10 3.15 • 10-2 9.79 • 10-3
4 1/10 4.66 • 10-3 1.13 • 10-3
2 1/20 7.71 • 10-3 2.35 • 10-3
4 1/20 4.02 • 10-4 7.92 • 10-5
2 1/40 1.95 • 10-3 5.90 • 10-4
4 1/40 2.54 • 10-5 5.60 • 10-6
2 1/80 4.87 • 10-4 1.49 • 10-4
4 1/80 1.59 • 10-6 3.69 • 10-7
2 1/160 1.22 • 10-4 3.73 • 10-5
4 1/160 1.00 • 10-7 2.35 • 10-8
Таблица 3. Погрешность разностного решения Я(Н) на последовательности сеток при в = 10 3
Порядок схемы к Шаг сетки Н ця|и 1|Я»2
2 1/10 6.31 • 10-1 1.54 • 10-1
4 1/10 1.62 3.71 • 10-1
2 1/20 1.52 • 10-1 3.27 • 10-2
4 1/20 3.65 • 10-2 7.36 • 10-3
2 1/40 4.28 • 10-2 7.46 • 10-3
4 1/40 2.54 • 10-3 3.87 • 10-4
2 1/80 1.11 • 10-2 1.78 • 10-3
4 1/80 1.59 • 10-4 2.21 • 10-5
2 1/160 2.76 • 10-3 4.40 • 10-4
4 1/160 1.00 • 10-5 1.16 • 10-6
Таблица 4. Погрешность разностного решения Д(Л) на последовательности сеток при в = 10-4
Порядок схемы k Шаг сетки h IIRU IIRII2
2 1/20 1.67 4.59 ■ 10-1
4 1/20 14.2 1.90
2 1/40 3.57 ■ 10-1 7.28 ■ 10-2
4 1/40 1.90 4.08 ■ 10-1
2 1/80 1.83 ■ 10-1 2.49 ■ 10-2
4 1/80 3.32 ■ 10-2 3.61 ■ 10-3
2 1/160 7.43 ■ 10-2 7.23 ■ 10-3
4 1/160 3.22 ■ 10-3 2.62 ■ 10-4
В работе [3] показаны хорошие возможности метода коллокаций и наименьших квадратов для решения подобных задач. В частности, в ней для ß = 10-4 на адаптивной сетке достигнута точность 10-2. На равномерной сетке с шагом h = 1/160 установлено, что схема второго порядка имеет в равномерной норме погрешность 7.43 ■ 10-2, а схема четвертого порядка — погрешность 3.22 ■ 10-3. В среднеквадратичной норме они имеют погрешности соответственно 7.23 ■ 10-3 и 2.62 ■ 10-4. То есть схема четвертого порядка имеет существенные преимущества перед схемой второго порядка и при относительно умеренном шаге сетки дает решение задачи с высокой точностью. Результаты численного эксперимента на последовательности сеток при ß = 10-4 приведены в табл. 4. В этом случае для наблюдения сходимости с порядком, соответствующим порядку аппроксимации, требуется дальнейшее мельчение шага. Но сходимость с меньшим порядком начинается раньше: с момента, когда шаг сетки укладывается в толщину погранслоя.
Очевидно, что сложность построения программной реализации схем высокого порядка аппроксимации не является их достоинством по сравнению со схемами низкого порядка аппроксимации. Однако из приведенных примеров следует, что при решении задач для эллиптических уравнений в условиях ограниченности ресурсов (память, время) схемы высокого порядка имеют преимущества перед схемами низкого порядка как при решении задач без особенностей, так и при решении задач с особенностями. Они могут дополнять, а иногда и заменять другие приемы достижения большей точности (применение адаптивных сеток, многосеточных алгоритмов и др.). Сложности, которые возникают при построении и исследовании схем повышенного порядка точности [4,5], можно преодолевать с помощью систем компьютерной алгебры [6-8].
Список литературы
[1] Вллиуллин А.Н., Ганжа В.Г., Ильин В.П. и др. Задача автоматического построения и исследования на ЭВМ разностных схем в аналитическом виде // Докл. АН СССР. 1984. Т. 275, № 3. С. 528-532.
[2] Шапеев А.В., Шапеев В.П. Разностные схемы повышенной точности для решения эллиптических уравнений в области с криволинейной границей // Журн. вычисл. математики и мат. физики. 2000. Т. 40, № 2. С. 223-232.
[3] Слепцов А.Г., Шокин Ю.И. Адаптивный проекционно-сеточный метод для эллиптических задач // Журн. вычисл. математики и мат. физики. 1997. Т. 37, № 5. С. 572-586.
[4] Вллиуллин А.Н. Схемы повышенной точности для задач математической физики. Новосибирск: Изд-во НГУ, 1973. 138 с.
[5] Lele S.K. Compact finite difference schemes with spectral-like resolution //J. Comput. Phys. 1992. Vol. 102, N 1. P. 16-42.
[6] Мазурик С.И., Шапеев В.П. Применение символьных преобразований на ЭВМ для исследования аппроксимации и устойчивости разностных схем // Журн. вычисл. математики и мат. физики. 1986. Т. 26, № 4. С. 586-600.
[7] Steinberg S. A problem solving environment for numerical partial differential equations // 6th IMACS Intern. Conf. on Applications of Computer Algebra. Abstracts. St. Petersburg, Russia, June 25-28, 2000. P. 98-99.
[8] Shapeev A.V. Application of computer algebra systems to construct high-order difference schemes // 6th IMACS Intern. Conf. on Applications of Computer Algebra. Abstracts. St. Petersburg, Russia, June 25-28, 2000. P. 92-93.
Поступила в редакцию 4 апреля 2006 г.
