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

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

Том 12, № 2, 2007

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

В. И. ПААСОНЕН

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: paas@ict.nsc.ru

In this article we study a finite difference method of the arbitrary order of accuracy applied to solution of the interpolation problem using stretched hyperbolic splines. The method relies on the high-order compact difference schemes combined with one-sided multi-point difference approximations for the first derivatives in the smoothness conditions. This method is a direct generalization of the modification developed earlier by the author on the classical second order accuracy interpolation method. The proposed method allows both usual and the parallel realizations.

Введение

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

Данная работа является продолжением [4] и содержит описание прямого обобщения разработанного метода на схемы любого порядка аппроксимации. Метод основан на компактных аппроксимациях дифференциального уравнения в дополнительных узлах сетки и многоточечных односторонних аппроксимациях производных в условиях их гладкого сопряжения в узлах интерполяции. Метод формулируется универсально для граничных условий любой точности по аналогии с технологией [5] для расчета краевых задач в неоднородных областях. Характерная черта предлагаемого подхода — возможность параллельной реализации на основе метода [6].

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 05-01-00146-а и № 06-01-00030).

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

1. Система уравнений и ее разностная аппроксимация

Пусть (а = Хо, £1,..., Хк = Ь} — строго возрастающая совокупность значений аргумента х, которой соответствует множество значений функции {/о, /1, •••, /к}• Длину подотрезков шк = (хк-1 < х < хк} (к = 1,К) обозначим через Ик =

хк хк-1.

Интерполяционным гиперболическим сплайном £ с множеством параметров натяжения (Рк, к = 1, } будем называть (см. [3, 4]) решение з(х) дифференциальной многоточечной краевой задачи:

— - Якт = 0, Як = (Рк/Ик)2 Ух € ^, к =1,...,К; (1)

в2 8

—2 = т Ух € , к = 1,К (2)

вх2

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

з'(х — 0) = з'(х + 0), т(х — 0) = т(х + 0), х €(х1,...,хк-1}, (3)

условиями интерполяции

5(хк ) = /к, к = 0,...,К (4)

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

т(а) = ф0, т(Ь) = ф1. (5)

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

Сформулируем теперь конечно-разностную аппроксимацию многоточечной краевой задачи (1)-(5).

Построим на отрезке [а, Ь] кусочно-равномерную сетку, согласованную с границами подотрезков ^к и равномерную на каждом из них. Пусть кк = Ик/пк — шаг сетки, а п = Пк — число таких шагов на подотрезке ^к.

Пусть О = (0,...,Ж} — глобальная нумерация узлов этой сетки } и узлам

интерполяции соответствуют номера Г = (10,11,}, т. е.

£/к = хк, к = 0,...К.

Наряду с глобальной будем рассматривать также локальную нумерацию Ок = (0,..., п}, (п = пк) узлов на .

Во внутренних узлах подотрезков систему уравнений (1), (2) аппроксимируем разностной схемой

Лт, — ПЯт, = 0, г € (О \ Г); (6)

Лз, = Птг, г € (О \ Г). (7)

Здесь Л — оператор обычной трехточечной аппроксимации второй производной, числа Р, И, к принимают значения Рк, Ик, кк, если узел сетки лежит в к-м подотрезке, а постоянная П имеет вид

П = £ (20! Я '

где Ь — произвольное натуральное число. Такой выбор для коэффициента П обеспечивает аппроксимацию исходной системы (1), (2) с погрешностью высокого порядка 0(Н2Ь).

Условия интерполяции (4) и внешние краевые условия (5) задаются точно, непрерывность второй производной выполняется автоматически благодаря выбору ее в качестве искомой функции, а условия гладкости первых производных аппроксимируются непосредственно как равенство левых и правых многоточечных односторонних разностных аналогов первых производных произвольного /-го порядка аппроксимации:

1 J 1 J

1тУ2аз+ 1-У2аз=0, « = 1к, к =1,--,К - 1. (8)

Ьк 3=0 Ьк+1 3=0

При этом система условий аппроксимации первых производных односторонними разностями в (8) с порядком 3 имеет вид

J

к = , к = 0,...,3. (9)

3=0

Ее решением является набор коэффициентов

(-1)3+1

J

С3, 3 = 0, ао = -5] аз. (10)

3 3=1

Таким образом, сформулированный разностный аналог многоточечной краевой задачи представляет собой семейство разностных задач с любым порядком аппроксимации шт(3, 2Ь), произвольным ввиду произвольности натуральных чисел 3 и Ь. В частном случае Ь =1 и 3 = 2 имеем простейший метод второго порядка точности, предложенный и исследованный автором в [4]. Формально предлагаемое обобщение отличается от указанного метода тем, что в условии (8) вместо двойки в качестве верхнего предела суммирования появилось произвольное 3 (т.е. вместо трехточечных теперь используются (3 + 1)-точечные односторонние разности), а в системе уравнений (6), (7) появилась константа П, которая при частном значении Ь =1 равна единице.

Разностная задача имеет специальную ленточную структуру: в узлах интерполяции получаем равенство односторонних многоточечных аппроксимаций первых производных, использующее в качестве шаблона 23 +1 узлов, а во всех дополнительных узлах сетки — систему двух трехточечных соотношений (6), (7) для двух искомых функций. Будем предполагать, что выполнено естественное ограничение на число шагов сетки п = Пк > 3, с тем чтобы шаблоны разностных аппроксимаций в условиях гладкости не простирались за пределы ближайших подотрезков слева и справа.

2. Алгоритм численного решения разностной задачи

Предположим, что значения вторых производных в узлах интерполяции уже известны. Обозначим их Мк = т1к (к = 0,...,К). Используя локальную нумерацию узлов, сформулируем две вспомогательные задачи Дирихле на подотрезке Шк. Обобщая рассуждения работы [4], имеем

Лр - Пф = 0, ро = 1, Рп = 0; (11)

Лд - П^д = 0, до = 0, д„ = 1. (12)

Ради простоты записи здесь и в дальнейшем условимся опускать индекс к, подразумевая при этом, что речь идет о к-м подотрезке.

Тогда линейная комбинация У = Мк-1р + Мкд векторов р, д является решением однородного разностного уравнения (6) с граничными условиями У0 = Мк_1, УП = М&, т.е. вектор У на данном подотрезке совпадает с т. Аналогично, комбинация

£ = М^У + 7 + # (13)

векторов V, 17 и Ж, которые являются решениями задач Дирихле

ЛУ = Пр, Уо = К = 0, (14)

Л7 = Пд, 7о = = 0, (15)

Л# = 0, Жо = Д_ь Ж = Д, (16)

удовлетворяет разностному уравнению (7) и граничным условиям интерполяции. Поэтому значения сплайна в совпадают с компонентами векторов £, вычисленных по формуле (13), на соответствующих подотрезках.

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

1 ^ ^ 1 J ^

Г £ + Г £ ^ = 0- (17)

г=1 + г=1

Здесь знаками минус и плюс помечены величины слева и справа от текущего узла интерполяции, т.е. £_ — вектор (13) и Н_ — шаг сетки на промежутке ^, а — вектор (13) и Л,+ — шаг сетки на промежутке Подстановка (13) в (17) дает систему связей для

значений второй производной в трех последовательных узлах интерполяции:

АМк_1 + Бк+ СМк+1 = Ек, к = 1,..., К - 1, (18)

где

1 J _ 1 J

А = £ с = ^ £ ,

_ г=0 + г=0

1 1 J J

Бк = —Б_ + — Б+, Б_ = £ аг7га"_г, В + = £

+

г=0 г=0

Система (18) в соответствии с (5) замыкается граничными условиями Дирихле М0 = т(а), Мк = т(6).

Таким образом, алгоритм приближенного решения задачи интерполяции с любым порядком аппроксимации сводится к такой последовательности шагов:

1) решение на каждом подотрезке вспомогательных краевых задач (11), (12);

2) решение на каждом подотрезке вспомогательных краевых задач (14)-(16);

3) расчет второй производной в узлах интерполяции из системы (18);

4) вычисление по подотрезкам значений сплайна по явной формуле (13).

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

Заметим также, что р£ = дп-% и V = ип-% для всех г = 0,..., п. Это вытекает из того, что вектор с компонентами р£ — дп-% является решением однородного уравнения (6) с нулевыми граничными условиями, а вектор V — ип-% — решением однородного уравнения ЛУ = 0, также с нулевыми граничными условиями. Это обстоятельство несколько сокращает объем вычислений, избавляя от необходимости решать задачи (12) и (15).

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

Предполагая совпадение знаков двух слагаемых В- и В+ выражения В& (ниже это предположение будет доказано) и действуя по аналогии с [4], заключаем, что достаточным условием диагонального преобладания матрицы системы для решений и и V на каждом подотрезке является выполнение неравенств

а(V — и^ а(и + ^ > 0, и = к-. (19)

Этот критерий отличается от сформулированного в [4] только произвольностью верхнего предела суммирования 3 > 3.

3. Устойчивость алгоритма

Исследуем неравенство (19), обеспечивающее диагональное преобладание.

Заметим, что преобразование к = соответствующее линейному растяжению ко-

ординаты х, преобразует систему разностных уравнений (6), (7) к той же системе, но с константой П = 1, а именно такая система исследована в [4]. Поэтому формально решения V и и вспомогательных задач выглядят так же, как в упомянутой работе, с точностью до замены к на к. По этой причине условимся опускать подробности, касающиеся вывода точных решений вспомогательных задач.

Согласно [4]

/ ' \ *

V = к2 ОД — -Сп(р) ) , ОД = ^ £Рз. (20)

V - / ¿=1 з = 1

Рассмотрим два варианта решения в зависимости от значения параметра натяжения.

При значении параметра натяжения Q = 0 решение (20) приводится к виду

гк2(- — г)(2п — г) V£ =--.

6П

Представляя построенное решение в канонической форме полинома третьей степени от г и учитывая соотношения (9) и тождество и = ^-г, легко прямым вычислением убедиться в отрицательности обоих сомножителей в левой части неравенства (19). Конкретно

оно приводится к истинному неравенству

h2n \ I h2n

--i i--

6 / \ 2

Аналогично можно убедиться в совпадении знаков слагаемых B , B+ в выражении Bk, так как они оба отрицательны ввиду того, что в данном случае

v^ т. ^ 2n

]>>*V = - — < 0.

Диагональное преобладание в первом случае доказано.

Во втором случае Q = const > 0. Заметим, что в силу симметрии задач (11), (12) и вследствие теоремы Виета произведение корней р+ и р- соответствующего характеристического уравнения равно единице, причем в нашем случае оба корня положительны. Следовательно, один из них р = р+ > 1, тогда второй — р- = 1/р < 1. В этих обозначениях решение V вспомогательной задачи имеет вид

V = Fra (п(рга-г - рг-га) - (n - г)(рга - р-га)) , =

й2рга+1

п(р2 - 1)(р2п - 1)'

Следовательно,

Ui = Vra-i = (n(рi - р-i) - г(рга - р-п))

с тем же коэффициентом

Рассмотрим полином вида

= £ a . (21)

j=0

Подставляя в (21) значения коэффициентов (10) и проведя ряд преобразований [6], получим представление полинома в интегральной форме

г ^

Ф(*,3)= У*1 - (1 - ж) ¿ж. (22)

3 ж

1

С учетом равенств (9) для сомножителей в левой части исследуемого неравенства (19) получим выражения

У>(У* + и) = П(рП- 1) К1(р, 3,п),

£ - и) = П(ррга+ 1) К2(р, 3,п),

где

К1(р,7,п) = +Ф(Р,^),

К2(р, 7, п) = рпФ (Р, з) - Ф(р, 3) + П(рп - 1). \р У п

Так как p > 1 и > 0, все множители перед функциями Ki, K2 в обоих выражениях положительны. Опуская их, приведем неравенство (19) к эквивалентному неравенству

Ki(p,J,n) K2(p,J,n) > 0. (23)

Отрицательность обоих сомножителей в (23) при J =2 доказана V n > J и V р> 1 в работе [4]. Методом индукции распространим результат на все J > 3. С этой целью, используя интегральную форму (22) и определения функций Ki, K2, для произвольного l > J > 3 оценим знаки разностей:

Ki(p, l, n) - Ki(p, l - 1, n) = -(pn-1 + (-1)1) < 0,

K2(p, l, n) - K2(p, l - 1, n) = -(pn-1 - (-1)1) < 0,

причем оба неравенства являются строгими при l < n. Следовательно, в (23) оба сомножителя — убывающие функции второго аргумента. Учитывая сказанное выше относительно их знаков при J = 2, заключаем, что неравенство (23) справедливо для любого натурального J = 2, 3, . . . , n.

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

Для этого значения коэффициентов (10) и компоненты вектора решения V подставим в выражения для B- и B+. После сокращения на несущественные положительные множители получим выражение K2(p, J, 2n), которое отличается от второго сомножителя в неравенстве (23) только удвоенным третьим аргументом. Поскольку отрицательность K2 имеет место для любых n, утверждение доказано.

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

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

[1] ЯнЕнко Н.Н., КвАсов Б.И. Итерационный метод построения поликубических сплайн-функций // Докл. АН СССР. 1970. Т. 195. С. 1055-1057.

[2] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.

[3] КвАсов Б.И. Разностные методы построения изогеометрических сплайнов. Новосибирск: Изд-во НГУ, 2004. 48 с.

[4] Паасонен В.И. Параллельный алгоритм построения гиперболических сплайнов // Вычисл. технологии. 2006. Т. 11, № 6. С. 88-97.

[5] Paasonen V.I. Compact difference schemes for inhomogeneous boundary value problems // Russ. J. Numer. Anal. Math. Modell. 2004. Vol. 19, N 1. P. 65-81.

[6] Паасонен В.И. Сходимость параллельного алгоритма для компактных схем в неоднородных областях // Вычисл. технологии. 2005. Т. 10, № 3. С. 81-89.

Поступила в редакцию 15 ноября 2006 г., в переработанном виде — 1 декабря 2006 г.