CC BY
8Вестник Сыктывкарского университета. Сер. 1.Вын.2.1996
УДК 519.652
Дискретные периодические сплайны и решение
задачи о бесконечной цилиндрической оболочке 1
А. В. Певный
Изучаются дискретные периодические сплайны на равномерной сетке и их В-сплайновое представление. Решается задача интерполяции . Предлагаемый метод, в основе которого лежит дискретное преобразование Фурье, позволяет избежать решения систем уравненийи приводит к явным выражениям для коэффициентов сплайнов. * Метод применяется для решения задачи о бесконечной цилиндрической оболочке. . .
Введение
В статье В.А.Желудева [1] развивается своеобразное сплайн-операционное исчисление (СОЙ), основанное на дискретном преобразовании Фурье(ДПФ) и на представлении периодического сплайна S(x) в виде линейной комбинации В-сплайнбв. В данной работе исследуются дискретные сплайны S(j),j £ Z, имеющие период N. Начало изучению дискретных сплайнов положили работы [2],[3]. Систематическое изложение теории дискретных сплайнов см.
в И,[5].
В основе дискретного варианта СОИ лежит представление дискретного сплайна S(j) в виде линейной комбинации дискретных В-сплайнов
m—1
S{j) = ^2ciQr(j - ji), 1=0
где jt = In - узлы крупной сетки на Z. Далее от последовательности коэффициентов {q} переходим к ее ДПФ, и задача переносится в
'Работа выполнена при поддержке Российского фонда фундаментальных исследований (проект 96-01-00269).
© Певный A.B., 199*6.
пространство Фурье-образов. Использование ДПФ позволяет для сокращения объема вычислений применять алгоритмы быстрого преобразования Фурье. Необходимая информация из области дискретного гармонического анализа имеется, например, в [6].[7].
1. Рекуррентное определение дискретных периодических В-сплайнов
1. Дискретные сплайны рассматриваем на множестве целых чисел Z. Зафиксируем натуральные г, п, /V, такие, что Лг > 2m, п > 2.
Определим дискретные В-сплайны Q\ (_у),..., Qr(j) {j G Z) следующим образом. Положим
\n-h j € 0 :п— 1,
- Qi(j) = < 0, jen:N-n, ' (1.1)
I j-(N-n), j G N — n + 1 : Дг — 1.
Определенную на 0 : Д* - 1 функцию Q\ продолжим на всю целочисленную ось Z с периодом N (так, чтобы выполнялось условие Qi(j + N) =г Çi(j) Vj € Z). Функция Qi п])едставляет собой совокупность " домиков1' с вершинами в точках О, ±N. ±2N, — Носитель Q\ есть
\J[kN-n + l : kN + rï- 1].
kez ,
Далее используем рекуррентное определение
Qk = Qy* Qk-i, к = 2,... ,т, (1.2)
или ■ . ' ■
.V-1
= £ Qi {p)Qk~]{j - р), J € Z, А: = 2,... , г.
Лемма 1.1. В-сплайн Qr имеет период N и обладает следующими свойствами: ,
Qr(—j) = Qr(j) Для любого целого j; (1.3)
Qr{j) > 0 при kN - r(n - 1) < j < kN + r(n - 1) . . и Qr(j) = 0 при остальных j;
r(n-l)
£ = (1.5)
аи/<-1)) = 1. (1.6)
Доказательство п1)оводится индукцией по г. При г = 1 свойства (1.3)-(1.6) легко проверяются. Пусть имеет период N и обладает свойствами (1.3)-(1.6). Тогда также имеет период-.N. Имеем по индуктивному предположению
N-1
ЯЛ-}) = =
р=о
N-1
= Е ял-^яг-м -») = ялп у* € г.
■ и=0 •
При 3 £ Ъ справедлива формула
71-1
ЛЛз)= Е (1.7)
р--» + 1
При ] 6 —г(п — 1) : г(п — 1) в сумме (1.7) все слагаемые неотрицательны и хотя бы одно положительно, поэтому 01г(]) > 0-.
При г(п — 1) < < ДГ — г(п — 1) будем иметь {] — р) — 0 при всех р Е — п + 1 : п — 1, поэтому (},■{]) = 0. С учетом периодичности выполнено (1.4).
Из (1.7) и индуктивного предположения следует, что
г(п—1) „_1 (г-1)(п-1)
Е ялз) = Е ЯЛР) Е:
Наконец, при j = г(п — 1) в сумме (1.7) ровно одно положительное слагаемое: (}г(г(п — 1)) = (}\{п — 1)ф,._1(г(п — 1) — (п — 1}) = 1 -1. Лемма доказана.
2. Пусть теперь N делится нацело на п:
■ - - • - . " N — гпп.' .
В силу сделанных предположений п > 2, гп > 2г > 2.
Дискретным А^-периодическнм сплайном 5 будем называть линейную комбинацию В-сплайнов
7П.-1
= О'ег), (18)
/=о
где j| = 1п, I £ 0 : т — 1, называются узлами сплайна. Такое представление дискретных сплайнов будет использоваться на протяжении всей статьи. „
Остановимся на вычислении для всех ] 6 0 : Лг — 1. Пусть 3 € [)У,.7*+])• Тогда, считая последовательность {с/} периодически продолженной влево и вправо с периодом т (с/+т == С( V/), можно записать равенство
к+г
Ё (1-9)
1=к-г+\
Всего 2г умножений с/ на С<).
Модено показать, что сплайн на каждом целочисленном промежутке ^ : совпадает с некоторым полиномом Р* степени не выше 2г — 1. Это можно использовать для быстрого вычисления сплайна 5. Поясним идею алгоритма на примере г — 2. В этом случае на : функция 5 является полиномом З-ей степени
где 2 € 0 : п, = 5(л---(- 1) - 5(.д-)- Используя (1.9), можно
вычислить конечные разности
50 5(л-) = а1 * + аО * с* + а1 *
= а2* ск-1 - п*ск + аЗ*сш,
52 := Д250У) = (п - 1)^-1 + (3 - 2п)с* + (п - 3)с*+1 + сш,
53 = -с^г-НЗс* - ЗсЛ+, + с*+2, - . . • .
где константы аО = (2п3 + п)/3, а1 == (п3 — п)/С, а2 = —п(п — 1)/2, аЗ == п(п + 1)/2 можно вычислить заранее. Значения полинома
Р{г) = 50 + 51 • 2 + ^52 • г(г - 1) + ^53 • ф - 1)(г - 2)
можно вычислить в точках 2 € 0 : п — 1, не используя операцию умножения. Воспользуемся рекуррентными соотношениями
Р(г + 1) = Р(г) + 0(г + 1),
д(г + 1)=<?(г) + Г>(2 + 1),
D(z + l) = D{z) + SZ,
с начальными условиями
Р(0) = SO, Q{0) = 51 - 52 + S3, D{0) = 52 - 53 - 53.
Здесь используется разность R — 52 — 53, для которой справедлива формула R n(ck-1 — 2ск + ct+i).
Окончательно вычисление сплайна S(j) на промежутке [jk,ik+\) реализуется следующим фрагментом программы
скО с[к - 1]; ск := ф]; ск\ := с[к + 1]; ск2 := с[к + 2];
50 := al * (скО + ск1) + «0 * ск;
51 := а2 * скО - п * ск + аЗ * cfcl; А:— ск - сА:1; В := А — скО;
S3:=A + A + B + ck2\R:=-n*{B + cky, Р ;= 50; Q := 51 - 7?; D := R - S3; for 2 := 1 to n — 1 do
begin D :— D + 53; Q := Q + D; P := P + <5 end
Переменная P последовательно принимает значения P(0),P(1), ..., P(n — 1). Всего в приведенном фрагменте 6 умножений и « Зп сложений, а для вычисления значений сплайна на т отрезках потребуется бгп умножений и ~ Ътп = 3/V сложений.
2. ДПФ и интерполяция дискретными сплайнами
Будем использовать аппарат дискретного преобразования Фурье (ДПФ). Пусть сих = е2;г,/Л'. ДПФ вектора а = {а,}^-1 имеет вид
j=o
Последовательность {/4/.} является iV-периодической. Введем норму вектора
Fk(a) := Ак := £
ДПФ обладает свойствами
dj — — Y] AkUift, j = 0,... ,N — 1 (формула обращения)
/V
N-1 /V-1
у=о ' ¿=0
1
||а||2 = |А,ь|2 (равенство Парсеваля). (2.1)
" ¿=о
Здесь черта обозначает комплексное сопряжение.
Найдем ДПФ от В-сплайна фг(/), введенного в п.1. Сверточное свойство (1.2) позволяет легко вычислить ДПФ от <2Г. Действительно,
- П(0г-1)П«?1) = ....= [П(с?1)]г, А- = о,1...../V-1.
Осталось вычислить ДПФ от С)):
Д72-1. .
(2-2)
Учтем, что на —N/2 : N/2 — 1 сплайн — это исходный В-сплайн (1.1). Формула (2.2) принимает-вид
1 1 / ¡¡.— 1 71 — 1
Fk(Q{) n + J2(n - j)i^NkJ+4ч -»+ 2E^n -^cos
2nkj N
j=г , j=l
(здесь важно, что N¡2 — I > и — 1. Последнее выполняется ввиду условия N = ттг > 2п). Используем тождество Фейера
п—1 . /0 «-,
.. . /sin nx/z\z п +2} {п - j) cosjх = I —------t~ .
' \ SllU'/Z )
j-1 '
При ж ='27rk/N получаем
-\smirk/N / ибо- N — mn. При к = О
Для ДПФ от Qr будем использовать обозначение urk ~ Fk(Qr). В силу (2.1)
и,-к =
п
2 г
(т\(ктт/т)/зт(ктг/М))2г при к = 1,..., А' - 1.
■при к = О,
Но ¡¡формуле обращения 1
Жз) =
N
п
2г ( (§\п1ж/тл^г ^
Е/ ЬШ 4 /1 / III \
ш
N
, зеж.
(2 3)
Эта формула получена в [9].
Из (2.3) получаем значения В-сплайна (¿г{з) в узлах д. =:кп:
1
т— 1
■ШМ) = - п^-Ч^Л^/т^
т . А—'
^=1
где
н-1
п ^ \5т(5-га 4- Лл/М
ч 2г
) > 3 = 0,...,т- 1.
Таким образом, последовательность {^-(л-)}/-!^1 является обратным ДПФ (длины т) от положительной последовательности
= тг2'^1. 7) = Ау ( яп 2', ^ = 1,.'.., т - 1.
(2-4)
Положительность чисел Т) позволяет легко решить интерполяционную задачу для дискретных сплайнов степени 2г — 1:
771-1
£(.?'*) = Е С^г{3к - 31) = ¿ь Ж = о,..., т - 1.
(2.5)
1=0
Здесь числа даны, требуется найти коэффициенты с/.
Отметим, что матрица <3 с элементами С1к1 — Яг{3к ~ 31) симметрична и положительно определена. Симметричность следует из (1.2), а положительная определенность - из соотношения
т— 1 т—1 ^ т— 1
((¿с, с) := У ЯнСкС, = V - V =
■к.,1=0
т— 1 то—1
к=0 193
для любого ненулевого вектора с = {с*}^1-
Решать систему = .г не нужно. Надо сделать ДПФ над обеими частями (2.5), Придем к эквивалентной системе уравнений
т—1
£ $(3к)^ = 4 3 = 0,..., т - 1. (2.6)
¿•=о -
Но из (2.5) видно, что {б'Сд)} есть циклическая свертка последовательностей {с/} и {Я 1-й к)}- По теореме о свертке система (2.6) принимает вид .,'..'.
"'ОД-= 4 (2.7)
где
¡и--1 /=о
После перехода к Фурье-образам получили -простейшую (диагональную) систему (2.7), которая ввиду положительности чисел 7) легко решается:
= ,'; .42-8)
Для нахождения неизвестных коэффициентов с/ сплайна 5 осталось сделать обратное ДПФ. Если учесть, что числа Ту можно вычислить заранее, то можно сказать, что решение интерполяционной задачи (2.5) сводится к. выполнению двух ДПФ длины т.
Полученный результат можно сформулировать в матричной форме (это частный случай общего результата о факторизации пра-воциркулянтных матриц [8]).
Теорема 2.1. Для матрицы — {ЯЛЗк — здУи-о справедливы равенства
где Т = (Над{7}),.. .,Тш-\}, Р = >х"7-1о " матрица ДПФ. Замечание.
Пользуясь (2.8), выведем формулу для суммы В-еплайнов. Рассмотрим интерполяционную задачу (2.5), в которой все ¿к = 1- Тогда по формуле (2.8) получим
п ~ 7 !Т _ / т/п2г~', ] = 0,
С помощью обратного ДПФ найдем коэффициенты сплайна: Ск = 1/п2'"-1, к = 0,..., т — 1. Итак, сплайн 8(]) = 1 записывается в виде
»i-i
t=0
Отсюда получаем полезное тождество »i-i '
£ От0' - *п) ^ n2-1 Vj 6 Z. (2.9)
к=.0
3. Решение задачи о бесконечной цилиндрической оболочке
Метод ДПФ может применяться для решения периодических задач математической физики. Проиллюстрируем на примере одной задачи теории упругости.
Рассмотрим круговую цилиндрическую оболочку радиуса Я, толщины с/, имеющую бесконечную длину, причем й/Я« 1. Предположим, что нормальная р\{<р) и касательная р-^'-р) компоненты внешней нагрузки не зависят от координаты х вдоль оболочки и таковы, что перемещения <¿(9), и>{<р) также не зависят от х и удовлетворяют системе уравнений ' .
Г «' - 2а, > + а; + а,[>> - 2ш"] = \ «'' + ш' - 2а0ч,'" = -ЛаРаС-У»),
Где ч;; \
1 2ч ^ а0
= тгт; -4а = ^1(1 - ^ ), а0 = ттгБ9' а1 ~
ЯЛ' ' и " 12В?" 1 I-!/2'
Здесь 0 < // < 1/2—коэффициент Пуассона, Е > 0—модуль упругости.
Требуется найти 27г-периодическое решение системы (3.1). Если ( и, ги)— решение, то (и + С, ш)—также решение. Константа С определяется из какого- либо дополнительного условия. Кроме того, интеграл по отрезку [0,2тт] от левой части второго уравнения (3.1) равен нулю, поэтому естественно предположить, что
/•2тг
7 р2{ч>)Лф - 0- (3.2)
Пусть N — mit. Рассмотрим схему, в которой вычисляются значения функций « {у), <'(у>) на мелкой сетке у>у — jh, j = О----, iV — 1,
где h = 27г/Лг, но при этом делаются ДПФ длины то.
Аппроксимируем u(çj) дискретным УУ-периодическим сплайном Si(j), a w(ipj) —сплайном S2(j)- Сплайны S\, S? разложим по В-сплайн&м:
m—1 m—1
Зд) = Е c'QrU - ji), s2(j) = diQr(3-ji)\
1=0 1=0
где ji = In, r - фиксированное натуральное число.
Введем h = nh\ = 2тт/т- -шаг к])упной сетки. Для гладкой функции / обозначим Д = f{kh) значения на крупной сеткё. Рассмотрим операторы!)1, D1. D'\ D4, аппроксимирующие 1-ю, 2-ю, З^ю и 4-ю производные соответственно с погрешностью порядка h'2:
D'h = (Д+1 - A-i)/(2/>). D*fk = (Д+1 - 2/, + /,. ,)//;.2,
- [Д+2--. fk—ï - 2( Д+1 - Д-г)]/(27г*3), J^'A - (Л+24Д+г + 6Д - 4Д_, + Д_2)А4.
Коэффициенты ci,di находятся из системы линейных уравнений
= А\Р1(кк), к ~ 0,...., т — 1,
(п) + Я'ЗД*) - 2а0Я3ЗД*) = ~А2р2(кк), к = 0,...,т~ 1.
Возьмем ДПФ длины гп от обеих частей этих равенств. Последовательность {3[(}1г)} есть циклическая свертка последовательностей с = {с,} и {д,.0'/)}, поэтому ДПФ {5,(Д)} есть {СкТк}^0\ где Ск = /ч-(с)- Аналогично ДПФ {ЗД*)} есть {ПкТк}™^, где = Рк{й). Запишем теперь обратные ДПФ
771—1 , 771—1
ад) =4Еад) = ~Y1
m *—' m '
/=о /=о
. m—1 - m—1
pi (**) = - E Pii^h) =-Y,
m « rn. —*
m ' ГП ,
/=() 1=0
Применим операторы I)1, £>2, £>3,£?4 к (3.4). Имеем т • 1 '_7'~1 • .
где я/ == 5ш(2/7г/тп). Аналогично,
т-1
1=0
где VI — (2 эт ~)2,
1=0 4 7
- . »1-1 2
т А—' /г . ..
; ■ • /=о
Подставив в (3.3), получим для каждого /. = О,..., т — 1 систему двух уравнений с неизвестными С/, Д:
а,с,+/?;!>, = л^Ри/п, ;;
7/С, + ¿/А = -Л,/ут,,
где
>7 Л , 2а ^а /и,2 2г>,\
7/
и/
7 г / 2а0г;/\
'= гд V "Ж")'
При I = 0 имеем а-о = 7о = ¿о = 0 и ввиду условия (3.2) можно считать, что
т— 1
Рю :=£р2(М)=0.
¿=о
Поэтому £>о = А\Р\ъ!Тъ, а Со можно взять произвольным. Если определители Д/ = а/6/ — /?/7/ отличны от нуля при / = 1,... ,т — 1, то системы имеют единственные решения
С, = (AiPuS, + A-tPuM/biT,. Dt = -(AiPuii + A2P.2la,)/
Делая обратное ДПФ длины m, находим коэффициенты q, di сплайнов Si(ji), S2U). При этом легко проследить, как входит произвольная постоянная С'о в сплайн S\(7):
1 1
<, = -Со + с-;, где с/ CW*', 1 = 0.,..., т - 1;
т т £—J
Лг — 1
Отсюда
1 III— I Ml—1
■ U) = -СоУ Qr(j - ji) + V c:\il-U - 3l) -
III ' *—'
/=0 /=0
1 m~ 1
^-cw^ + Vcmj-jih
in
/=0
ибо, как установлено в п.2, первая сумма равна к2' -1.
Приведем результаты расчетов при Л = 50, г/ 1. В этом случае параметр «о = 3.3 • 10~° весьма мал, что осложняет вычисления. В качестве точного решения возьмем и(<р)' = —5cos2<рч »'(у) = sin2<р, подставим в систему (3.1) и определим pi(<p) и pii'-p)- Для иллюстрации характера сходимости Si, 5*2 к u,w при т —> оо приведем таблицу отклонений
ш I'Sl Ь') ~ w( ) I 1О') ~ ш ('-Pi) I ''г 1
256 1.1 ■ ю-' 2.2 ■ ■ Ю-1
512 5.6 ■ ю-2 1.1 • ■Ю-1
1024 1.8 ■ ю-2 3,6' Ю-2
2048 5.2 • 10~3 1.0 ■ Ю-'2
4096 1.6 • 10"3 2.8 • кг3
8192 4.3 : 10-1 7.6 • 10~4
Время счета на РС АТ 386 при т = 8192 не превосходило одной минуты.
Благодарю В.Л.Никитенкова за постановку задачи о бесконечной цилиндрической оболочке и ВА.Кирушева га разработку модуля БНА по дискретному гармоническому анализу.
Литература ,
1. Желудев В,А. Периодические сплайны и быстрое преобразование Фурье// Ж. вычисл. машем, и матем. фаз. 1992. Т.32. т. С.1296-1309.
2. Whit taker Е.Т. On a new method of graduation//Proc. Edinb. Math. Soc. 1923. V.41. P.63-75.
3. Mangasarian O.L., Schumaker L.L. Discrete splines via mathematical programming//SI AM J. on Control. 1971. V.9. №2. P.174-183.
4. Schumaker LX. Spline functions: basic theory. New York: Wiley, 1981.
5. Вер M.F. Дискретные сплайны и задача восстановления дискретных данных: Дне... канд. физ.-мат. наук. Л.: Изд-во Ле-нннгр. ун-та, 1990. 140 с.
6. Макклеллан Дж.Х., Рейдер Ч.М. Применение теории чисел в цифровой обработке сигналов. М.: Радио и связь, 1983.
7. Влейхут Р. Быстрые алгоритмы цифровой обработки сигналов. М.: Мир. 1989.
8. Воеводин В.В., Тыртышников Б.Б. Вычислительные процессы с теплицевыми матрицами. М.: Наука, 1987.
9. Belir M.G., Malozemov V.N., Pevnyi А.В. A discrete version of spline-operational calculus/Оптимизация конечно-элементных аппроксимаций. Тезисы докладов межд, конф. СПб, 1995.
10. Вер М.Г., Малоземов В.Н. Наилучшие формулы для приближенного вычисления ДПФ//Ж. вычисл. матем. и матем. физ. 1992-T.32.JYsU .С.1709-1719.
11. Вер М.Г., Малоземов В.Н. О восстановлении дискретных периодических данных // Вестник Ленингр. ун-та. Сер.1. 1990. Вып.З. С.8-13.
Summary
Pevnyi A. B. Discrete periodic splines and solution of the problem concerning infinite cylindrical shell
The definition of the discrete N-periodic B-spline is given. Every discrete N-periodic spline is the linear combination of shifts of B-spline. The problem of interpolation by discrete N-perioclic splines is solved. A method for the solution of the problem concerning the infinite cylindrical shell is presented. The method uses the Discrete Fourier Transformation. It permits to get rid of the solution of the systems of algebraic equations and gives explicit formulae for coefficients of splines.
Сыктывкарский университет. Поступила 18.08.90
