Отмечая расхождения в трактовке en.Wikipedia и ру.Рувики, даётся комментарий к истории появления сплайн-функций. Рассматривается на конкретных примерах, как методы расчёта изгиба балок «предвосхитили» некоторые из достижений математической теории сплайнов.
Вики – Рувики
С середины января 2024 года пользователям интернета стала доступна отечественная онлайн-энциклопедия «Рувики». По словам разработчиков, будет сделан акцент на достоверность, качество и нейтральность подачи информации. Для чего вводится понятие «выверенная» и «рецензированная» статья.
Любопытно было сравнить статьи из en.Wikipedia и первые статьи ру.Рувики по некоторым неоднозначно трактуемым научно-техническим разделам. Интересный подход использован в статье из ру.Рувики о изобретении радио: «Американцы считают, что радио изобрели Д.Хьюз и Т.Эдисон, немцы − Г.Герц, русские – А.С.Попов, западноевропейцы – Г.Маркони». Но наряду с этим приведён обширный список дат с описаниями вклада конкретного автора в изобретение радио. Другими словами, использован подход: «Наши, конечно, лучшие, а вот кто первый – судите сами … по фактам». В других разделах тоже есть, что обсуждать, но остановлюсь на вопросе, о котором имею некоторое представление. Речь пойдёт о сплайнах и сплайн-функциях. Сплайн-функции получили своё название от гибких деревянных или металлических реек, которые использовались чертёжниками при проведении плавных линий через заданные точки при разработке (на плазе) чертежей обводов судов, самолётов и автомобильных кузовов.
В разделе «История создания» статьи в en.Wikipedia.org о сплайнах не понравилось безапелляционное утверждение:
«Общепринято, что первой математической ссылкой на сплайны является статья Шёнберга 1946 года…» Вряд ли стоит отрицать, что И.Шёнберг впервые в 1946 году ссылался на функции (которые он стал называть сплайнами). Однако, математическая запись этих функций, их возникновение, использование и связанные с этим идеи появились задолго до этой даты. В статье [13] Ю.С.Волкова и Ю.Н.Субботина «50 лет задаче Шёнберга о сходимости сплайн-интерполяции» отмечено, что сплайны, без употребления самого термина, применялись и раньше. В 1938 году их использовали В.Кваде и Л.Коллац. Ссылка на это есть в книге [12, 1981 г.]. А в книге Дж.Альберга, Э.Нильсона, Дж.Уолша «Теория сплайнов и ее приложения» [7, 1967 г.] (или Из-во «Мир» Москва, 1972) сказано: «существует очень тесная связь между теорией сплайнов и расчётом балок» и далее − в методах расчёта «предвосхищены некоторые из недавних достижений теории сплайнов». Покажем, на конкретных примерах, как возникало это «предвосхищение». Постараемся наметить в хронологическом порядке сетку математических ссылок на появление и применение функций (которые впоследствии стали называть сплайнами).
Балки и сплайны в технической механике
Математическое представление сплайн-функций идёт из теории расчёта изгиба балок, основы которого заложили [~1750 г.] Л.Эйлер и Д.Бернулли. Эта теория стала востребованной перед строительством Эйфелевой башни в конце 19 века.
Изгиб балок постоянного поперечного сечения при малых прогибах и малых углах поворота описывается дифференциальным уравнением второго порядка (1).
Решение этого уравнения для случая, когда на балку перпендикулярно к ней действует несколько, разнесённых по длине, поперечных сил, является кусочно-полиномиальная функция третьей степени y(x). Она представляет форму упругой лини изогнутой балки (её среднего волокна). В дальнейшем эту функцию стали называть кубическими сплайнами. В настоящее время используют несколько видов математического представления таких сплайнов. Одно из «современных» представлений кусочно-полиномиальной функции (сплайна) можно найти в книге [1] профессора Ивана Григорьевича Бубнова «Строительная механика корабля», изданной в 1912 году. Эта функция задаёт форму упругой линии нагруженной балки. Для случая двухопорной балки, показанной на рисунке 1 (вместе с эпюрой изгибающих моментов), запись функции даётся формулой (2).
В записи (2) прерыватели Бубнова указывают, что следующее за ним выражение добавляется лишь при x> xi, а иначе оно равно нулю. При EI =1 и большом числе n точек xi приложения сил структура аналитической зависимости y(x) упругой линии имеет вид (3).
Это и есть представление сплайна с прерывателями Бубнова. В изданиях книги А.А.Попова «Сопротивление материалов» [14] функции (3) с прерывателями Бубнова использованы в расчётах упругих линий балок.
Aвгуст Фёппль [3, 1944 г.] ввёл обозначение с подстрочным знаком «+», которое имеет такой же смысл, что и прерыватель Бубнова. Это обозначение («преобразование Фёппля») задаётся выражением (4).
При таком обозначении прерывателя Бубнова сплайн (3) при EI =1 записывается в виде формулы (5).
В «современной» записи [10, 1976 г.] сплайна И.Г.Бубнова прерыватель Бубнова обозначают подстрочным знаком "+", а определяют его функцией max(см. формулы (6)-(7))
Стоит отметить, что В. Лангнер в статье [9, 1970 г.], ссылаясь [6, 1961 г.] на связанный с судостроением отчет Ф.Тайльхаймера и В.Стракветера, называет выражение (5) − «кусочно-полиномиальная функция в форме записи (представления) Тайльхаймера» − своего соотечественника. Однако, когда публиковалась книга И.Г.Бубнова [1], содержащая запись (2) с прерывателями Бубнова, Тайльхаймеру было два года. А в 1933 году он эмигрировал из Германии в США.
Интересно также отметить, что в упоминаемой выше статье [9] используется для гибкой рейки чертёжника термин из судостроения − «штракен» (Straken, Straklatte – штрак-рейка=страк-рейка), а не сплайн.
Математика изгиба балки и задача сглаживания данных
Рассмотрим как «математика» изгиба балки под воздействием сил, передаваемых через упругие элементы (пружины – на рисунке 2), используется для сглаживания точечных данных.
Проиллюстрируем это на примере расчёта балки (см. рисунок 3) с заделкой на левом конце, где y(x0)=0; y'(x0)=0, и подвижной опорой справа, где y(x2)=0; y''(x2)=0. Единственная пружина тянет балку вверх силой P к точке (x1, у1).
Раскроем статическую неопределённость и определим силу R2 из условия y(x2)=0 отсутствия перемещения на правом конце балки: R2δ22 – Δ2P = 0, где δ22 – перемещение от единичной силы, действующей в сечении x2 по направлению R2. Δ2P – перемещение в сечении x2 от действия силы P. Перемножая эпюры моментов (см. рисунок 3.в) по правилу Верещагина или вычисляя интегралы, получим силу R2 (см. формулы (8)-(9)).
Из условий равновесия для сил и моментов находим R0, M0, M1,M2 (см. формулы (10)-(13)).
Подстановка в формулу (5) даёт выражения (14) и (15) для упругой линии изгиба балки.
Если принять жёсткость c пружины бесконечно большой, то упругая линия пройдёт (см. рисунок 3.б) через точку (x1,у1)=(0.5,0.5) закрепления пружины. Это соответствует интерполяции, а сила P при этом определяется из выражений (16).
Если снизить жёсткость пружины, то «упругая рейка» (балка), стремясь распрямиться, сместится вниз. Например, при какой-то жёсткости c* ордината y1* станет равной 0.3. Эта жёсткость c* определяется соотношениями (17). Нижняя кривая на рисунке (3.б) соответствует изогнутой «гибкой рейке» (балке), нагруженной силой P от растянутой пружины. При другой трактовке задачи можно принять, что допустимо отклонение (погрешность) в 0.2 от полученной из эксперимента точки (0.5, 0.5), и определяется наиболее гладкая линия (функция y*(x) ), которая учитывает данные эксперимента.
Рассмотрение изгиба балки под воздействием сил, передаваемых через упругие элементы (пружины), приводит к задаче, которую можно также решать минимизацией суммы потенциальной энергии изгиба и энергии деформации пружин. Эта сумма выражена формулой (18) или формулой (19) – в форме с безразмерным параметром сглаживания λ.
Записав необходимые условия экстремума, можно получить систему линейных алгебраических уравнений для определения неизвестных коэффициентов функции y*(x) в представлении (5). При безразмерном параметре λ, определяемом (при c* из (17)) формулой (20), получаются значения (21) и кривая (см. рисунок 3б), рассчитанные выше.
В случаях с несколькими пружинами, воздействующими на гибкую рейку, параметр можно подбирать (с учетом их количества) ориентируясь на среднее квадратичное отклонение упругой линии от опорных точек. Такой подбор выполняется решением нелинейного уравнения с относительно неизвестного параметра λ.
Е.Уиттекер [2, 1923г.] предложил при аппроксимации (со сглаживанием) точечные данных использовать подобный подход. Там вторая производная в выражении энергии изгиба была представлена второй разделённой разностью. Эту же идею использовал И.Шёнберг [4,1946 г.], но там уже в выражении энергии изгиба рассматривалась вторая производная.
К.Райнш опубликовал [8, 1967г.] алгоритм (на языке Алгол-60) аппроксимации, регулирующий подбор (итерационным решением нелинейного уравнения) степени сглаживания кубическим сплайном по заданной величине отклонений от заданных точек (с учётом их весов).
Минимум потенциальной энергии изгибной деформации балки математики обосновали (Дж. Холлидей [5, 1957 г.]) и «часто называют свойством минимальной кривизны» [7] сплайна. Заметим, что сплайны, хотя и связаны с теорией изгиба балок, но не придерживаются допущений для балок. Для гибкой сплайн-рейки кривизна и энергия деформации, вообще говоря, имеют более сложные выражения (22).
От сил к точкам – интерполяция
Добавим количество сил, действующих на балку с опорами на концах (см. рисунок 4), и рассмотрим расчёты, порождающие задачу интерполяции.
Сначала распределим силы равномерно x0=0, x1=1, x2=2, x3=3, x4=4 и примем P4= P0,
P3= P1=P, то есть рассмотрим симметричный случай (см. рисунок 4). А силами на концах обеспечим нулевые углы поворота (как при заделках).
При этом на левом конце балки: y(x0) = y0 = 0, y'(x0) = 0, y''(x0) = 0. В таком случае выражение (5) для упругой линии балки примет вид (23). Запишем в виде системы (24) равенство нулю суммы сил и условие равенства нулю угла поворота (y'(x0) = 0) на правом конце балки.
Потребуем, чтобы упругая линия имела равную единице ординату в точке x2, то есть y(x2) =y2 =1. Из этого условия определим в (25) значение силы P. А зная величину силы P=6, можно определить другие (см. (26)) недостающие значения и перемещение y1=y(x1).
Если связать перемещения yi = y(xi) не с силами, а с коэффициентами ai, то для данного случая (с «нулевыми условиями» на левом конце балки) эти коэффициенты можно вычислять по формуле (27).
Если же рассматривать зависимости для изгибающих моментов (29), то можно получить их связь с координатами для общего случая.
Проинтегрировав дважды от xi до xi+1 соответствующее выражение для изгибающего момента, получим зависимость (30). А чтобы избавиться от производной y'i проинтегрируем дважды это же выражение для изгибающего момента от xi до xi-1 и получим зависимость (31).
После почленного вычитания из (30) выражения (31) получим известное соотношение (32), связывающее ординаты (перемещения) с моментами в точках приложения сил. По моментам, выражающих вторые производные, и ординатам в узловых точках можно определить сегменты кубических полиномов, а также вычислить коэффициенты an по формулам (33).
Вообще, рассмотренная в этом примере упругая линия y(x) представляет локальный сплайн B3(x) третьей степени на равномерной сетке узлов.
Рассмотрим далее общий случай расположения узлов (точек приложения сил). В этом случае система условий для определения коэффициентов
aT = [a0,…, a4] в выражении упругой линии балки примет вид (34).
Локальный кубический сплайн B3(x) с коэффициентами, которые определяются из системы (34), имеет максимальное значение равное единице в средней точке x2. Такой локальный сплайн удобно использовать как базис для представления интерполяционного сплайна s(x). Но возможно и иное нормирование сплайна.
При взгляде на матрицу в (34) напрашивается желание записать в её последнюю строку соответствующие длины отрезков hij в четвёртой степени. Такую запись даёт условие (см. формулы (35)), связанное с интегралом от функции y(x).
Это условие с интегралом означает, что площадь под кривой y(x) (до оси x) будет равна одной четвёртой. И естественно, что система уравнений (35) с такой «красивой» матрицей тоже должна иметь «красивое» решение. Записав алгебраическое дополнение к единице в правом верхнем углу матрицы (35) и раскрыв его как определитель Вандермонда, получаем выражение (36).
Последовательно заменяя столбцы в матрице вектором правых частей и раскрывая определители как определитель Вандермонда, получаем выражение (37) для коэффициентов ai сплайна B3(x). А сам сплайн B3(x) с такими коэффициентами можно теперь записывать по-разному (см. (38)-(40)) и в том числе определить его, как в [11], через разделённые разности (40).
Кусочно-линейные и базисные функции
В ряду первых применений сплайнов первой степени стоит метод Эйлера (метод ломаных) для численного решения дифференциальных уравнений первого порядка. Само решение представляет не что иное как сплайн 1-ой степени. Безразрывные эпюры моментов также являют собой примеры непрерывных кусочно-линейных функций. А эпюру на рисунке 1 можно рассматривать ещё и как пример базисной функции 1-ой степени – базиса, который удобно применять для представления более сложных составных эпюр. По существу базис в неявном виде используется при переходе, например, от эпюры перерезывающих сил к эпюре моментов (см. рисунок 6) и позволяет вычислять изгибающий момент в произвольной точке. Базисные функции можно использовать и при переходе от момента к углу поворота, и далее – к прогибу, то есть к упругой линии балки.
Рассмотрим использование базисных функций нулевой степени для представления эпюры перерезывающих сил Q(x) и перехода от него к эпюре моментов M(x). Сделаем это на примере балки, показанной на рисунке 6.а.
Базис 0-ой степени (1-го порядка) образует функция, заданная на сетке узлов t0,…,tn формулой (41). Её график показан на рисунке 5.
Индекс j в обозначении функции привязывает её к узлу tj сетки, от которого вправо следуют ненулевые значения функции. А «степень» 0 у символов N0 и B0 функции указывает, что это функция нулевой степени.
Эпюра перерезывающих сил Q(x) с рисунка 6.б может быть представлена формулой (42), где θ(u) − функция Хевисайда, равная 1 при u≥0 и 0 при u<0.
Через базисные функции Q(x) записывается в виде выражения (43), поскольку справедливо соотношение баланса сил (44), обращающее в ноль слагаемое при N20(x).
Показанная на рисунке 6.в эпюра изгибающих моментов M(x) подобна базисной функции
1-ой степени (2-го порядка – по числу отрезков-носителей). Рассмотрим её получение на основании известной формулы Д.И.Журавского: dM(x)/dx=Q(x).
Интегрируя выражение (42), получаем выражение (45) для момента.
Момент равен нулю при x ≤ t0 и x ≥ t2. Для интервала t0 < x < t2 последнее слагаемое в (45) можно отбросить, а для подинтервалов выражения для изгибающего момента даются формулами (46) и (47).
При выводе формулы (47) учитывалось соотношение (48), выражающее равенство моментов сил Va и Vb относительно точки t1. .
Полезно отметить, что этот, вполне очевидный результат в (47), можно получить, если интегрировать Q(x) по x в направлении от t2 к t1 (см. формулу (49)).
Используя выражение (45) и (47), а также базисные функции 0-ой степени, получаем выражение (50) для изгибающего момента.
Вернёмся к стандартным (хотя и менее наглядным, чем использованным на рисунках 6.а-6.в) обозначениям: (51) и (52).
Напомним, что для определения P0, P1, P2 используется такая же как (35) система уравнений (53) и формулы (54) для вычисления, аналогичные формулам (37). А множитель перед интегралом должен на единицу превосходить степень базисной функции.
Коэффициенты в (52) для показанных на рисунке 6г базисных функций (эпюр изгибающих моментов) вычисляются по формулам (54).
Для значений x, принадлежащих отрезку [t1, t2], получаем соотношения (55) и (56). Из них видно, что на отрезке [t1, t2] сумма этих базисных функций не равна единице.
Но если нормировать B01(x) и B11(x), как записано в формулах (57), то сумма получаемых нормализованных сплайнов N01(x) и N11(x) окажется постоянной и равной единице (см. (58)). Тогда для отрезка [t1, t2] сумму изгибающих моментов M(x)=M0(x)+M1(x) можно выразить через нормализованные сплайны. А приняв M0(t1)=M1=y1, M1(t2)=M2=y2, записать её в виде формулы (58).
Отмеченным важным свойством функций Ni1(x) и объясняется необходимость умножения на размер носителя базисного сплайна при нормировании.
Подстановка в (50) выражения для Va =P0 и Vb =P2 даёт (59).
А с учётом нормирования получаем известное в теории сплайнов соотношение (60) для нормализованных базисных сплайнов. Оно в общем виде записано в (61), и определяет базисные сплайны, показанные на рисунке 7.
«Обратное преобразование Фёппля» и «abs-интерполяция»
Рассмотрев примеры расчётов балок и потренировавшись в «предвосхищении» достижений теории сплайнов, постараемся тоже что-нибудь открыть. Запишем сплайн Ɲi1(x), используя соотношения (45) и (54), в виде формулы (62), где индекс i приписан к узлу под вершиной сплайна.
Применим к формуле (62) «обратное преобразование Фёппля», то есть перейдём обратно от скобок с подстрочным плюсом (см. формулу (4)) к выражению (64) с абсолютными величинами (abs-функциями).
Перепишем теперь с учетом (64) выражение (58) в виде интерполяционной формулы (65).
Обойдем неудобство в назначении двух дополнительных узлов приемлемым для интерполирования способом − сохраним за границами отрезка интерполирования [x1, xn] граничные значения, то есть y1 слева и yn справа. Для этого расположим левый и правый дополнительные узлы сетки соответственно в –∞ и в +∞. В таком случае в формуле следует принять допущение (66) и окончательно использовать для кусочно-линейной интерполяции формулу (67).
Для примера применения формулы (67) интерполируем точки данных, приведенные в таблице на рисунке 8. Полученная зависимость s(x) представлена формулой (68). Отметим, что вне отрезка интерполяции [t1, tn] при x ≤ t1 функция s(x)=0, а при x ≥ tn функция s(x)=1. То есть за пределами этого отрезка функция сохраняет значениям граничных ординат.
Итак, состоялось «открытие»! Была открыта … старая тетрадка, и из неё переписана формула с абсолютными величинами. Она позволяет на [t1, tn] задавать по точкам ломаную (для набивки в один клик сразу двух модульных скобок!). Конечно, такое «открытие» достойно нового термина: «abs-интерполяция».
Но вернёмся от шуточек к основной теме и …
Подведём итоги
Было показано, что математические записи функций, которые впоследствии стали называть сплайнами, появились задолго до статьи И.Шёнберга 1946 года. И что методы расчета изгиба балок позволили «предвосхитить» первые достижения в теории сплайнов. Идеи этих методов перешли в теорию сплайнов.
А в связи с цитатой из en.Wikipedia.org: «Общепринято, что первой математической ссылкой на сплайны является статья Шёнберга 1946 года, в которой, вероятно, впервые слово "сплайн" используется …», и указанным за ней внушительным перечнем дальнейших разработок, вспоминаются слова Менделя Маранца: «Что такое идея? – Яйцо. Все зависит от того, кто на нём сидит».
Более взвешенным представляется текст из Ruviki.ru: «Теория интерполяции сплайнами и сам термин сплайн ведут свой отсчёт со статьи Иссака Шёнберга 1946 года. Особенно интенсивное её развитие произошло в 50-70 годы». Эта констатация не отрицает предшествующих результатов и здравого смысла, подсказывающего, что в случае с яйцом также важна генетика. А список того, что было сделано и заложено раньше, (до статьи И.Шёнберга) также важен. И, наверное, не лишне знать: «Кто же первым записал сплайн-функцию?»
Литература
1. Бубнов И.Г. Строительная механика корабля. Часть I. CПб, Тип. мор. м-ва, 1912.
2. Whittaker E.T. On new method of graduation. Proc. Edinburgh Math. Soc. 41 (1923), 63- 75; p.236
3. Föppl A. Vorlesung über Technische Mechanik, Bd. 3. Berlin, 1944.
4. Schoenberg I.J. Contributions to problem of approximation of equidictant data by analitic functions, Quart. Appl. Math., 4 (1946) 45-99, 112-141.
5. Holladay J.C. Smoothest curve approximation. Math.Tables Aids Comput., 11 (1957), 233-243.
6. Theilheimer F. und Strakweather W. The Fairing of Ship Lines on a High-Speed-Computer. David Tylor Model Basin, Report 1474, 1961.
7. Alberg J.H., Nilson E.N., Walsh J.L. The Theory of Splines and Their Applications, Academic Press, New York, 1967.
8. Reinsch C.H. Smoothing by spline functions. Numer. Math. 10, 1967, 177-183.
9. Langner W. Die Lösung Strakproblem bei empirische Funktionen mittels stückweiser kubischer Polynome. Elektron. Rechenanl. 12 (1970), H. 5, S. 262-269.
10. Стечкин С.Б., Субботин Ю.Н. Сплайны в вычислительной математике, М. «Наука», 1976. 248 с.
11. De Boor C. A practical guide to splines. New York–Heidelberg–Berlin: Springer Verlag. 1978. 326 p.
12. Schumaker L.L. Spline function: basic theory. New York: Wiley, 1981, 553 p.
13. Волков Ю.С., Субботин Ю.Н. 50 лет задаче Шёнберга о сходимости сплайн-интерполяции, Труды ин-та механики и математики УрО РАН, 2014, т.20, №1
14. Попов А.А. Сопротивление материалов. Машгиз, М. 1958.
