Экономная интерполяция в ℝ Текст научной статьи по специальности «Математика»

УДК 518.5

n

ЭКОНОМНАЯ ИНТЕРПОЛЯЦИЯ В Ж И. А. Пахнутов

n

COST-EFFECTIVE INTERPOLATION IN Ж I. A. Pakhnutov

Рассматриваются некоторые вычислительные аспекты интерполирования в

n

пространствах многих переменных, точнее, конструирования функций f: D& ^ D& , n > 1, по конечному множеству известных значений. Основным инструментом для этой цели служат варианты напряженных, натуральных и полиномиальных сплайнов. Отмечается существенное упрощение алгоритмизации и вычислительных затрат при использовании равномерных сеток исходных данных. многомерная интерполяция, численные методы интерполяции

The paper considers some computation aspects of interpolation in multidimen-

n

sional spaces, such as constructing functions f: E ^ E , n > 1, according to finite set of known values. The computational technique is based mainly on strained, natural and polynomial splines. Massive simplification of both algorithmic and computational processes has been noted, when uniform grids of initial data are used.

multidimensional interpolation, numeric methods of interpolation

Основная задача, рассматриваемая ниже, состоит в конструировании неко-

п

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

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

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

х, у) = I к Ск фО: - Хк),

где коэффициенты с = {ск} выбираются из условий интерполяции Б(х^ х, у) = у, V 1 х = {х^ - дискретное множество точек пространства, в которых известны значения у = у(х^ исследуемой поверхности. Это - линейная задача относительно коэффициентов с, и она решается до конца (см., напр. [3]). Вопрос о погрешности в таком случае, как правило, не возникает - не с чем сравнивать.

С другой стороны, если искомая поверхность достаточно гладкая, то традиционно в окрестности каждой точки х пространства аргументов можно записать формулу Тейлора в виде

п

А(х) = Рк(х) + К|к|(х), хейс К ,

где Рк(х) = I а (х - х) - полином степени к (здесь и далее использованы

0 < 1 < к

правила мультииндексной нотации: j = {]8}, | ] | = I 1 = П]8|, х1 = П Ба =

Б 8 8

Я|а| Я|а| \ ~ к

-), Я|к|(х) = о(|х - х| ) - остаточный член. Тогда естественно возника-

дх П(дхк)ак'

ет вопрос о представлении поверхности в виде полиномов подходящей степени или полиномиальных сплайнов [4], гладко склеенных из многочленов фиксированной степени.

Если множество х {хк} (назовем его сеткой) невелико (с вычислительной точки зрения), то часто бывает достаточно ограничиться многочленами Лагранжа в форме

т

х, у) = I к ук П - , 1 (хк - х1) . к 1 Ф к ||хк - х)||

Здесь все выписывается явно и не требуется никакой вычислительной "изворот-ливости". Кроме того, сетка в данном случае может быть совершенно произвольной. Если упорядочить узлы сетки по номерам: к = 0, 1, ..., б ({ук} -соответствующее множество сеточных значений функции), то интерполяцию по узлам хт, ..., хт+г, легко можно выполнить рекурсивно (А. Аккеп, [5]):

х, у, т, г) =

ут, При Г = 0,

фГр(и,у,т,г-1), Р(и,у,т+1,г-1), " Хт) (х"1+г "2Хт)1 при г > 0,

. ^ ||Хт+г-Хт|| )

где ф(а, b, s) = (1-s)a + sb. При больших объемах исходных данных этот путь оказывается (из-за больших колебаний полинома и быстрого накопления погрешностей) совершенно непригодным даже в одномерном случае.

Некоторой альтернативой многочлену Лагранжа служит неполиномиальная линейная модель (охотно используемая в геодезии) взвешенных расстояний вида

™ ч V (||t - xk|| + 5)-в

F(t, x, У) = L k yk v ^ k L j (||t - j + 5)

при достаточно малом 5 > 0 (для избежания сингулярности) и в > 1. Очевидно, точность интерполяции и дифференциальные свойства модели зависят от выбора 5 и в (часто полагают в = 2, 5 е (0.001, 0.1)).

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

2

пластинки, то поверхность S = f(x, y), (x, y) е Q с Е , можно искать из условий минимума функционала энергии T(S) = i (S^,, + S + S) dQ и интерполяции

Q y yy

S(xk, yk) = zk, V k [6], а в общем случае рассматривать функционал

Tm(S)=1у | Da S(t) |2 dt,

q |a| = m Ш

Q с En, m > n . Решение этой вариационной задачи должно быть точным на всех многочленах Pm-1 степени не выше m-1 и может быть записано в виде

F(t) = Pm-1(t) + L Ck z(m, ||t - xk||), k

. ч [x2m"n 1п(т2), если n четное, < , n

где z(m, t) = i 2mn а множители {ck} и полином Pm-i

It , если п нечетное,

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

fF(xk) = yk, Vk, {ZkCkPm.i(xk) = 0

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

Некоторым упрощением является конструкция так называемых натуральных сплайнов [8] нечетной степени I = 2m - 1, m е N, вида

I

8£(1) = Рт-1(:) + I Ск ||: - хк|| к

при условиях Sl(xk) = ук, V к, где Рт-1 - некоторый полином степени т-1. Условие I Ск Р(хк) = 0 для всех таких полиномов должно добавляться (как и в предыдущем случае) к условиям интерполяции для корректного отыскания {ск}. Уменьшение степени полинома Рт-1 в сравнении с напряженными сплайнами несколько снижает вычислительные затраты при решении линейных систем уравнений, определяющих параметры модели.

Существенное снижение вычислительных затрат можно получить при использовании равномерных сеток данных: хк = х0 + к-И (для каждой координаты аргумента может быть выбран свой шаг сетки). Поскольку шаг по любому направлению фиксирован, можно исходные данные рассматривать на целочисленной сетке, а при окончательном вычислении учитывать размер шага. В этом случае полиномиальный сплайн удобно представить произведением одномерных сплайнов выбранной степени. Так, например, на основе кубического базисного сплайна (В-сплайна) на целочисленной сетке

1 4 з

(В(1:) > 0 при : е (0, 4)), где ф(1:)+ = тах{ф(1), 0}, можно в двумерном случае, например, взять в качестве базисной функции интерполяции Р(х, у) = В(х)-В(у), перенумеровать узлы целочисленной сетки по строкам и представить поверхность в виде функции

Б(и, V) = I.. су Р(ц - х„ у - у) (2)

с линейными параметрами {су}. При этом возникает естественная в таких случаях неопределенность: поскольку базисные сплайны определены во всем пространстве, следует и сеточные значения функции доопределить вне области интерполирования.

При использовании кубических сплайнов (1) ситуация с продолжением сеточных значений существенно упрощается: поскольку В(1:) > 0 на целочисленной сетке только в узлах : = 1, 2, 3, то достаточно по каждой переменной сеточные значения функции продолжить на один узел слева и один узел справа, например, кубическим же полиномом, полагая (в одномерном случае) у 1 = 4у - 6у + 4у - у

(на левом конце) и ук+1 = 4ук - 6ук-1 + 4ук-2 - ук-3 (на правом конце).

Пусть, например, данные для интерполирования на прямоугольной равно-

2

мерной сетке {(х;, у)) £ К : х; = х0 + гЬх, у) = у0 + 1 •Иу, 1 = 0, ..., М, 1 = 0, ..., К} представлены в виде прямоугольной (М+1)х(К+1)-матрицы. Продолжив ее указанным выше способом, получим (М+3)х(К+3)-матрицу данных. Теперь эти данные мож-

(М+3)х(№3) - 1

но организовать в виде вектора с (М+3)х(К+3) координатами X = {Ху}

(М+З) (N+3) - 1

, где у = у(1, 1) = 1(К+3) + 1 V 1, Аналогично и матрицу дан-

"у}у = 0

Г!

\х2,уЛ1 у = 0

(М+3)х(№3) - 1

ных можно записать в виде вектора У = {Уу} 0 . Тогда условия интерпо-

ляции представят систему линейных уравнений

к

I Су Р(р - 1, я - Л = Уи

V = 0

и = 0, ..., к = (М+3)(№3) -1, р =

и

N+3

, Я = ишоа N+з, 1

V

N+3

, ] = УШоа N+3, [ • ] - це-

лая часть числа. Матрица А этой системы имеет блочный трехдиагональный вид [9]:

А =

В 0.25В 0

0.25В

В 0.25В

0

0.25В В

0 0

0.25В

0 0 0 0 где матрица В имеет такую же структуру:

0 0 0

0 0 0

0 0 0

0 0.25В В ■

В=9

1 0.25 0 0 0.25 1 0.25 0 0 0.25 1 0.25

0

0

0

0

0 0 0

0 0 0

0 0 0

0 0.25 1

Линейные системы уравнений с такими матрицами легко решаются прогонкой с вычислительными затратами, пропорциональными числу уравнений. При такой регулярной структуре, как у приведенных матриц, нет необходимости в хранении последних. Таким образом, достаточно выполнить прогонку с матрицей А со столбцами матрицы исходных данных в качестве правой части, затем (относительно результата) с матрицей В, и в силу их идентичной структуры численная трудоемкость решения этих систем пропорциональна лишь объему исходных данных, так что последняя модель является наиболее экономной с вычислительной точки зрения. Гладкость результата (непрерывность вторых производных) также вполне удовлетворительна для многих практических приложений. В Я11, П > 1, трудоемкость вычислений растет пропорционально п (т.е. достаточно выполнить п процедур прогонки, сами матрицы вычислять и хранить не нужно).

Рассмотрим следующий модельный пример функции трех переменных (см. подробности в [9]):

ъ - х

^х,у,ъ) = еу2 + 1 еов(х+0.5у-0.2ъ),

заданной на равномерной сетке параллелепипеда [0, 1]х[-1, 1]х[-0.5, 1] с шагом Ьх = 1 Ц = 0.2, = 0.3, полагая х1 = гЬх, у = _рЬу - 1, ък = к-И2 - 0.5, 1, ], к = 0, 1, ... Запишем в соответствии с предыдущим

F(u, V, w) = £. кci.i.kß

u-xo

hx

i +2.

v-yo

hy

- j +2,

w-zo

hz

k +2

(3)

где теперь Р(х, у, z) = B(x)•B(y)•B(z). В результате получим функцию, графически не отличающуюся на указанном параллелепипеде от данной (с максимальной абсолютной погрешностью в заданной области 8 = 0.054 при w = 0). На рисунке представлены графики сечения функции F при w = -0.5, w = 0 и w = 0.5.

Рис. Сечения F(u,v,w) при w = -0.5, w = 0 и w = 0.5 Fig. Sections of F(u,v,w) at w = -0.5, w = 0 and w = 0.5

Если коэффициенты (с^к) в (3) вычисляются по расширенной таблице данных, то в силу симметрии В-сплайнов сдвигать аргументы ß() следует не на две, а на три единицы.

Вернемся еще раз к вычислению функции F в (2). Заметим, что В-сплайны Bn(t) степени n, формирующие базис ß(), допускают простое рекурсивное вычисление [8]

B„(t) = <

1 при 0 < t < 1, п

F - 'при n = 0, 0 иначе

1

n

[t-Bn-i(t) + (n+1 -t)Bn-1(t-1) при n > 0.

Однако эта простота полностью поглощается вычислительными затратами: рекурсивно В-сплайн вычисляется в два раза медленнее, чем при его явном представлении. В многомерном пространстве это играет далеко не последнюю роль.

СПИСОК ИСПОЛЬЗОВАННЫХ ЛИТЕРАТУРНЫХ ИСТОЧНИКОВ

1. Стейн, И. М. Сингулярные интегралы и дифференциальные свойства функций / И. М. Стейн. - Москва: МИР, 1973. - 342 с.

2. Кравченко, В. Ф. Алгебра логики, атомарные функции и вейвлеты в физических приложениях / В. Ф. Кравченко, В. Л. Рвачев. - Москва: Физматлит, 2006. - 416 с.

3. Уилкинсон, Дж. Х. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра /Дж. Х. Уилкинсон, К. Райнш. - Москва: Машиностроение, 1976. - 390 с.

4. Квасов, Б. И. Методы изогеометрической аппроксимации сплайнами / Б. И. Квасов. - Москва: Физматлит, 2006. - 360 с.

5. Hildebrand, F. B. Introduction to NUMERICAL ANALYSIS / F. B. Hildebrand. - New York: Dover Publication, 1987. - 686 p.

6. Василенко, В. А. Сплайн-функции: теория, алгоритмы, программы / В. А. Василенко. - Новосибирск: Наука (СО), 1983. - 215 с.

7. Ашкеназы, В. О. Сплайн-поверхности. Основы теории и вычислительные алгоритмы / В. О. Ашкеназы. - Тверь: ТГУ, 2003. - 82 с.

8. Бор, К. де. Практическое руководство по сплайнам / К. де Бор. - Москва: Радио и связь, 1986. - 305 с.

9. Пахнутов, И. А. Интерполяция (алгоритмический аспект) / И. А. Пахнутов. - Калининград: LAP publishing, 2017. - 121 с.

REFERENCES

1. Stein I. M. Singulyarnye integraly i differentsialnye svoystva funktsiy [Singular integrals and differential properties of functions]. Moscow, MIR, 1973, 342 p.

2. Kravchenko V. F. Algebra logiki, atomarnye funktsii I veyvlety v fizicheskih prilozheniyah [Algebra of logic, atom functions and wavelets with applications to physics]. Moscow, Fizmatlit, 2006, 416 p.

3. Wilkinson J. H. Spravochnik algorytmov na yazyke ALGOL. Lineynaya algebra [Handbook for automatic computation. Linear algebra]. Moscow, Mashinostroenie, 1976, 390 p.

4. Kvasov B. I. Metody izogeometricheckoy approksimatsii splaynami [Isogeometrical approximation methods with splines]. Moscow, Fizmatlit, 2006, 360 p.

5. Hildebrand F. B. Introduction to NUMERICAL ANALYSIS. New York, Dover Publication, 1987, 686 p.

6. Vasilenko V. A. Splayn-funktsii: teoriya, algoritmy, programmy [Spline-functions: theory, algorithms, programmes]. Novosibirsk, Nauka (SO), 1983, 215 p.

7. Ashkenazy V. O. Splayn-poverhnosti. Osnovy teorii i vychislitelnye algoritmy [Spline-surfaces. Basic theory and computational algorithms]. Tver', TGU, 2003, 82 p.

8. Bor, K. de. Prakticheskoe rukovodstvo po splaynam [A practical guide to splines]. Moscow, Radio i svyaz', 1986, 305 p.

9. Pakhnutov I. A. Interpolyatsiya (algoritmicheskiy aspect) [Interpolation in algorithmic aspect]. Kaliningrad, LAP publishing, 2017, 121 p.

ИНФОРМАЦИЯ ОБ АВТОРЕ

Пахнутов Игорь Александрович - Калининградский государственный технический университет; кандидат физико-математических наук, доцент кафедры высшей математики; E-mail: IA-Pa2010@yandex.ru

Pakhnutov Igor Alexandrovich - Kaliningrad State Technical University; cand. phys.-math. sci., Associate Professor, Department of Higher Mathematics;

E-mail: IA-Pa2010@yandex.ru