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

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

Том 11, № 6, 2006

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

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

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

e-mail: [email protected]

A modification of the known method for construction of hyperbolic splines with stretch is proposed. In differs from the original as in the proposed method we use simpler and universal approximation for smoothness conditions at interpolation points. Therefore the proposed algorithm is more flexible and allows using of the high-order schemes and parallel methods in computations.

Введение

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

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

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

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

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

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

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

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

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

Уравнение (1) будем рассматривать также в эквивалентной форме в виде системы

Сформулируем теперь конечно-разностную аппроксимацию дифференциальной задачи

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

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

интерполяции соответствуют номера Г = {ї0, її, ... , ІК}, т. е.

(1)

где Qk = (Рк/Нк)2;

в Є С2[а, Ь],

(2)

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

в(хк) = !к, к = 0,... ,К,

(3)

и краевыми условиями

в" (а) = фо, в"(Ь) = фї.

(4)

в- = в+, в- = в"+, х Є {хї,... ,хК-1}.

(5)

Qk т = 0 Ух Є шк, к = !,..., К;

(6)

йх2

(7)

(1)-(4) или, точнее, системы (6), (7) с условиями (3)-(5).

= хк, к = 0,...К.

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

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

Лт1 — Qmi = 0, г £ (Б \ Г); (8)

Лвг = т-1, г £ (Б \ Г), (9)

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

Условия интерполяции (3) задаются точно

$1к = /к, к = 0,...,К, (10)

краевые условия (4) задаются очевидным образом как

то = фо, тм = Ф1, (11)

первое из условий сопряжения производных (5) аппроксимируется непосредственно как равенство левых и правых трехточечных односторонних разностных аналогов первых производных второго порядка аппроксимации 1 2 1 2

-~У^аэв— + ---------ав+ = 0, г = 1к, к =1,...,К — 1, (12)

Нк п кк+1 п 3=0 3=0

где

3 о 1

ао = — 2, а1 = 2 а2 = — ^,

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

Итак, разностная схема неоднородна: во внутренних узлах подотрезков имеем систему двух трехточечных соотношений, а на внутренних границах — равенство односторонних трехточечных аппроксимаций потоков. Заметим, что в методе [3] производные склеиваются иначе, а именно симметричными аппроксимациями первых и вторых производных в узлах интерполяции с использованием механизма исключения значений решения в фиктивных законтурных узлах для каждого подотрезка. Это обстоятельство принципиально ограничивает точность метода лишь вторым порядком. Предлагаемый подход отличается возможностью непосредственного обобщения на схемы более высокого порядка точности. Для этого в (12) необходимо использовать более “длинные” односторонние аппроксимации производных, а внутри подотрезков — компактную схему. Заметим, что повышение порядка позволяет получать необходимую точность интерполяции при меньшем числе расчетных узлов, благодаря чему разностный метод интерполяции может стать намного экономичнее алгебраических методов.

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

Предположим, что значения вторых производных в узлах интерполяции уже известны. Обозначим их Мк = тгк, к = 0,... ,К. Используя локальную нумерацию узлов, решим две вспомогательные задачи Дирихле на подотрезке Шк:

Лр - С^р = 0, Ро = 1, Рп = 0;

(13)

Лд - Яд = 0, до = 0, дп = 1. (14)

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

У = Ык-1'Р + Мк д (15)

векторов р, д является решением однородной задачи

ЛУ - ЯУ = 0, Уо = Мк-и Уп = Мк. (16)

Таким образом, вектор У совпадает с функцией т на соответствующем подотрезке. Так как задачи вида (13), (14) могут быть решены корректным образом (независимо и даже параллельно) в каждом подотрезке, при известных значениях второй производной в узлах интерполяции глобальный вектор решения т может быть составлен с помощью формулы (15) непосредственно из компонентов У по всем последовательным подотрезкам. Аналогично комбинация

5 = Мк-гУ + Мк и + W (17)

векторов У, и и W, являющихся решениями задач Дирихле

ЛУ = Р, Уо = Уп = 0, (18)

Ли = д, ио = ип = 0, (19)

ЛW = 0, Wо = !к-1, Wn = 1к, (20)

решаемых независимо (или параллельно) на каждом подотрезке, удовлетворяет разностному уравнению (9) и граничным условиям интерполяции (10). Поэтому компоненты глобального вектора в составляются из компонентов векторов Б на подотрезках, вычисляемых

по явной формуле (17), если только известны значения Мк второй производной в узлах ин-

терполяции. При этом устойчивость сформулированных выше трехточечных разностных задач Дирихле (13), (14), (18)—(20) сомнений не вызывает.

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

1 2 1 2

= 0. (21)

г=1 + г=і

Здесь Б- и Н- — вектор (17) и шаг сетки на Шк, а Б + и к+ — вектор (17) и шаг сетки на Шк+і. Подстановка (17) в (21) дает систему связей для значений второй производной в трех последовательных узлах интерполяции

АкМк-і + ВкМк + СкМк+і = Ек, к = 1,..., К — 1. (22)

Коэффициенты системы (22) имеют вид

Ак = ь аіУп-г, Ск = аи+,

а ее правая часть —

Система (22) в силу (11) замыкается граничными условиями

Мо = Фо, Мк = Фі-

(23)

Размерность полученной системы уравнений (22), (23) равна числу узлов интерполяции. Устойчивость алгоритма ее решения обеспечивала бы устойчивость процедуры построения сплайна в целом. Исследуем поэтому, при каких условиях матрица системы (22), (23) имеет диагональное преобладание.

Для диагонального преобладания при одинаковых знаках слагаемых в выражении Бк достаточно выполнения неравенств

Если записать такие пары условий последовательно для каждого узла интерполяции, а затем сгруппировать полученные неравенства иначе — не по узлам интерполяции, а по подотрезкам (второе из выписанных условий для данного узла и первое — для следующего узла), получим систему неравенств, сформулированную для решений и и V на одном и том же подотрезке (поэтому ниже сняты обозначения индексов — знаки “ +” и “ —”)

Ясно, что если система (22) замкнута краевыми условиями (11), то на каждом приграничном подотрезке по одному из двух условий (24), (25) следует снять. Если же на границе задана первая производная искомой функции в, то для соответствующего крайнего подотрезка необходимо требовать выполнения обоих сформулированных неравенств.

»ип_,| > |£

|£агУ+1 > |£а<и?1'

\ ^ ^ агип-г \ > \ ^ ^ агУп-г \;

| агУг\ > | ^ агиг\.

(24)

(25)

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

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

Решая формально задачи (18), (19), находим

(26)

где

г—і і г—і і

і=і з=і і=і З=і

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

В первом случае Я = 0. Подставляя решения задач (13), (14)

Подставляя эти выражения в (27), а затем полученный результат в (26), вычислим решения задач (18) и (19):

Отсюда, в частности, вытекают тождества ип-г = У для всех г = 0,... ,п, что означает тождественное совпадение неравенств (24), (25), каждое из которых приводится к эквивалентному виду

Непосредственная подстановка сюда значений аг, г = 0,1, 2, и полученных выражений для У и иг с точностью до положительного коэффициента приводит к верному неравенству п2 — 4 > 0, строгому при п > 2.

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

Покажем, что они отрицательны. Ввиду равенства ип-г = Уг , установленного выше для любого подотрезка, проверка сводится к анализу неравенства

Непосредственная подстановка в выражение для Б значений коэффициентов аг и решения У приводит к неравенству

которое является истинным Vn > 1.

Во втором случае Q = const > 0. Воспользуемся тем, что решения p и q вспомога-

найдены точно в виде линейных комбинаций степеней корней р+ и р- характеристического уравнения. Для задачи (13) с учетом граничных условий имеем

Pj = 1 — 3/п> qj = j/n в (27) и проводя упрощения, с учетом тождества

1 ■ 2 + 2 ■ 3 + ••• + (i — 1)i

(i — 1)i(i + 1) з

находим

Gi(p)

(i — 1)i (i — 1)i(i + 1)

2

6n

Vi

6n i 6n

(28)

B = ^ a:iVi < О.

тельных задач Дирихле (13), (14) для однородного уравнения в каждом слое могут быть

р—

Пусть Т — решение вспомогательной задачи

ЛТг = рг, г = 1,...,п — 1, Т0 = Тп = 0.

Явное выражение решения нетрудно найти:

Тг (р) = ^ Р ,2 (прг — грп + г — п). (30)

п(р — 1)2

В силу линейности решение задачи (18) представляет собой линейную комбинацию

Уг = С-Тг(р-) + с+Тг(р+),

где константы взяты из (29).

Заметим, что в силу симметрии задач (13), (14) и вследствие теоремы Виета произведение корней р+ и р- равно единице. Легко убедиться в том, что оба они положительны. Ввиду равноправности корней будем считать, что именно корень р = р+ превосходит единицу. Тогда второй корень, который меньше единицы, равен р- = 1/р. Учитывая это и подставляя константы из (29) и решения вида (30) в выражение для У, получим окончательно

У = К (п(рп-г — рг-п) — (п — г)(рп — р-п)) , К = п(р —'^ _ 1) • (31)

Аналогично получается решение задачи (14)

Чг = С(Р+ — P-), С = рп 1 рп

р+ — ри решение задачи (19) ( )

иг = Рп {п(р% — р г) — г(рП — р п)) (32)

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

Из (31), (32) следуют тождества ип-г = Уг для всех г = 0,...,п, поэтому и в этом

случае неравенства в системе (24), (25) совпадают, причем каждое из них приводится к

виду (28).

Вводя обозначение Ф(г) = ^ агг% и учитывая тождества аг = 0, ^ агг = 1, получим отдельно выражения для двух сомножителей в неравенстве (28):

£ аг(Уг + ГЛ) = РпП(Р"~ ^ (рпФ (1) + Ф(р) £ аг(уг — К) = Рп(рпф (р) — ф(р) + 2(рп — 1)

Так как р > 1 и коэффициент Рп положителен, множители перед скобками в обоих выражениях следует опустить. Более того, выражения в скобках можно дополнительно сократить на положительный множитель (р — 1)/2, так как в нашем случае

= (р — 1)(3 — р) ф {\ (р — 1)(1 — 3р)

р) 2р2

В результате получим эквивалентное (28) неравенство

(з— р + рП 2(1 — Зр)) ^рП 2(1 — Зр) — (з— р) + п р _ 1 ^ — 0. (33)

Очевидно, 1 — Зр < 0, кроме того, для всех п — 2 справедливо неравенство рп-2 — 1, отсюда вытекает оценка

3 — р + рп-2(1 — Зр) < 3 — р + 1 — Зр = 4(1 — р) < 0,

свидетельствующая об отрицательности выражения в первой скобке. Следовательно, условие диагонального преобладания (28) сводится к условию неположительности выражения во второй скобке, которое умножением на п приводится к полиномиальному неравенству

Зп рп-1 — п рп-2 — пр + Зп — 4(1 + р + ■ ■ ■ + рп-1) — 0. (34)

Докажем, что неравенство (34) выполняется безусловно для всех р > 1 и п — 2 и что при п > 2 выполняется соответствующее строгое неравенство.

При п = 2 в левой части неравенства (34), очевидно, тождественный нуль и неравенство оказывается равенством.

При п = З неравенство (34) превращается в истинное неравенство, причем строгое 5(р2 — 2р+ 1) > 0.

При п = 4 в левой части неравенства снова имеем положительное выражение

2(р3 — р2 — р + 1) = 2(р — 1)2(р +1) > 0.

При п = 5 требуется специальная группировка слагаемых

11 р4 — 9р3 — 4р2 — 9рр + 11 = 2(р — 1)2(р — 1)2 + 9(р3 — 1)(р — 1) > 0,

которая получается на основе разбиения коэффициентов 11 = 9 + 2, 4 = 2 + 2.

При п = 6 поступаем аналогично:

7р5 — 5р4 — 2р3 — 2р2 — 5рр + 7 = 2(р3 — 1)2(р2 — 1) + 5(р4 — 1)(р — 1) > 0,

представляя 7 = 5 + 2 и группируя отдельно члены с коэффициентами, равными по модулю.

В общем случае при п > 5 левая часть неравенства (34), записанного в симметричном виде

(Зп — 4) рп-1 — (п + 4) рп-2 — 4( рп-3 + рп-4 +-+ р3 + р2) — (п + 4)р + (Зп — 4) — 0,

также сводится к сумме положительных слагаемых специальной группировкой слагаемых. Для этого представим первый и последний одночлены в виде суммы двух слагаемых на основе представления коэффициента Зп — 4 = (п + 4) + (2п — 8). Тогда сумма четырех одночленов с коэффициентами ±(п + 4) легко разлагается на положительные множители (п + 4)(рП 1 — рП 2 — р + 1) = (п + 4)(рП 2 — 1)(р — 1).

Для оценки суммы оставшихся слагаемых рассмотрим отдельно случаи четного и нечетного п.

При четном n коэффициент (2n — 8) при старшей степени и число (2n — 8), оставшееся от свободного члена, представим как сумму четверок, а центральный полином с многоточием, имеющий перед скобкой коэффициент —4, разобьем на две группы одночленов посередине так, чтобы в первую группу попали члены со старшими степенями, а во вторую — с младшими. Каждый одночлен первой группы объединим с соответствующим слагаемым степени n — 1 с коэффициентом 4, а каждый одночлен второй группы — с числом 4 как составной частью свободного члена. Тогда каждая сумма симметричным образом сгруппированных пар слагаемых положительна:

4(pn-1 — pn-k) — 4(pk-1 — 1) = 4(pn-k — 1)(pk-1 — 1) > 0, k = 3,4,..., n/2.

При нечетном n отличие в доказательстве состоит лишь в том, что коэффициент (2n—8) представляется в виде суммы многих четверок и одной двойки, а в центральном выражении с многоточием оказывается нечетное число слагаемых и при разбиении их на две группы остается центральный член с коэффициентом —4. С двумя группами поступаем так же, как в случае четного n, а центральный член разбиваем на два равных слагаемых, группируя одно с 2pn-1, а другое — с числом 2, суммируя затем их и разлагая эту сумму на положительные множители.

Для завершения исследования диагонального преобладания необходимо показать отрицательность слагаемых коэффициента Bk в разностном уравнении (22). Подставим, как это было сделано в предельном случае Q = 0, в выражение для B значения коэффициентов и решения V. После сокращения на несущественные положительные множители в результате получим неравенство

2 p2n — 1

р2п 2(1 — 3p) — (3 — р) + 1“ < °

n p — 1

Замена 2n на n преобразует левую часть данного неравенства в точности к выражению во второй скобке неравенства (33), отрицательному по доказанному выше при n > 2. Ввиду характера замены наше неравенство выполняется Vn > 1.

Таким образом, устойчивость метода доказана при любом n > 2 и любых параметрах натяжения. Заметим, что Vn > 2 в условиях диагонального преобладания неравенства строгие. Рассмотрение предельного случая Q = то опускается ввиду его тривиальности.

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

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

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

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

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

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

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