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

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

Том 13, № 3, 2008

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

А.И. Задорин Омский филиал Института математики СО РАН, Россия e-mail: [email protected]

For a function of two variables, corresponding to a solution of the elliptic problem with a boundary layer, the interpolation function is constructed. It is shown, that interpolation error is small for any relations between a small parameter and mesh steps. It is proved, that in case of the Shishkin mesh polinomial interpolation the method is uniformly accurate with regard to a small parameter. Formulas for calculation of derivatives, that use mesh values of the solution with a uniformly small error are proposed.

Введение

Для эллиптических задач с пограничным слоем разностные схемы со свойством равномерной сходимости, как известно, могут быть построены либо на неравномерной сетке, сгущающейся в пограничном слое [1, 2], либо на равномерной — с подгонкой схемы к погранслойной составляющей решения [3]. Если решение в узлах сетки найдено, то актуален вопрос интерполяции сеточного решения. Как будет показано ниже, полиномиальная интерполяция может приводить к значительным погрешностям. В данной работе в случае произвольной прямоугольной сетки строится и обосновывается интерполяция, учитывающая погранслойную составляющую в решении. Показано, что в случае сгущающейся сетки [4] полиномиальная интерполяция дает равномерно малую погрешность. Рассмотрены и вопросы нахождения производных в произвольной точке области через значения функции в узлах сетки.

Итак, пусть функция и(х,у), определенная в прямоугольной области

В = {(х,у) : 0 < х < Ь1, 0 < у < Ь2],

представима в виде

u(x,y) = y(y) exp(—a0e -x) + P(x,y), где |Y(j)(y)| < C, 0 < j < 3, e > 0, a0 > 0,

Qi+j

dxldy-

P (x,y)

<C

1

—- exp(—ae x) + 1

i-1

dj

—-P (x, y) dyj v 'y>

< C, j < 3.

(1)

0 < i < 3, 0 < j < 1,

*Работа выполнена при финансовой поддержке программы 1.3 ОМН РАН и Национального фонда Болгарии, проект HS-MI-106/2005.

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

Под С и С, подразумеваем положительные постоянные, не зависящие от е и шагов сетки. Предполагаем, что а > 0. Функция вида (1), с экспоненциальной погранслойной составляющей, соответствует решению эллиптической задачи с регулярным пограничным слоем у границы х = 0 в двумерной прямоугольной области [2, 5]. В частности, и(х,у) может быть решением краевой задачи:

д2и д2и ди е-тг^ + + а(х) —--Ь(х, у)и = f (х, у), (х, у) € Б \ Г,

дх2 ду

дх

и(х,у) = ф(х,у), (х,у) € Г,

где а(х) > а > 0, Ь(х,у) > 0, а0 = а(0), выполнены условия согласования краевых условий и функция и(х,у) является достаточно гладкой, Г — граница области Б.

Предполагаем, что для функции и(х,у) известны значения = и(х;,у,) во всех узлах некоторой прямоугольной сетки П области Б:

П = {(х;,Уj), 0 < г < М, 0 < ] < М},

и исследуем вопрос интерполяции сеточного решения в произвольную точку области Б.

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

Д;,, = {х; < х < хт, у, < у < У,+1}, К;+1 = хт - х;, Т,+1 = Уj+l - Уj.

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

иь(х,у) = ах + Ьу + сху +

Зададим функцию линейной интерполяции по х для узла у, :

и(х,у,)

К;

[(хг+1 - х)иг, + (х - хг)иг+1,,] .

(2)

;+1

Теперь определим

иь(х,у)

г?'+1

[(у,+1 - У)u(x, у,) + (у - у,)u(x, у,+1)] •

(3)

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

,'К1 \ /К1 \ 11 /-КЛ (-Ьц

иМ у,у -и -2,у) = 2 + 2ехр — -ехр

Отсюда получаем, что при е = К1

К1

иь[т,у -ит,у

К1

1 1 -1 2 + 2 ехр(-1) - ехр ( — ) « 0.077.

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

1

1

1. Экспоненциальная интерполяция

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

Yj exp(-ao£ x») + P (x», yj) = u

Yj exp(-ao£ 1x»+i) + P(x¿+b y,-) = Ui+ij, где Yj = Y (yj )• Из условия ограниченности P£ (x,y) следует

P(xi+i,yj) - P(xi,yj) = pihi+i, |pi| < Ci.

Из системы получим

= ui+i,j — ui,j — Pihi+i

Yj exp(-a0e-ixi+i) — exp(-a0e-ix»)

Учитывая соотношения

P (x,yj ) = P (x»,yj) + P2 (x — x»), |p2 | < Ci,

P (x»,yj) = Ui,j — Yj exp(—ao£-iXi), получим выражение для u(x,yj) на интервале [xi,xi+i]:

u(x,yj) = Yj [exp(—ao^-ix) — exp(—ao£-ixi)] + Uij + p2(x — x»). (4)

Пренебрегая в (4) членами порядка O(hi+i), получим интерполяционную по x формулу для интервала [xi,xi+i]:

, г ! exp(—ao£ ix) — exp(—ao£ ix») , u(x, yj) = [Ui+i,j — Ui,j ]--í-Г--í—V + u»,j • (5)

exp(—aoe-ixi+i) — exp(—aoe-ix»)

Вычитая из (5) равенство (4), получим

|и(ж, Уу ) - )| < . (6)

Итак, для каждого узла уу точность построенной интерполяционной формулы (5) порядка О(^г+1) для интервала [ж^ж^] равномерна по е. Теперь построим интерполяционную формулу для ячейки Д^-:

У) = (и) = — [(Уу+1 - У)и(ж, Уу) + (У - Уу)и(ж Уу+1)] • (7)

г?'+1

i,j

2

Лемма 1. При всех (ж, у) € Д^

К (ж, у) - и(ж,у)| < С(^¿+1 + т2+1). (8)

Доказательство. Пусть (ж, у) € • Тогда

|и(ж,у) - и^(ж, у) | < |и(ж,у) - /п^(и(ж,у))| + (и - и)| < 1 тах |т2+1 + С2^+1.

4 А;^

Учитывая ограниченность производных по у, получим утверждение леммы. □

Интерполяционная формула (7) устойчива к возмущению ui;j. Если u^(x,y) соответствует вычислениям по формуле (7) в случае возмущенных значений , когда для каждого узла (x^yj) |ui;j — | < r, то при всех (x,y) G D

K(x,y) — Uh(x,y)| < 3 r.

Таким образом, если интерпретировать, что u(x, у) — решение краевой задачи с регулярным пограничным слоем, r — максимальная по узлам сетки погрешность разностной схемы, то

|u(x, у) — Uh(x, у)| < Cmax(h + т2) + 3r.

i,j

Вычисление производной по y. Рассмотрим, как можно приближенно найти значение производной u'y(x,y) в произвольной точке области D на основе значений функции u(x,y) в узлах сетки. Для этого сначала найдем значения производной во всех узлах сетки, а затем осуществим интерполяцию на основе построенной формулы (7). Пусть Vi,j — найденное приближенное значение производной u'y (x^yj). В силу ограниченности производных функции u(x,y) по переменной y в узлах сетки на основе формул численного дифференцирования [6] можно получить приближение производной с некоторой заданной точностью. Например, в случае постоянного шага т по y для всех внутренних узлов можно определить vi;j = (ui;j+i — ui)j-i)/(2T) с точностью 0(тj). Рассмотрим вопрос интерполяции вычисленной производной в рассматриваемой двумерной области D. Воспользуемся тем, что функция v(x,y) = u'y(x,y) представима в виде, аналогичном (1):

v(x,y) = Y(y)exp(—aoe-1x) + P'y (x,y),

и поэтому можно применить формулу (7) для интерполяции значений {vi;j}. Пусть Vh(x,y) — вычисленное в соответствии с (7) значение производной в точке (x,y) G D. Тогда с учетом леммы 1

|vh(x,y) — uy(x,y)| < Cmax(hi + т2) + 3r,

i,j

где r = max |vid — u'y(xi,yj)|.

i,j

Вычисление производной по x. Покажем, что формулы численного дифференцирования, построенные на приближении функций многочленами, в случае функций с погранслойной составляющей могут приводить к значительным погрешностям. Для этого рассмотрим функцию одной переменной u(x) = exp(—e-1x). Несложно убедиться, что при е = h1

u(hi) — u(0) v --u (0)

= e-1.

Следовательно, относительная погрешность приведенной формулы численного дифференцирования не уменьшается с уменьшением шага сетки Оценивается относительная погрешность, так как производная и'(0) не ограничена при малых значениях е. Построим равномерно точную формулу для нахождения производной.

Дифференцируя построенный интерполянт (7) по х, получим приближенную формулу:

<9 ao

—uh(x,y) =--

<x е

yj+1 — yf \ . y — yj / N

-(ui+1,j — ui,j) +-(ui+1,j+1 — ui,j+1)

Tj+1 Tj+1

х-

ехр(—аое

ехр(—аое 1жг+1) — ехр(—аое

(9)

Лемма 2. При всех (ж, у) € Р

< С шах(кг + Т,- ).

Доказательство. В соответствии с (1) функция и(ж,у) представима в виде и (ж, у) = V (ж, у) + Р (ж, у), V (ж, у) = 7(у) ехр(—аое - 1ж). Пусть у) = Р(и,ж,у). Тогда у) = Р(V,ж,у) + Р(Р,ж,у),

джи^(ж,у) — дЖu(ж'У)

<

< е

5 5

_Р (^у) — -V (.„у)

+ е

дд 5жР (P'ж'У) — 5жР ^

Пусть (ж, у) € Д,". Оценим первое слагаемое в (10)

дд 5жР (V'ж'У»— ^

ао

= ао ехр(—аое ж) (7) — 7(у)| <

< ^тшах (у» |г?2+1"

4 у

Осталось оценить второе слагаемое в (10). Сначала докажем, что /Р (жг+1,у) — Р (жг,у)\ 5

Имеем

и 1 — ?Т Р (ж,у»

/Р(ж,+ 1, у) — Р(жг,у)\ 5

< С [г/+1 + ^¿+1 ].

<е

+е

\ кг+1

( Р(жг+1,у) — Р(жг,у) 5

4

5

— 5жР (ж,у)

<

к

г+1

— 5жР (ж,у)

+

/пМ 5жР(ж,у)) — 5жР(ж,у)

5

< С [к+1 + т-+1],

что доказывает (12). Несложно убедиться, что

5 j-.tr> \ т + (Р(ж¿+1,у) — Р(ж у» -р (Р,ж,у»—^ ^-—-

(10)

(11)

(12)

= (Р(ж¿+1, у) — Р(жг, у)) Фг(ж),

где

Фг(ж) =

ао ехр(—аое 1 ж)

ехр(—аое 1 ж,) — ехр(—аое 1 жг+1» кг+1

Докажем ограниченность Фг(ж). Пусть ж = ж г+1. Тогда

|Фг(жг+1) |

к

г+1

1-

ехр(т)— 1

, т = аое 1 кг+1.

е

е

е

е

е

е

е

Пусть т > 1. Тогда

0 < 1 -

т

ехр(т)— 1

< 1,

из этого следует, что |ФДх^)| < а0.

Пусть т < 1. В этом случае используем разложения по т и получим |ФДх^)| < а0/2. Пусть х = Xi. В этом случае

(Xi)|

к

i+1

т

1 — ехр(—т)

— 1

По аналогии с предыдущим случаем можно убедиться в ограниченности ф (х^|. Итак, функция ФДх) убывает и по модулю ограничена на концах интервала [х^ Хг+1]. Из этого следует ее ограниченность на всем интервале, равномерная по е и Л^+ь

Из ограниченности Р'х(х,у) следует |Р(xi+1, у) — Р(х^у)| < Ск^, следовательно,

д ^е» л т + (Р(xi+l ,У) — Р(х^У)

< С^+ь

(13)

Из неравенств (12), (13) следует оценка второго слагаемого в (10), соответствующая утверждению леммы. □

е

е

2. Случай сгущающейся сетки

Как известно, классические разностные схемы, например, монотонная схема Самарского [6] и направленных разностей, обладают свойством равномерной сходимости при решении задач с пограничным слоем, если использовать сетки, специальным образом сгущающиеся в пограничном слое. Подход, основанный на использовании сгущающейся сетки для задачи с пограничными слоями, был предложен Н.С. Бахваловым [1], позднее строились и другие сетки (см., например, [5]). Остановимся на случае сетки, предложенной Г.И. Шишкиным [4]. Известно, что ряд разностных схем на такой сетке обладает свойством равномерной сходимости в случае эллиптической задачи с регулярным пограничным или параболическим пограничным слоем. Исследуем возможности полиномиальной интерполяции для вычисления функции и производных в случае такой сетки.

Пусть для функции и(х, у) справедлива оценка производных:

&+3

"и(х,у)

<С

1 -1 ч

— ехр(—ае х) + 1

0 < г < 3, 0 < 3 < 1,

_о_ ду

-и

< С, 3 < 3.

(14)

Оценка (14) имеет место, если и(х,у) — решение эллиптической задачи с регулярным или параболическим пограничным слоем [2] у границы х = 0. В частности, функция вида (1) соответствует регулярному пограничному слою и удовлетворяет условиям (14). Пусть Пх — кусочно-равномерная сетка из [4]:

^^ х — I х.

N

xi = гк, 0 < г < —, xi

а +

г — < г < N

а = т1п{£,1 !„(^} , Н =£, Н = , N > 4. (15)

Исследуем точность формулы полиномиальной интерполяции (3) в случае сетки (15). Лемма 3. При всех (х,у) € В

|ui(x,y) - u(x,y)| < С

ln2 N П 2'

max< , Ni + max Tj

(16)

Доказательство. Сначала докажем, что при всех j и x G [x^x^]

fln2 N 1 ]

|u(x,yj) - u(x,yj)| < Cmax j , n f, (17)

где U(x,yj) соответствует формуле (2) линейной интерполяции по x при y = y-. В соответствии с формулой для погрешности линейной интерполяции на интервале [xi,xi+1]:

|u(x,yj) - u(x,yj)| < M2'x4hi+1, (18)

где М2'Х — верхняя оценка для второй производной:

д2

u(x,yj) < M2'X, x G [xi,xi+i]. (19)

дх2

Пусть г < N/2. В соответствии с (14) М2,х = С/е2. Учитывая формулу для шага сетки по х, получим

|и(х,у) - и(х,у)|< N2 1„2 N, (20)

что соответствует (17).

Пусть г > N/2. В этом случае

и(х,у^') - и(х,у^') = - и(х,у^') - х^Г1-х (иг+1,^' - ) =

Нг+1

хг+1 х

uX(s,yj) - %h j uX(s,yj)

Si+1

Из этого соотношения следует

|и(х,у^-) - и(х,у^)| < 2 у К)|

Учитывая оценки (14) и условие х^ > а, получим

Сз

|u(x,yj) — u(x,yj)| < С2 exp(-ae-1xi) [1 — exp(-ае-1Я)] < -—

Остается рассмотреть случай равномерной по x сетки, когда

L < - ln N. 2 " а

Учитывая оценку (18), несложно показать, что в этом случае справедлива оценка (20). Оценка (17) доказана.

Учитывая оценку погрешности линейной интерполяции по у и оценку (17), получим

у) - и (ж, у)| < (и(ж, у) - и(ж, у))| + (и(ж, у)) - и(ж, у)| <

< С

тах

1п2 N

N2 ,

Лемма доказана. □

Интерполяционная формула (3) устойчива к возмущениям {иг;,}, поэтому погрешность разностной схемы, используемой на сетке из [4], не накапливается при интерполяции.

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

Вычисление производной по х. Пусть (ж,у) € Дг,. Дифференцируя (3), полу-

чим

^«с(ж,у) = —

(у+1 - у) + (у - у,)

иг+1,,+1 - иг,,+1

к

г+1

к

г+1

(21)

Лемма 4. При всех (ж, у) € Д

< С

1п N

+ тах т,

Доказательство. Остановимся на случае ячейки Дг,. Справедливы неравенства

(ж,у) - ¿и(ж,у)

<

т,+1

(у,+1 - у)

Иг+1',' и (ж у )

к.

г+1

+

+ (у - у,)

иг+1,,+1 - иг,,+1

к

г+1

- иХ(ж,у,+1)

+ £

Несложно показать, что

иг+1,,' Иг>.,' 7/ (ж у )

кг

г+1

к

г+1

Жг+1 «

<

<

к

г+1

^г+1 «

В соответствии с условиями (14)

(22)

Ю*,у,-)| < С

— ехр(-а£ ¿) + 1

(23)

£

,

£

£

£

£

X ж

£

ж ж

Пусть £¿+1 < а. Учитывая соотношения (15) для шагов сетки Пж, неравенства (22) и (23), получим

ui+1,j ui,j

h

i+1

- uX(x,y)

< C In N

В случае xi > а усилим неравенство (22) с учетом (23), интегрируем и получим

ui+1,j ui,j

h

i+1

- uX(x,y)

<

C

h

i+1

[exp(—ae 1s) — exp(-ae ds+

^i+i

+ J [exp(—ae 1x) — exp(—ae ds

X

Учитывая, что при xi > а

(24)

| exp(—ae s) — exp(—ae x)| < —

N

из (24) получим

Ui+1,j Ui,j , (x y ) — ux(X,yj)

h

i+1

< C.

- N

Используя свойства линейной интерполяции по у, получим утверждение леммы. □

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

[1] Бахвалов Н.С. К оптимизации методов решения краевых задач при наличии пограничного слоя // Журн вычисл. математики и мат. физики. 1969. Т. 9, № 4. С. 841-890.

[2] Miller J.J.H., O'Riordan E., Shishkin G.I. Fitted Numerical Methods for Singular Perturbation Problems. Singapore: World Scientific, 1996. 163 p.

[3] Ильин А.М. Разностная схема для дифференциального уравнения с малым параметром при старшей производной // Матем. заметки. 1969. Т. 6, № 2. С. 237-248.

[4] Шишкин Г.И. Сеточные аппроксимации сингулярно возмущенных эллиптических и параболических уравнений. Екатеринбург: УрО РАН, 1992. 232 с.

[5] ЛисЕйкин В.Д., Петренко В.Е. Адаптивно-инвариантный метод численного решения задач с пограничными и внутренними слоями. Новосибирск: ВЦ СО АН СССР, 1989. 258 с.

[6] Самарский А.А. Теория разностных схем. М.: Наука, 1977. 656 с.

e

X

e

X

e

Поступила в редакцию 5 марта 2007 г.