3.6. Интерполяция сплайнами

Полиномиальная интерполяция не всегда дает удовлетворительные результаты при аппроксимации зависимостей. Так, например, при представлении полиномами резонансных кривых колебательных систем большая погрешность возникает на концах ("крыльях”) этих кривых. Несмотря на выполнение условий Лагранжа в узлах, интерполяционная функция может иметь значительное отклонение от аппроксимируемой кривой между узлами. При этом повышение степени интерполяционного полинома приводит не к уменьшению, а к увеличению погрешности. Возникает так называемое явление волнистости [42].

Для проведения гладких кривых через узловые значения функции чертежники используют упругую металлическую линейку, совмещая ее с узловыми точками. Математическая теория подобной аппроксимации развита за последние двадцать лет [43-44] и называется теорией сплайн-функций (от английского слова spline - рейка, линейка). Разработано и обширное программное обеспечение для практического применения сплайнов в различных областях науки и техники [45].

Рассмотрим один из наиболее распространенных вариантов

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

Функцию и будем использовать для аппроксимации зависимости заданной в узлах значениями

Если в качестве функции выбрать полином, то в соответствии с уравнением (3.19) степень полинома должна быть не выше третьей. Этот полином называют кубическим сплайном, который на интервале записывают в виде

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

В отличие от полиномиальной интерполяции, когда вся аппроксимируемая зависимость описывается одним полиномом, при сплайновой интерполяции на каждом интервале строится отдельный полином третьей степени (3.20) со своими коэффициентами.

Коэффициенты сплайнов определяются из услозий сшивания соседних сплайнов в узловых точках;

1) равенство значений сплайнов и аппроксимируемой функции в узлах - условия Лагранжа

2) непрерывность первой и второй производных от сплайнов в узлах

Кроме перечисленных условий необходимо задать условия на концах, т. е. в точках В общем случае эти условия зависят от конкретной задачи. Довольно часто используются условия свободных концов сплайнов. Если линейка не закреплена в точках вне интервала то там она описывается уравнением прямой, т.е. полиномом первой степени. Следовательно, исходя из условий (3.23) непрерывности вторых производных сплайнов на концах интервала, запишем соотношения

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

Получим алгоритм определения коэффициентов кубических сплайнов из условий (3.21)-(3.25). Условия (3.21) в узлах после подстановки сплайна принимают вид

где

Продифференцируем дважды сплайн (3.20) по переменной

Из условий непрерывности производных (3.22) и (3.23) при переходе в точке от сплайну с учетом выражений (3.28) и (3.29) получим следующие соотношения:

И, наконец, из граничных условий (3.24) и (3.25) на основании выражения для второй производной (3.29) получим, что

Соотношения (3.26), (3.27), (3.30)-(3.33) представляют собой полную систему линейных алгебраических уравнений относительно коэффициентов сплайнов Но прежде чем решать эту систему, выгодно преобразовать ее так, чтобы неизвестными была только одна группа коэффициентов

Из уравнения (3.31) коэффициенты выразим через коэффициенты

Объединяя уравнения (3.26), (3.27) с соотношением (3.34), представим коэффициенты Ь также через коэффициенты

После подстановки выражений (3.34) и (3.35) в соотношение (3.30) получим уравнение, в которое входят только неизвестные коэффициенты Для симметричности записи в полученном уравнении уменьшим значение индекса на единицу

где

При учитывая условие свободного конца сплайна, в уравнении (3.36) следует положить

Таким обрезом, уравнение вида (3.36) вместе с условиями (3.32) и (3.37) образует систему линейных алгебраических уравнений для определения коэффициентов Коэффициенты и вычисляются после нахождения - по формулам (3.34) и (3.35), коэффициенты равны значениям аппроксимируемой функции в узлах в соответствии с формулой (3.26).

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

где

Так же, как и метод Гаусса, метод прогонки разделяется на два этапа - прямой и обратный ходы. В процессе прямого хода метода прогонки вычисляют значения коэффициентов линейной связи каждого предыдущего неизвестного с последующим

Из уравнения (3.38) при с учетом граничного условия (3.32) установим связь коэффициента с коэффициентом

где и прогоночные коэффициенты,

Затем, подставляя выражение (3.40) в уравнение (3.38) при получим линейную связь коэффициента с коэффициентом

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

В процессе выполнения прямого хода метода прогонки необходимо вычислить значения всех прогоночных коэффициентов для которых получим рекуррентные соотношения. Подставим формулу связи коэффициентов

в уравнение (3.38), в результате получим

Сравнение последнего соотношения с формулой (3.41) позволяет записать рекуррентные формулы для определения прогоночных коэффициентов

Учитывая граничное условие (3.32) и соотношение

а также полагая задаем начальные коэффициенты Затем по формуле (3.42) вычислим все пар прогоночных коэффициентов

На основании соотношения

и граничного условия (3.37) получим, что

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

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

После определения всех коэффициентов - другие коэффициенты сплайнов вычисляются по формулам (3.26), (3.34) и (3.35), после чего аппроксимирующую функцию можно рассчитать с помощью соотношения (3.20) в любой точке х на интервале

Реализация рассмотренных алгоритмов интерполяции кубическими сплайнами со свободными концами осуществлена в программах 3.6.

В основном блоке программы описаны массивы для узлов X, значений аппроксимирующей функции в узлах коэффициентов сплайнов С и прогоночных коэффициентов К (строка 10). Для размещения прогоночных коэффициентов не вводится отдельного массива, чтобы сэкономить память, эти коэффициенты располагаем в массиве С. С помощью операторов в строке 20 осуществляется ввод в диалоговом режиме числа узлов начального конечного значений и шага изменения текущего значения аргумента где необходимо вычислить аппроксимирующие сплайны. Из основной программы происходит обращение к блокам формирования таблицы аппроксимируемой функции и вычисления коэффициентов сплайнов (строки 30-40). В цикле по текущему значению аргумента обращаемся к подпрограмме вычисления значений сплайнов и их производных, затем результаты выводим на дисплей (строки 50-80).

Подпрограмма формирования таблицы аппроксимируемой функции (строки 100-190) аналогична соответствующим блокам программ интерполяции другими методами, здесь осуществляется последовательный ввод с клавиатуры узлов и значений функции.

В строках 200-290 реализована подпрограмма прямого и обратного ходов метода прогонки для определения коэффициентов сплайнов. Прогоночные коэффициенты размещаются в процессе прямого хода (строки 200-250) в массиве С, при выполнении обратного хода (строки 260-270) элементы этого массива замещаются результирующими коэффициентами Арифметические операторы подпрограммы соответствуют формулам метода прогонки.

В виде отдельной подпрограммы (строки 300-390) оформлено вычисление сплайнов и производных. В зависимости от текущего значения аргумента определяется номер сплайна (переменная 1), который необходимо вычислять (строки 300-320). Затем по формулам (326), (3.34) и (335) вычисляются коэффициенты (строки 330-350). Сплайн его первая

и вторая производные вычисляются по схеме Горнера (строка 360).

В программах на языках Фортран и Паскаль вычисление коэффициентов сплайнов осуществляется с помощью подпрограмм с именем входными параметрами которых являются переменные - количество сплайнов, массивы узлов и значений аппроксимируемой функции соответственно, выходной параметр - массив С для коэффициентов сплайнов Значения аппроксимирующей функции, ее первая и вторая производные вычисляются в подпрограммах с именем со входными параметрами и выходными - совпадающими по смыслу с соответствующими переменными программы на Бейсике.

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

(см. скан)

(см. скан)

(см. скан)