Ошибка интерполяции кубическими сплайнами

Материал из MachineLearning.

Перейти к: навигация, поиск

Введение

Постановка математической задачи

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

( 1 )

Основная цель интерполяции — получить быстрый (экономичный) алгоритм вычисления значений для значений , не содержащихся в таблице данных.

Интерполируюшие функции , как правило строятся в виде линейных комбинаций некоторых элементарных функций:

где — фиксированный линейно независимые функции, — не определенные пока коэффициенты.

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

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

.

Тогда по заданным однозначно определяются коэффициенты .

Изложение метода

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

Кроме того, на границе при и ставятся условия

( 2 )

Будем искать кубический полином в виде

( 3 )

Из условия имеем

( 4 )

Вычислим производные:

и потребуем их непрерывности при :

( 5 )

Общее число неизвестных коэффициентов, очевидно, равно , число уравнений (4) и (5) равно . Недостающие два уравнения получаем из условия (2) при и :

Выражение из (5) , подставляя это выражение в (4) и исключая , получим

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

( 6 )

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

( 7 )

Условие эквивалентно условию и уравнению . Разностное уравнение (6) с условиями (7) можно решить методом прогонки, представив в виде системы линейных алгебраических уравнений вида , где вектор соответствует вектору , вектор поэлементно равен правой части уравнения (6), а матрица имеет следующий вид:

где и .

Метод прогонки

Метод прогонки, основан на предположении, что искомые неизвестные связаны рекуррентным соотношением:

( 8 )

Используя это соотношение, выразим и через и подставим в i-e уравнение:

,

где — правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать

Отсюда следует:

Из первого уравнения получим:

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

Пример: интерполирование неизвестной функции

Построим интерполянту для для функции , заданной следующим образом:

1	1.0002
2	1.0341
3	0.6
4	0.40105
5	0.1
6	0.23975

В результате интерполяции были рассчитаны следующие коэффициенты интерполянты:

				Отрезок
1,0002	-0,140113846	0,440979231	-0,266965385
1,0341	-0,291901538	-0,359916923	0,217718462
0,6	-0,22553	0,293238462	-0,266658462
0,40105	-0,100328462	-0,506736923	0,306015385
0,1	-0,134456154	0,411309231	-0,137103077

Ошибка интерполяции

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

от шага h, где .

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

причем константа в этом неравенстве является наилучшей из возможных

Пример: интерполяция синуса

Постром интерполянту функции на отрезке , взяв равномерно отстоящие узлы с шагом 0.5 и шагом 0.25, и сравним полученные результаты.

	Ошибка интерполяции	Оценка ошибки	Иллюстрация
	0.429685	3.(3)	Результат интерполяции sin(4x) с шагом 0.5
	0.005167	0.208(3)	Результат интерполяции sin(4x) с шагом 0.25

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

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

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

• А.А.Самарский, А.В.Гулин.  Численные методы М.: Наука, 1989.
• А.А.Самарский.  Введение в численные методы М.: Наука, 1982.
• Костомаров Д.П., Фаворский А.П.  Вводные лекции по численным методам
• Н.Н.Калиткин.  Численные методы М.: Наука, 1978.

См. также

• Интерполяция каноническим полиномом
• Тригонометрическая интерполяция
• Рациональная интерполяция
• Проблема выбора узлов для интерполяции
• Практикум ММП ВМК, 4й курс, осень 2008

Время на прочтение
11 мин

Количество просмотров 79K

Доброго времени, Хабр!

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

Начну с небольшой вводной. Будучи студентом 4-го, на тот момент, курса бакалавриата, я изучал курс «Компьютерная графика». Много там было разных интересных (и не очень) заданий, но одно прямо особо запало мне в душу: интерполяция кубическими сплайнами с заданными первыми производными на концах интервала. Пользователь должен был задавать значения первых производных, а программа — считать и выводить на экран интерполяционную кривую. Особенность и основная сложность задания заключена в том, что задаются именно первые производные, а не вторые, как в классической постановке сплайн-интерполяции.
Как я ее решал, и к чему оно в итоге пришло, я как раз и изложу в этой статье. И да, если по описанию задачи вы не поняли ни в чем ее смысл, ни в чем сложность, не переживайте, все это я также постараюсь раскрыть. Итак, поехали.

А, нет, погодите один момент. Вот вам два числовых ряда:
a) 2, 4, 6, 8, ?
b) 1, 3, ?, 7, 9

Какие числа должны стоять на месте вопросов и почему? Вы действительно уверены в своем ответе?

Интерполяция

Интерполяция, интерполирование (от лат. inter-polis – «разглаженный, подновлённый, обновлённый; преобразованный») – в вычислительной математике способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений. (с) Википедия

Поясню на примерах. Существуют задачи, когда нам требуется узнать, условно, «закон распределения» (взял в кавычки, так как это, вообще говоря, термин из другой области математики) некого параметра по нескольким известным его значениям. Чаще всего речь идет об изменении некого параметра во времени: координаты движущегося тела, температуры объекта, колебания курса валюты, etc. При этом в силу каких-либо обстоятельств у нас не было возможности наблюдать за этим параметром непрерывно, мы могли узнавать его значения лишь в какие-то отдельные моменты времени. Исходными данными в таком случае у нас является множество точек вида value(time), а целью задачи – восстановить кривую, проходящую через эти точки и непрерывно описывающую изменение этого параметра.

Следует понимать, что невозможность постоянного наблюдения за соответствующим параметром – это обычно какого-либо рода технологическое ограничение. С развитием техники таких ситуаций становится все меньше и меньше. Из современных задач такого плана – траектория движения, например, марсохода. Поддерживать непрерывный сеанс связи (пока что) все еще не представляется возможным, а контролировать его перемещение и рисовать красивые траектории хочется. Получается, что конкретные координаты можно узнать только в те моменты, когда связь все-таки налажена, а траекторию целиком приходится восстанавливать по полученным таким образом время от времени точкам.
Другой вариант применения интерполяции. Некоторые современные телевизоры показывают изображение с частотой обновления картинки до >=1000Гц (хотя это все еще запредельные значения). Большинство телевизоров так не умеет, но даже так многие отображают картинку на частоте 100Гц — такая величина уже вполне себе классика. А если верить википедии, то в кинематографе «частота 24 кадра в секунду является общемировым стандартом». Для того, чтобы превратить 24 кадра в секунду исходного видеопотока в 100 кадров в секунду результата, телевизор использует интерполяцию. А именно какие-нибудь алгоритмы в стиле «взять два соседних кадра 1 и 2, посчитать разницу между ними и сформировать из нее 3 дополнительных кадра, которые надо впихнуть между теми двумя изначальными» -> получаются кадры 1, 1_1, 1_2, 1_3, 2

Для дальнейших рассуждений возьмем более простой пример. Представим себе, например, лабораторную работу по географии в каком-нибудь 6-ом классе (кстати, у меня когда-то и правда была такая). Необходимо каждые 3 часа измерять температуру воздуха и записывать данные, а потом сдать учителю график изменения температуры от времени суток. Допустим, по результатам измерений у нас получилась вот такая табличка (данные придуманы случайным образом и никак не претендуют на какую-либо правдоподобность):

Отобразим полученные данные на графике:

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

Количество условий и степень интерполирующего полинома

Можем ли мы вообще гарантировать, что такая функция, которая соединяет все заданные точки, вообще существует?

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

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

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

Возвращаясь к нашей задаче с температурой – в ней мы определили 6 точек, значит, для того, чтобы провести полином единственным образом, он должен быть 5-ой степени

Интерполирующий полином тогда будет выглядеть так:

А сейчас следует сделать важное замечание и пояснить, что я имел ввиду под «условием». Полином можно задать не только координатами точек, через которые он проходит, условиями могут быть любые параметры этого полинома. В простейшем случае это действительно координаты точек. Но в качестве условия можно взять, например, первую производную этого полинома в какой-либо из точек. Вторую производную. Третью производную. В общем, любую возможную производную в любой из точек, в которой этот полином существует. Поясню на примере:
Прямую можно задать однозначно, как я уже говорил, двумя точками:

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

С полиномами более высоких степеней можно использовать и более сложные условия (вторая производная, третья производная, etc.), и каждый такой параметр будет идти в общий счет количества условий, которые однозначным образом определят этот полином. Чтобы не быть голословным, вот еще пример:

Пусть нам заданы такие три условия:

Условий три, значит, мы хотим получить полином второй степени:

Подставляем

Считаем первую производную и считаем

Считаем вторую производную и считаем

Отсюда получаем, что наш полином выглядит так:

Интерполяция кубическими сплайнами

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

Принципиальное отличие идеи сплайн-интерполяции от интерполяции полиномом состоит в том, что полином один, а сплайн состоит из нескольких полиномов, а именно их количество равно количеству инервалов, внутри которых мы производим интерполяцию. В примере с нашей температурой воздуха, в которой у нас определено 6 точек, у нас будет 5 интервалов – соответственно, у нас будут 5 полиномов, каждый на своем интервале.

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

Возвращаясь к рассказанному в предыдущем пункте, для того, чтобы однозначно задать один полином 3-ей степени, необходимо 4 условия. В этой задаче у нас 5 полиномов, то есть, чтобы задать их все, нам нужно суммарно 5∙4=20 условий. И вот как они получаются:

1) Первый полином определен на первой и второй точках – это два условия. Второй полином определен на второй и третьей точках – еще два условия. Третий полином, четвертый, пятый – каждый из них определен на 2-х точках – суммарно это дает 10 условий.

2) Для каждой промежуточной точки из множества (а это 4 точки с временами 12:00, 15:00, 18:00, 21:00) должно выполняться условие, что первые и вторые производные для левого и правого полиномов должны совпадать. Формульно:

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

3) Остаются два условия, которые пока еще не определены. Это так называемые «граничные условия», от задания которых и зависит, какой именно сплайн получится. Обычно задают вторые производные на концах интервала равными 0:

Если сделать так, то мы получим так называемый «естественный сплайн». Для вычисления таких сплайнов написано уже огромное количество библиотек, бери и используй любую.

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

И вот мы подошли к условию моей задачи. Преподаватель придумал такое задание, что задаваться должны первые производные и на левом и правом концах интервала, а программа должна считать интерполирующую кривую. А для такого требования готовых алгоритмов я не нашел…
Я, разумеется, не стану описывать весь твой «творческий» путь от момента, когда я услышал задание, до того, как я его сдал. Расскажу лишь саму идею и покажу ее реализацию.

Сложность задания состоит в том, что, задавая первые производные на концах интервала, да, мы задаем этот сплайн. Теоретически. А вот посчитать его на практике – задача довольно сложная и совершенно неочевидная (желающие могут посмотреть код нахождения естественного сплайна на Вики – ru.wikipedia.org/wiki/Кубический_сплайн – и попробовать его понять хотя бы). Разумеется, я совершенно не хотел провести кучу времени, закопавшись в матан и пытаясь вывести нужные мне формулы. Я хотел более простое и элегантное решение. И я его нашел.
Рассмотрим наш сплайн и возьмем первый из его интервалов. На этом интервале уже заданы 3 условия:

— задается пользователем

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

— ничем не обоснованное предположение

Таким образом, зная эти 4 условия, мы полностью определяем этот полином. Зная все параметры этого полинома, мы можем вычислить значения первой и второй производных на второй точке, и поскольку они совпадают со значениями первой и второй производной для полинома на втором интервале, это приводит к тому, что мы также определяем и второй полином:

— вычисляется из

— вычисляется из

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

А поскольку такой сплайн, который бы удовлетворял заданным условиям, гарантированно существует, и существует в единственном экземпляре, это значит, что мы можем рассмотреть разность:

и попытаться найти такое значение , при котором обращалась бы в 0 – и это будет тем самым правильным значением , которое строит искомый пользователем сплайн:

Самое замечательное в моей идее то, что эта зависимость оказалась линейной (вне зависимости от количества точек, через которые мы проводим сплайн. Этот факт доказан теоретическими подсчетами), а значит можно случайным образом взять любые два начальные значения и , посчитать и , и сразу же посчитать то самое верное значение, которое построит нам искомый сплайн:

Итого, мы гарантированно находим искомый сплайн за 3 прогонки таких вычислений.

Немного кода и скриншотов программы

class CPoint
{
    public int X { get; }
    public int Y { get; }

    public double Df { get; set; }
    public double Ddf { get; set; }

    public CPoint(int x, int y)
    {
        X = x;
        Y = y;
    }
}

class CSplineSubinterval
{
    public double A { get; }
    public double B { get; }
    public double C { get; }
    public double D { get; }

    private readonly CPoint _p1;
    private readonly CPoint _p2;

    public CSplineSubinterval(CPoint p1, CPoint p2, double df, double ddf)
    {
        _p1 = p1;
        _p2 = p2;

        B = ddf;
        C = df;
        D = p1.Y;
        A = (_p2.Y - B * Math.Pow(_p2.X - _p1.X, 2) - C * (_p2.X - _p1.X) - D) / Math.Pow(_p2.X - _p1.X, 3);
    }

    public double F(int x)
    {
        return A * Math.Pow(x - _p1.X, 3) + B * Math.Pow(x - _p1.X, 2) + C * (x - _p1.X) + D;
    }

    public double Df(int x)
    {
        return 3 * A * Math.Pow(x - _p1.X, 2) + 2 * B * (x - _p1.X) + C;
    }

    public double Ddf(int x)
    {
        return 6 * A * (x - _p1.X) + 2 * B;
    }
}

class CSpline
{
    private readonly CPoint[] _points;
    private readonly CSplineSubinterval[] _splines;

    public double Df1
    {
        get { return _points[0].Df; }
        set { _points[0].Df = value; }
    }
    public double Ddf1
    {
        get { return _points[0].Ddf; }
        set { _points[0].Ddf = value; }
    }
    public double Dfn
    {
        get { return _points[_points.Length - 1].Df; }
        set { _points[_points.Length - 1].Df = value; }
    }
    public double Ddfn
    {
        get { return _points[_points.Length - 1].Ddf; }
        set { _points[_points.Length - 1].Ddf = value; }
    }

    public CSpline(CPoint[] points)
    {
        _points = points;
        _splines = new CSplineSubinterval[points.Length - 1];
    }

    public void GenerateSplines()
    {
        const double x1 = 0;
        var y1 = BuildSplines(x1);
        const double x2 = 10;
        var y2 = BuildSplines(x2);

        _points[0].Ddf = -y1 * (x2 - x1) / (y2 - y1);

        BuildSplines(_points[0].Ddf);

        _points[_points.Length - 1].Ddf = _splines[_splines.Length - 1].Ddf(_points[_points.Length - 1].X);
    }

    private double BuildSplines(double ddf1)
    {
        double df = _points[0].Df, ddf = ddf1;
        for (var i = 0; i < _splines.Length; i++)
        {
            _splines[i] = new CSplineSubinterval(_points[i], _points[i + 1], df, ddf);

            df = _splines[i].Df(_points[i + 1].X);
            ddf = _splines[i].Ddf(_points[i + 1].X);

            if (i < _splines.Length - 1)
            {
                _points[i + 1].Df = df;
                _points[i + 1].Ddf = ddf;
            }
        }
        return df - Dfn;
    }
}

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

Достоинства и недостатки алгоритма

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

Навскидку можно сказать, что сложность алгоритма — O(N), так как, как я уже говорил, вне зависимости от количества точек, достаточно двух прогонов вычислений, чтобы получить правильное значение второй производной на левом конце интервала, и еще одного, чтобы построить сплайн.

Впрочем, если кому-то захочется покопаться в коде и провести какой-нибудь более подробный анализ этого алгоритма, я буду только рад. Напишите мне разве что о результатах, мне было бы интересно.

Так а в чем провинились тесты IQ?

В самом начале статьи я написал два числовых ряда и попросил их продолжить. Это довольно частый вопрос во всяких IQ тестах. В принципе, вопрос как вопрос, но если копнуть чуть глубже, окажется, что он довольно бредовый, потому что при некотором желании можно доказать, что «правильного» ответа на него не имеется.

Рассмотрим для начала ряд «2, 4, 6, 8, ?»
Представим себе этот числовой ряд как множество пар значений :

, где в качестве мы берем само число, а в качестве – порядковый номер этого числа. Какое значение должно быть на месте ?

Мысль, к которой я стараюсь плавно подвести – это то, что мы можем подставить абсолютно любое значение. Ведь что по факту проверяют такие задачи? Способность человека найти некое правило, которое связывает все имеющиеся числа, и по этому правилу вывести следующее число в последовательности. Говоря научным языком, здесь стоит задача экстраполяции (задача интерполяции состоит в том, чтобы найти кривую, проходящую через все точки внутри некоторого интервала, а задача экстраполяции – продолжить эту кривую за пределы интервала, «предсказав» таким образом поведение кривой в дальнейшем). Так вот, экстраполяция не имеет однозначного решения. Вообще. Никогда. Если бы было иначе, люди давным-давно бы предсказали прогноз погоды на всю историю человечества вперед, а скачки курса рубля никогда не были бы неожиданностью.

Разумеется, предполагается, что верный ответ в этой задаче все-таки есть и он равен 10, и тогда «закон», связывающий все эти числа, – это

Однако возьмем любое другое значение – и мы также сможем найти закон, который бы обосновывал именно его:

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

Я считаю, верный ответ . Кто сможет оспорить?

Git-репозиторий

В прошлый раз меня ругали за то, что я выложил проект в виде архива в облаке, а не в виде кодов в репозиторий, поэтому в этот раз я исправляю эту свою ошибку: github.com/WieRuindl/Splines

Материал из MachineLearning.

Перейти к: навигация, поиск

Введение

Постановка математической задачи

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

( 1 )

Основная цель интерполяции — получить быстрый (экономичный) алгоритм вычисления значений для значений , не содержащихся в таблице данных.

Интерполируюшие функции , как правило строятся в виде линейных комбинаций некоторых элементарных функций:

где — фиксированный линейно независимые функции, — не определенные пока коэффициенты.

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

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

.

Тогда по заданным однозначно определяются коэффициенты .

Изложение метода

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

Кроме того, на границе при и ставятся условия

( 2 )

Будем искать кубический полином в виде

( 3 )

Из условия имеем

( 4 )

Вычислим производные:

и потребуем их непрерывности при :

( 5 )

Общее число неизвестных коэффициентов, очевидно, равно , число уравнений (4) и (5) равно . Недостающие два уравнения получаем из условия (2) при и :

Выражение из (5) , подставляя это выражение в (4) и исключая , получим

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

( 6 )

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

( 7 )

Условие эквивалентно условию и уравнению . Разностное уравнение (6) с условиями (7) можно решить методом прогонки, представив в виде системы линейных алгебраических уравнений вида , где вектор соответствует вектору , вектор поэлементно равен правой части уравнения (6), а матрица имеет следующий вид:

где и .

Метод прогонки

Метод прогонки, основан на предположении, что искомые неизвестные связаны рекуррентным соотношением:

( 8 )

Используя это соотношение, выразим и через и подставим в i-e уравнение:

,

где — правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать

Отсюда следует:

Из первого уравнения получим:

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

Пример: интерполирование неизвестной функции

Построим интерполянту для для функции , заданной следующим образом:

1	1.0002
2	1.0341
3	0.6
4	0.40105
5	0.1
6	0.23975

В результате интерполяции были рассчитаны следующие коэффициенты интерполянты:

				Отрезок
1,0002	-0,140113846	0,440979231	-0,266965385
1,0341	-0,291901538	-0,359916923	0,217718462
0,6	-0,22553	0,293238462	-0,266658462
0,40105	-0,100328462	-0,506736923	0,306015385
0,1	-0,134456154	0,411309231	-0,137103077

Ошибка интерполяции

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

от шага h, где .

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

причем константа в этом неравенстве является наилучшей из возможных

Пример: интерполяция синуса

Постром интерполянту функции на отрезке , взяв равномерно отстоящие узлы с шагом 0.5 и шагом 0.25, и сравним полученные результаты.

	Ошибка интерполяции	Оценка ошибки	Иллюстрация
	0.429685	3.(3)	Результат интерполяции sin(4x) с шагом 0.5
	0.005167	0.208(3)	Результат интерполяции sin(4x) с шагом 0.25

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

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

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

• А.А.Самарский, А.В.Гулин.  Численные методы М.: Наука, 1989.
• А.А.Самарский.  Введение в численные методы М.: Наука, 1982.
• Костомаров Д.П., Фаворский А.П.  Вводные лекции по численным методам
• Н.Н.Калиткин.  Численные методы М.: Наука, 1978.

См. также

• Интерполяция каноническим полиномом
• Тригонометрическая интерполяция
• Рациональная интерполяция
• Проблема выбора узлов для интерполяции
• Практикум ММП ВМК, 4й курс, осень 2008

Материал из MachineLearning.

Перейти к: навигация, поиск

Содержание

• 1 Введение
  • 1.1 Постановка математической задачи
• 2 Изложение метода
  • 2.1 Метод прогонки
  • 2.2 Пример: интерполирование неизвестной функции
• 3 Ошибка интерполяции
  • 3.1 Пример: интерполяция синуса
• 4 Список литературы
• 5 См. также

Введение

Постановка математической задачи

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

( 1 )

Основная цель интерполяции — получить быстрый (экономичный) алгоритм вычисления значений для значений , не содержащихся в таблице данных.

Интерполируюшие функции , как правило строятся в виде линейных комбинаций некоторых элементарных функций:

где — фиксированный линейно независимые функции, — не определенные пока коэффициенты.

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

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

.

Тогда по заданным однозначно определяются коэффициенты .

Изложение метода

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

Кроме того, на границе при и ставятся условия

( 2 )

Будем искать кубический полином в виде

( 3 )

Из условия имеем

( 4 )

Вычислим производные:

и потребуем их непрерывности при :

( 5 )

Общее число неизвестных коэффициентов, очевидно, равно , число уравнений (4) и (5) равно . Недостающие два уравнения получаем из условия (2) при и :

Выражение из (5) , подставляя это выражение в (4) и исключая , получим

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

( 6 )

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

( 7 )

Условие эквивалентно условию и уравнению . Разностное уравнение (6) с условиями (7) можно решить методом прогонки, представив в виде системы линейных алгебраических уравнений вида , где вектор соответствует вектору , вектор поэлементно равен правой части уравнения (6), а матрица имеет следующий вид:

где и .

Метод прогонки

Метод прогонки, основан на предположении, что искомые неизвестные связаны рекуррентным соотношением:

( 8 )

Используя это соотношение, выразим и через и подставим в i-e уравнение:

,

где — правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать

Отсюда следует:

Из первого уравнения получим:

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

Пример: интерполирование неизвестной функции

Построим интерполянту для для функции , заданной следующим образом:

Вводные значения для задачи интерполяции

1	1.0002
2	1.0341
3	0.6
4	0.40105
5	0.1
6	0.23975

В результате интерполяции были рассчитаны следующие коэффициенты интерполянты:

Результат интерполяции

				Отрезок
1,0002	-0,140113846	0,440979231	-0,266965385
1,0341	-0,291901538	-0,359916923	0,217718462
0,6	-0,22553	0,293238462	-0,266658462
0,40105	-0,100328462	-0,506736923	0,306015385
0,1	-0,134456154	0,411309231	-0,137103077

Ошибка интерполяции

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

от шага h, где .

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

причем константа в этом неравенстве является наилучшей из возможных

Пример: интерполяция синуса

Постром интерполянту функции на отрезке , взяв равномерно отстоящие узлы с шагом 0.5 и шагом 0.25, и сравним полученные результаты.

	Ошибка интерполяции	Оценка ошибки	Иллюстрация
	0.429685	3.(3)	Результат интерполяции sin(4x) с шагом 0.5
	0.005167	0.208(3)	Результат интерполяции sin(4x) с шагом 0.25

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

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

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

• А.А.Самарский, А.В.Гулин.  Численные методы М.: Наука, 1989.
• А.А.Самарский.  Введение в численные методы М.: Наука, 1982.
• Костомаров Д.П., Фаворский А.П.  Вводные лекции по численным методам
• Н.Н.Калиткин.  Численные методы М.: Наука, 1978.

См. также

• Интерполяция каноническим полиномом
• Тригонометрическая интерполяция
• Рациональная интерполяция
• Проблема выбора узлов для интерполяции
• Практикум ММП ВМК, 4й курс, осень 2008

bn = 2hn−3 (yn−1 − yn )− 2hn−−31 (yn−2 − yn−1 ). Такие системы быстро решаются специальными методами, с которыми мы познакомимся позднее.

Теорема 3.10. Пусть функция

y = f (x) имеет на отрезке [a, b] непрерывную про-

изводную четвертого порядка и M

4

= max

f (4)(x)

. Тогда для интерполяционного куби-

[a,b]

ческого сплайна S3 (x), удовлетворяющего граничным условиям типов 1, 2, 4 и 5, спра-

ведлива следующая оценка погрешности:

max

f

(x)− S

(x)

≤ C M

h4

, где C — некото-

рая константа.

[a, b]

3

4

max

не только сам аппроксимирует функцию y = f (x), но

Характерно, что сплайн S3 (x)

и его производные S3/ (x), S3// (x), S3/// (x)

приближают соответствующие производные функции

y = f (x).

Теорема 3.11. При выполнении условий теоремы 3.10 для указанных в ней

сплайнов справедливы неравенства max

f

(k )(x)− S (k )(x)

≤ C

k

M

4

h4−k , k = 1,2,3.

[a, b]

3

Пример. Пусть функция задана следующей таблицей:

	i				0				1				2				3				4
	xi			0				1				2				3				4
	yi				1.0				1.8				2.2				1.4				1.0

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

n = 4,

h = 1 = const . Запишем систему (3.10.4), реализующую эту схему в конкретном число-

вом

выражении.

Первое

уравнение

имеет

вид

h1−2 s0 + (h1−2 − h2−2 )s1

− h2−2 s2 =

= 2h−3

(y

− y

2

)− 2h−3 (y

0

− y

1

).

Вычисляя значения

h − h

−1

и подставляя

h = 1 ,

получим

2

1

1

i

i

s0 + 0

s1 − s2 = 2 (− 0.4)− 2 (− 0.8) s0 − s2

= 0.8. Аналогично,

последнее

уравнение дает

h−2 s

n−2

+ (h−2

− h−2 )s

n−1

− h−2 s

n

= 2h−3

(y

n−1

− y

n

)− 2h−3

(y

n−2

− y

n−1

), n = 4, n −1 = 3, n − 2 = 3.

n−1

n−1

n

n

n

n−1

Тогда

h3−2 s2

+ (h3−2

− h4−2 )s3

− h4−2 s4

= 2h4−3 (y3 − y4

)− 2h3−3 (y2

− y3 ) s2 − 0 s3 − s4 =

= 2 0.4 − 2 0.8 s2 − s4 = −0.8. Три внутренние узла дадут три следующих уравнения:

i = 1 : h1−1 s0 + 2 (h1−1 + h2−1 )s1 + h2−1 s 2 = 3[h1−2 (y1 − y 0 )+ h2−2 (y 2 − y1 )],
i = 2 : h2−1 s1 + 2	(h2−1	+ h3−1 )s2	+ h3−1 s3 = 3	[h2−2 (y 2	− y1 )+ h3−2 (y 3 − y 2	)],
i = 3 : h3−1 s 2 + 2	(h3−1	+ h4−1 )s3	+ h4−1 s 4 = 3	[h3−2 (y 3	− y 2 )+ h4−2 (y 4 − y 3	)].

Упрощая их аналогичным способом, что и выше, получим

s0 + 2 2 s1 + s2 = 3(0.8 + 0.4), s1 + 2 2 s2 + s3 = 3(0.4 + (− 0.8)), s2 + 2 2 s3 + s4 = 3(− 0.8 + (− 0.4)).

Тогда вся система имеет вид

88

s0

− s2

= 0.8,

s0

+ 4s1 + s2

= 3.6,

s1 + 4s2

+ s3

= −1.2,

s2

+ 4s3

+ s4

= −3.6,

s2

−

s4

= −0.8.

Решим эту систему методом Гаусса:

s0

− s2

= 0.8,

s0

− s2

= 0.8,

4s1 + 2s2

= 2.8,

s1 + 4s2

+ s3

= −1.2,

s1 + 4s2

+ s3

= −1.2,

−14s2

− 4s3

= 7.6,

s2 + 4s3 + s

= −3.6,

s2 + 4s3

+ s4 = −3.6,

4

s

2

−

s

4

= −0.8.

s

2

− s

4

= −0.8.

s0

− s2

= 0.8,

s0

− s2 = 0.8,

s1 + 4s2

+ s3

= −1.2,

s1 + 4s2 + s3

= −1.2,

s2

+ 4s3 + s4

= −3.6,

s2 + 4s3

+ s4

= −3.6,

− 4s3 − 2s4

= 2.8,

s3 + 0.5s4 = −0.7,

52s

3

+14s

4

= −42.8.

−12s

4

= −6.4.

Отсюда

s4 = 8 /15, s3

= −29 / 30,

s2

= −4 /15, s1

= 5 / 6,

s0

= 8 /15.

Запись формулы (3.10.1)

кубического интерполяционного многочлена Эрмита,

задавая yi−1 ,

yi , si−1 ,

si

однозначно

определяет сплайн на отрезке [xi−1 , xi ].

Эта формула уже учитывает, что интерполированная

функция

f (x) будет непрерывна и один раз дифференцируема везде на сетке

x0 , x1 ,…, xn .

Таким образом, на каждом из четырех отрезков

[0, 1], [1, 2], [2, 3], [3, 4] конкретный вид

сплайна

можно

установить

по

формуле

(3.10.1).

Например,

для

[0,1]

имеем

y0 = 1.0,

y1

= 1.8,

s0 = 8 / 15,

s1 = 5 / 6. Тогда

P

(x)

= y

0

(x1 − x)2 (2(x − x0 )+ h)

+ y /

(x1 − x)2 (x − x0 )

+

3,1

h3

0

h2

+ y

(x − x0 )2 (2(x1 − x)

+ h)

+ y /

(x − x0 )2 (x − x1 )

, h = x − x

0

.

1

h3

1

h2

1

P

(x)=1 (1 − x)2 (2(x − 0)+1)

+

8

(1 − x)2 (x − 0)

+1.8 (x − 0)2 (2(1 − x)+1)+

3,1

1

15

1

1

+

5

(x − 0)2 (x −1) = (x −1)2 (2x +1)+

8

(x −1)2 x +1.8x2 (3 − 2x)+

5

x2 (x −1)=

6

6

1

15

= −

7

x3 +

1

x2 +

8

x +1 = −0.233333x3 + 0.500000 x2

+ 0.533333x +1.000000.

30

15

2

Многочлены P3,i

на остальных отрезках рассчитываются аналогично. Именно:

x [1, 2],

x1

= 1,

x2

= 2,

y1

= 1.8,

y2 = 2.2,

s1

= 0.833333, s2 = −0.266667,

P

(x)= −0.233333x3 + 0.500000 x2

+ 0.533333x +1.000000.

3,2

x [2, 3],

x2

= 2, x3

= 3,

y2

= 2.2,

y3 =1.4, s2

= −0.266667, s3

= −0.966667,

P

(x)= 0.366667 x3 − 3.100000 x2

+ 7.700000 x − 3.800000.

3,3

x [3, 4],

x3

= 3,

x4

= 4,

y3

=1.4,

y4 =1.0,

s3

= −0.966667, s4

= 0.533333,

P

(x)= 0.366667 x3 − 3.100000 x2

+ 7.700000 x − 3.800000.

3,4

89

3.12. Лабораторная работа № 6. Аппроксимация функций кубическими сплайнами

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

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

Среди одиннадцати функций интерполяции и экстраполяции в среде Mathcad имеются средства линейной интерполяции (функция linterp(vx,vy,x) – значение в точке х линейного интерполяционного многочлена векторов vx и vy) и интерполяции сплайнами, осуществляемой функцией interp(K,X,Y,x) по значениям коэффициентов K векторов линейного, параболического или кубического сплайнов. Вектор коэффициентов соответствующего сплайна вычисляется подпрограммами lspline(vx,vy), pspline(vx,vy) и cspline(vx,vy).

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

Построим все три типа сплайнов и рассмотрим графическое представление данных с их помощью. Пусть

		0.593			0.53305
		0.598			0.53464
		0.605			0.54160
ORIGIN := 1		0.613			0.54324
X :=		Y :=	.
		0.619			0.54043
					0.55598
		0.627
		0.632			0.55843

Проведем сначала линейную интерполяцию и построим график результата. y(x):= linterp(X ,Y , x)

Для того чтобы исходные данные представлялись в виде квадратиков, надо войти в панель форматиро-

вания графиков (Formatting Currently Selected X-Y Plot), щелкнув дважды левой кнопкой мыши по любому месту графика. В открывшейся панели

для графика № 2 (trace 2) следует выбрать представление кривой в виде точек (points) квад-

ратной формы (box).

Сплайны соответствующей степени строятся аналогично. Наберем с клавиатуры

KLS := lspline(X ,Y ) KPS := pspline(X ,Y ) KCS := cspline(X ,Y )

y1(x):= interp(KLS, X ,Y , x) y2(x):= interp(KPS, X ,Y , x) y3(x):= interp(KCS, X ,Y , x).

На всем отрезке x [0.58,0.64] графики имеют вид

На концах отрезка вблизи точек x = 0.58 и x = 0.64 все четыре графика построим от-

дельно:
x := 0.56,0.5605…0.62	xx := 0.62,0.6205…0.64

Явно видимое различие трех графиков проявляется лишь при экстраполяции; внутри отрезка x [0.58,0.64] все три графика почти сливаются в один.

Из всех четырех графиков видно, что функция lspline(X,Y) строит вектор коэффициентов линейного сплайна дефекта нуль, то есть во всех внутренних узлах сетки состыкованы сама функция y = y(x) и ее первая производная. Линейный сплайн дефекта один – это ин-

терполяция ломаными, то есть функция linterp(X,Y,x), график которой приведен на самом первом рисунке.

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

Составим программу вычисления коэффициентов кубического сплайна с краевыми условиями четвертого типа – «отсутствие узла». Эти коэффициенты вычисляются решением системы линейных уравнений (3.10.4) методом прогонки. Система (3.10.4) отличается от классической трехдиагональной системы уравнений, используемой в методе прогонки, тем,

91

что в первом и последнем уравнении этой системы имеется по одному лишнему (ненулевому) коэффициенту. Поэтому формулы метода прогонки, приведенные в подразд. 5.4, необходимо несколько модернизировать. Модернизация касается способа вычисления x1 и x2 в

прямом ходе и x1 в обратном; все остальные переменные определяются стандартно по фор-

мулам (5.4.1) – (5.4.6).

Progonka1–подпрограмма, реализующая модернизированный метод прогонки. Параметры: вектор b содержит коэффициенты при неизвестных левой части системы уравнений, стоящие на главной диагонали, вектор a -под главной диагональю, вектор c -над главной

диагональю. Коэффициент p1 первого уравнения системы b1 x1 + c1 x2	+ p1 x3	= d1 находится в
cn , коэффициент p2 последнего уравнения p2 xn−2 + an xn−1 + bn xn = dn	в a1 .	Вектор d содер-

жит свободные члены (правые части системы уравнений). Подпрограмма выдает вектор x — вектор решений системы. Для правильной работы подпрограммы необходимо задание

ORIGIN :=1.

Следующая подпрограмма intspline вычисляет значение кубического сплайна с краевыми условиями четвертого типа для данного значения аргумента x . X и Y — векторы исходных данных, n — число точек сетки. Структура подпрограммы intspline совершенно прозрачна. Сначала вычисляются все коэффициенты диагональной системы (3.10.4), затем эта система решается модифицированным методом прогонки. Далее по значению аргумента x

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

при данном x по формуле (3.10.1) кубического интерполяционного многочлена Эрмита.

В заключение приведем графики полученных кубических сплайнов в одном масштабе:

KCS := cspline (X ,Y ) y3(x1):= interp(KCS , X ,Y , x1)

i := 1…100 x2i := X 1 + (X 7 − X 1 ) (i −1) 100

y1i := intspline (X, Y,7, x2i )

Видно, что оба графика на отрезке x [0.59,0.63] абсолютно совпадают.

В этом можно убедиться и сравнить полученные интерполированные значения f (x)

по двум разным программам:

x5 := 0.6 y5 := interp(KCS, X ,Y , x5) y5 = 0.536 y55 := intspline(X ,Y ,7, x5) y55 = 0.536

x6 := 0.63 y6 := interp(KCS, X ,Y , x6) y6 = 0.560 y56 := intspline(X ,Y ,7, x6) y56 = 0.560.

94

Задание № 1. Построить для таблично заданной функции y = f (x) кубический сплайн с граничными условиями любого из четырех типов на концах отрезка:

				Номера вариантов
	1		2		3		4		5
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
0.698	2.2234	0.100	1.1213	0.235	1.2080	0.095	1.0913	0.103	2.0128
0.706	2.2438	0.108	1.1316	0.240	1.2126	0.102	1.2349	0.108	2.0334
0.714	2.2645	0.119	1.1459	0.250	1.2217	0.104	1.2799	0.115	2.0607
0.727	2.2984	0.127	1.1565	0.255	1.2263	0.107	1.3514	0.120	2.0792
0.736	2.3222	0.135	1.1671	0.265	1.2355	0.110	1.4282	0.128	2.1072
0.747	2.3516	0.146	1.1819	0.280	1.2493	0.112	1.4826	0.136	2.1335
0.760	2.3869	0.157	1.1969	0.295	1.2633	0.116	1.6003	0.141	2.1492
0.769	2.4116	0.169	1.2134	0.300	1.2680	0.120	1.7321	0.150	2.1761
0.782	2.4478	0.175	1.2196	0.305	1.2726	0.125	1.8500	0.157	2.1805
				Номера вариантов
	6		7		8		9		10
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
0.296	3.2558	0.050	0.2079	0.902	1.2351	0.100	1.8378	0.400	1.6683
0.303	3.1764	0.052	0.2081	0.909	1.2369	0.104	1.8369	0.405	1.6664
0.310	3.1218	0.060	0.2090	0.919	1.2394	0.118	1.8335	0.410	1.6645
0.317	3.0482	0.065	0.2095	0.940	1.2448	0.139	1.8286	0.420	1.6607
0.323	2.9876	0.069	0.2099	0.944	1.2458	0.145	1.8272	0.429	1.6573
0.330	2.9195	0.075	0.2105	0.955	1.2486	0.158	1.8242	0.440	1.6532
0.339	2.8360	0.085	0.2116	0.965	1.2511	0.167	1.8221	0.449	1.6499
0.345	2.7771	0.090	0.2121	0.975	1.2537	0.185	1.8179	0.455	1.6476
0.352	2.6114	0.096	0.2127	1.010	1.2628	0.200	1.8145	0.465	1.6439
				Номера вариантов
	11		12		13		14		15
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
1.030	2.8011	0.200	0.0400	0.010	1.0101	0.100	0.8100	1.000	0.0000
1.080	2.9447	0.300	0.0899	0.105	1.1105	0.155	0.7140	1.510	0.2729
1.160	3.1899	0.350	0.1222	0.156	1.1681	0.220	0.6084	2.100	0.3533
1.230	3.4212	0.379	0.1431	0.200	1.2198	0.280	0.5184	2.750	0.3679
1.260	3.5254	0.415	0.1714	0.215	1.2378	0.375	0.3906	3.420	0.3595
1.330	3.7810	0.500	0.2474	0.289	1.3298	0.445	0.3080	3.915	0.3486
1.390	4.0149	0.596	0.3478	0.316	1.3645	0.510	0.2401	4.350	0.3380
1.450	4.1713	0.615	0.3693	0.390	1.4626	0.625	0.1406	4.800	0.3268
1.500	4.2390	0.700	0.4706	0.500	1.6152	0.770	0.0529	5.200	0.3171
				Номера вариантов
	16		17		18		19		20
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
0.100	0.0010	-2.100	-0.674	-1.000	1.368	1.000	0.000	1.000	0.6663
0.210	0.0441	-1.815	0.914	-0.010	1.990	5.000	0.999	1.300	0.5706
0.295	0.0868	-1.600	0.574	0.750	3.117	15.150	0.411	1.650	0.5429
0.348	0.1205	-1.420	0.827	1.900	7.686	27.900	-0.186	1.915	0.5887
0.419	0.1738	-1.290	0.931	2.815	17.692	35.100	-0.405	2.280	0.7256
0.475	0.2219	-0.750	0.983	3.100	23.198	44.000	-0.599	2.500	0.8262
					95

Соседние файлы в предмете Вычислительная математика

• #
• #

Макеты страниц

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

Рис. 11.8

Рис. 11.9

Наиболее широкое распространение на практике получили сплайны третьей степени (кубические сплайны) с дефектом, равным 1 или 2. Такие сплайны на каждом из частичных отрезков совпадают с кубическим многочленом:

и имеют на отрезке по крайней мере одну непрерывную производную

Термин «сплайн» происходит от английского слова (гибкая линейка, стержень) — названия приспособления, использовавшегося чертежниками для проведения гладких кривых через заданные точки. Бели гибкую стальную линейку поставить на ребро и, изогнув, зафиксировать ее положение в узловых точках (рис. 11.9), то получится механический аналог кубического сплайна. В самом деле, из курса сопротивления материалов известно, что уравнение свободного равновесия профиля линейки таково: Следовательно, в промежутке между двумя соседними узлами представляет собой многочлен третьей степени. В то же время отсутствие у линейки изломов свидетельствует о непрерывности касательной к графику функции и кривизны, т. е. производных

2. Интерполяционный сплайн.

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

Заметим, что на отрезке интерполяционный кубический сплайн однозначно определяется заданием значений самом деле, из равенства (11.31) вытекает следующая формула:

Здесь

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

3. Локальный сплайн.

Если в точках х, известны значения производной то естественно положить для всех Тогда на каждом частичном отрезке в соответствии с формулой (11.64) сплайн однозначно определяется значениями (поэтому его и называют локалъныл сплайном). Заметим, что он совпадает с кубическим интерполяционным многочленом Эрмита (11.31) для отрезка

Из неравенства (11.33) получается следующая оценка погрешности интерполяции локальным кубическим сплайном:

где Ащах максимальная из длин частичных отрезков.

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

Существуют и другие способы выбора коэффициентов а, приводящие к локальным сплайнам (кубический многочлен Бесселя, метод Акимы и др. [16]).

4. Глобальные способы построения кубических сплайнов.

Для того чтобы сплайн имел непрерывную на отрезке вторую производную необходимо выбирать наклоны а, так, чтобы в точках х, «стыка» многочленов совпадали значения их вторых производных:

Пользуясь формулой (11.64), найдем значение

Из подобной формулы, записанной

Таким образом, равенства (11.66) приводят к следующей системе уравнений относительно коэффициентов

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

1°. Если в граничных точках известны значения первой производной то естественно положить

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

2°. Если в граничных точках известны значения второй производной то можно наложить на сплайн граничные условия что приводит к следующим уравнениям:

(достаточно в равенстве (11.68) взять а в равенстве

3°. Полагая в уравнениях от того, выполнены ли эти условия для интерполируемой функции), придем к системе уравнений, определяющих так называемый естественный кубический сплайн.

4°. Часто нет, никакой дополнительной информации о значениях производных на концах отрезка. Один из применяемых в этой ситуациии подходов состоит в использовании условия «отсутствия узла». Выбор наклонов производят таким образом, чтобы для получаемого сплайна выполнялись условия Для этого достаточно потребовать совпадения в точках соответствующих третьих производных:

Эквивалентные алгебраические уравнения выглядят так:

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

5°. Если периодическая функция с периодом, равным о то систему (11-69) следует дополнить условиями

Существуют и другие подходы к заданию граничных условий (подробнее об этом см. [16]).

Пример 11.13. Для функции, заданной табл. 11.12, построй (естественный) кубический сплайн. В этом случае система уравнений для наклонов в точках записывается следующим образом:

Решал ее, получаем значения . Теперь на каждом частичном отрезке значения сплайна можно вычислить по формуле (11.64). Соответствующий график приведен на рис 11.10 (ср. с рис. 11.6 и 11.7).

Пример 11.14. Интерполируем функцию, заданную табл. 11.12, кубическим сплайном, используя условие «отсутствия узла». В этом случае уравнения останутся прежними, а уравнения (11 75) и (11.79) заменяются следующими:

Рис. 11.10

Решая систему получаем значения График соответствующего сплайна мало отличается от графика, изображенного на рис. 11.10.

5. Погрешность приближения кубическими сплайнами.

Теорема 11.9. Пусть функция имеет на отрезке непрерывную производную четвертою порядка и Тогда для интерполяционною кубическою сплайна удовлетворяющею граничным условиям типов 1°, 2°, 4° иди 5° (последнее — для случая периодической функции), справедлива следующая оценка погрешности:

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

Теорема 11.10. При выполнении условий теоремы 11.9 для указанных в ней сплайнов справедливы неравенства

Замечание. Благодаря большей простоте записи и благозвучному названию естественные сплайны получили значительное распространение. Однако искусственное наложение условий при интерполяции функций, которые этим условиям не удовлетворяют, приводит к значительной потере точности. Вместо четвертого порядка точности (как локальный кубический сплайн или кубические сплайны с граничными условиями типов 1°, 2°, 4°, 5°) естественный сплайн обладает лишь вторым порядком точности. Если использование естественного сплайна не вызвано какими-либо специальными причинами, то следует, по-видимому, отказаться от него в пользу кубического сплайна с граничным условием типа 4°.

1

Оглавление

• ПРЕДИСЛОВИЕ
• Глава 1. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И РЕШЕНИЕ ИНЖЕНЕРНЫХ ЗАДАЧ С ПРИМЕНЕНИЕМ ЭВМ
• § 1.2. Основные этапы решения инженерной задачи с применением ЭВМ
• § 1.3. Вычислительный эксперимент
• § 1.4. Дополнительные замечания
• Глава 2. ВВЕДЕНИЕ В ЭЛЕМЕНТАРНУЮ ТЕОРИЮ ПОГРЕШНОСТЕЙ
• § 2.1. Источники и классификация погрешностей результата численного решения задачи
• § 2.2. Приближенные числа. Абсолютная и относительная погрешности
• 2. Правила записи приближенных чисел.
• 3. Округление.
• § 2.3. Погрешности арифметических операций над приближенными числами
• § 2.4. Погрешность функции
• § 2.5. Особенности машинной арифметики
• 2. Представление целых чисел.
• 3. Представление вещественных чисел.
• 4. Арифметические операции над числами с плавающей точкой.
• 5. Удвоенная точность.
• 6. Вычисление машинного эпсилон.
• § 2.6. Дополнительные замечания
• Глава 3. ВЫЧИСЛИТЕЛЬНЫЕ ЗАДАЧИ, МЕТОДЫ И АЛГОРИТМЫ. ОСНОВНЫЕ ПОНЯТИЯ
• § 3.2. Обусловленность вычислительной задачи
• 2. Примеры плохо обусловленных задач.
• 3. Обусловленность задачи вычисления значения функции одной переменной.
• 4. Обусловленность задачи вычисления интеграла …
• 5. Обусловленность задачи вычисления суммы ряда.
• § 3.3. Вычислительные методы
• § 3.4. Корректность вычислительных алгоритмов
• § 3.5. Чувствительность вычислительных алгоритмов к ошибкам округления
• § 3.6. Различные подходы к анализу ошибок
• § 3.7. Требования, предъявляемые к вычислительным алгоритмам
• § 3.8. Дополнительные замечания
• Глава 4. МЕТОДЫ ОТЫСКАНИЯ РЕШЕНИЙ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
• § 4.2. Обусловленность задачи вычисления корня
• § 4.3. Метод бисекции
• § 4.4. Метод простой итерации
• § 4.5. Обусловленность метода простой итерации
• § 4.6. Метод Ньютона
• § 4.7. Модификации метода Ньютона
• § 4.8. Дополнительные замечания
• Глава 5. ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
• § 5.2. Нормы вектора и матрицы
• § 5.3. Типы используемых матриц
• § 5.4. Обусловленность задачи решения системы линейных алгебраических уравнений
• § 5.5 Метод Гаусса
• § 5.6. Метод Гаусса и решение систем уравнений с несколькими правыми частями, обращение матриц, вычисление определителей
• § 5.7. Метод Гаусса и разложение матрицы на множители. LU-разложение
• § 5.8. Метод Холецкого (метод квадратных корней)
• § 5.9. Метод прогонки
• § 5.10. QR-разложение матрицы. Методы вращений и отражений
• § 5.11. Итерационное уточнение
• § 5.12. Дополнительные замечания
• Глава 6. ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
• § 6.1. Метод простой итерации
• § 6.2. Метод Зейделя
• § 6.3. Метод релаксации
• § 6.4. Дополнительные замечания
• Глава 7. МЕТОДЫ ОТЫСКАНИЯ РЕШЕНИЙ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
• § 7.2. Метод простой итерации
• § 7.3. Метод Ньютона для решения систем нелинейных уравнений
• 7.4. Модификации метода Ньютона
• § 7.5. О некоторых подходах к решению задач локализации и отыскания решений систем нелинейных уравнений
• § 7.6. Дополнительные замечания
• Глава 8. МЕТОДЫ РЕШЕНИЯ ПРОБЛЕМЫ СОБСТВЕННЫХ ЗНАЧЕНИЙ
• § 8.2. Степенной метод
• § 8.3. Метод обратных итераций
• § 8.4. QR-алгоритм
• § 8.5. Дополнительные замечания
• Глава 9. МЕТОДЫ ОДНОМЕРНОЙ МИНИМИЗАЦИИ
• § 9.2. Обусловленность задачи минимизации
• § 9.3. Методы прямого поиска. Оптимальный пассивный поиск. Метод деления отрезка пополам. Методы Фибоначчи и золотого сечения
• § 9.4. Метод Ньютона и другие методы минимизация гладких функций
• § 9.5. Дополнительные замечания
• Глава 10. МЕТОДЫ МНОГОМЕРНОЙ МИНИМИЗАЦИИ
• § 10.1. Задача безусловной минимизации функции многих переменных
• § 10.2. Понятие о методах спуска. Покоординатный спуск
• § 10.3. Градиентный метод
• § 10.4. Метод Ньютона
• § 10.5. Метод сопряженных градиентов
• § 10.6. Метода минимизации без вычисления производных
• § 10.7. Дополнительные замечания
• Глава 11. ПРИБЛИЖЕНИЕ ФУНКЦИЙ И СМЕЖНЫЕ ВОПРОСЫ
• § 11.2. Интерполяция обобщенными многочленами
• § 11.3. Полиномиальная интерполяция. Многочлен Лагранжа
• § 11.4. Погрешность интерполяции
• § 11.5. Интерполяция с кратными узлами
• § 11.6. Минимизация оценки погрешности интерполяции. Многочлены Чебышева
• § 11.7. Конечные разности
• § 11.8. Разделенные разности
• § 11.9. Интерполяционный многочлен Ньютона. Схема Эйткена
• § 11.10. Обсуждение глобальной полиномиальной интерполяции. Понятие о кусочно-полиномиальной интерполяции
• § 11.11. Интерполяция сплайнами
• § 11.12. Понятие о дискретном преобразовании Фурье и тригонометрической интерполяции
• § 11.13. Метод наименьших квадратов
• § 11.14. Равномерное приближение функций
• § 11.15. Дробно-рациональные аппроксимации и вычисление элементарных функций
• § 11.16. Дополнительные замечания
• Глава 12. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ
• § 12.1. Простейшие формулы численного дифференцирования
• § 12.2. О выводе формул численного дифференцирования
• § 12.3. Обусловленность формул численного дифференцирования
• § 12.4. Дополнительные замечания
• Глава 13. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
• 13.2. Квадратурные формулы интерполяционного типа
• § 13.3. Квадратурные формулы Гаусса
• § 13.4. Апостериорные оценки погрешности. Понятие об адаптивных процедурах численного интегрирования
• § 13.5. Вычисление интегралов в нерегулярных случаях
• § 13.6. Дополнительные замечания
• Глава 14. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
• § 14.1. Задача Коши для дифференциального уравнения первого порядка
• § 14.2. Численные методы решения задачи Коши. Основные понятия и определения
• § 14.3. Использование формулы Тейлора
• § 14.4. Метод Эйлера
• § 14.5. Модификации метода Эйлера второго порядка точности
• § 14.6. Методы Рунге-Кутты
• § 14.7. Линейные многошаговые методы. Методы Адамса
• § 14.8. Устойчивость численных методов решения задачи Коши
• § 14.9. Неявный метод Эйлера
• § 14.10. Решение задачи Коши для систем обыкновенных дифференциальных уравнений и дифференциальных уравнений m-го порядка
• § 14.11. Жесткие задачи
• § 14.12. Дополнительные замечания
• Глава 15. РЕШЕНИЕ ДВУХТОЧЕЧНЫХ КРАЕВЫХ ЗАДАЧ
• § 15.1. Краевые задачи для одномерного стационарного уравнения теплопроводности
• § 15.2. Метод конечных разностей: основные понятия
• § 15.3. Метод конечных разностей: аппроксимации специального вида
• § 15.4. Понятие о проекционных и проекционно-разностных методах. Методы Ритца и Гадеркина. Метод конечных элементов
• § 15.5. Метод пристрелки
• § 15.6. Дополнительные замечания

This example shows how to use the csapi and csape commands from Curve Fitting Toolbox™ to construct cubic spline interpolants.

The CSAPI Command

The command

values = csapi(x,y,xx)

returns the values at xx of the cubic spline interpolant to the given data (x,y), using the not-a-knot end condition. This interpolant is a piecewise cubic function, with break sequence x, whose cubic pieces join together to form a function with two continuous derivatives. The «not-a-knot» end condition means that, at the first and last interior break, even the third derivative is continuous (up to round-off error).

Specifying only two data points results in a straight line interpolant.

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

Specifying three data points gives a parabola.

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

More generally, four or more data points give a cubic spline.

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

How to Check the Output of CSAPI

These look like nice interpolants, but how do we check that csapi performs as advertised?

We already saw that csapi interpolates, because we plotted the data points and the interpolant went right through those points. But to be sure that we get a cubic spline, it is best to start with data from a cubic spline of the expected sort and check whether csapi reproduces that cubic spline, i.e., gives back that cubic spline from which the data were taken.

Example: The Truncated Power Function

One simple example of a cubic spline function to check against is the truncated third power, i.e., the function

f(x)=((x-xi)+)3,

where xi is one of the breaks and the «+» subscript indicates the truncation function, provided by the command subplus:

 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

The truncated 3rd power is plotted below for the particular choice xi = 2. As expected, it is zero to the left of 2, and rises like (x-2)^3 to the right of 2.

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

Now we interpolate this particular cubic spline at the data sites 0:6, and plot the interpolant on top of the spline, in black.

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

The Interpolation Error

When comparing two functions, it is usually much more informative to plot their difference.

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

To put the size of their difference into context, you can also compute the maximum data value. This shows the error to be no worse than the inevitable round-off error.

A Truncated Power That Cannot be Reproduced

As a further test, we interpolate a truncated power whose csapi-produced interpolant at the sites 0:6 cannot coincide with it. For example, the first interior break of the interpolating spline is not really a knot since csapi uses the «not-a-knot» condition, hence the interpolant has three continuous derivatives at that site. This implies that we should not be able to reproduce the truncated 3rd power centered at that site since its third derivative is discontinuous across that site.

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

Since 1 is a first interior knot, it is not active for this interpolant.

The difference is as large as .18, but decays rapidly as we move away from 1. This illustrates that cubic spline interpolation is essentially local.

Using the ppform Instead of Values

It is possible to retain the interpolating cubic spline in a form suitable for subsequent evaluation, or for calculating its derivatives, or for other manipulations. This is done by calling csapi in the form

pp = csapi(x,y)

which returns the ppform of the interpolant. You can evaluate this form at some new points xx by the command

values = fnval(pp,xx)

You can differentiate the interpolant by the command

dpp = fnder(pp)

or integrate it by the command

ipp = fnint(pp)

which return the ppform of the derivative or the integral, respectively.

Example: Differentiating and Integrating the Interpolant

To show differentiation of an interpolant, we plot the derivative of this truncated power

f2′(x)=3((x-2)+)2,

(again in yellow) and then, on top of it, the derivative of our interpolant to the original truncated third power function (again in black).

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

Again, the more informative comparison is to plot their difference, and as before this is no bigger than the round-off error.

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

The second derivative of the truncated power is

f2′′(x)=6(x-2)+

A plot of the difference between this function and the second derivative of the interpolant to the original function shows that there are now jumps, but they are still within the round-off error.

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

The integral of the truncated power is

F2(x)=((x-2)+)4/4.

A plot of the difference between this function and the integral of the interpolant to the original function again shows that the errors are within the round-off error.

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

The CSAPE Command

Like csapi, the csape command provides a cubic spline interpolant to given data. However, it permits various additional end conditions. Its simplest version,

pp = csape(x,y)

uses the Lagrange end condition, which is a common alternative to the not-a-knot condition used by csapi. csape does not directly return values of the interpolant, but only its ppform.

For example, consider again interpolation to the function

f1(x)=((x-1)+)3,

which csapi fails to reproduce well. We plot the error of the not-a-knot interpolant returned by csapi (in black), along with the error of the interpolant obtained from csape (in red).

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

There is not much difference between the two interpolants in this case.

Other End Conditions: The ‘Natural’ Spline Interpolant

The csape command also provides ways to specify several other types of end conditions for an interpolating cubic spline. For example, the command

pp = csape(x,y,'variational')

uses the so-called ‘natural’ end conditions. This means that the second derivative is zero at the two extreme breaks.

This step shows how to apply ‘natural’ cubic spline interpolation to the function

f2(x)=((x-2)+)3,

and plot the error. The code below computes the ‘natural’ spline interpolant with an alternative argument syntax that is equivalent to the 'variational' argument: using 'second' specifies that csape should set the second derivative at the extreme data sites to the default value of 0.

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

Note the large error near the right end. This is due to the fact that the ‘natural’ end conditions implicitly insist on having a zero second derivative there.

Other End Conditions: Prescribing Second Derivatives

We can also explicitly use the correct second derivatives to get a small error. First, we compute the correct second derivative values of the truncated power at the endpoints.

endcond = 6*subplus(x([1 end])-2);

Then we create the interpolant, specifying that second derivatives at the endpoints are to be matched to the second derivative values we just computed. We do this by providing endcond(1) for the left endpoint condition, and endcond(2) for the right, along with the data values.

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

Other End Conditions: Prescribing Slopes

csape also permits specification of endpoint slopes. This is the clamped (or, complete) cubic spline interpolant. The statement

pp = csape(x,[sl,y,sr],'clamped')

creates the cubic spline interpolant to the data (x, y) that also has slope sl at the leftmost data site and slope sr at the rightmost data site.

Other End Conditions: Mixed End Conditions

It is even possible to mix these conditions. For example, our much-exercised truncated power function

f1(x)=((x-1)+)3

has slope 0 at x=0 and second derivative 30 at x=6 (the last data site).

Therefore, by matching the slope at the left end and the curvature at the right, we expect no error in the resulting interpolant.

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

Other End Conditions: Periodic Conditions

It is also possible to prescribe periodic end conditions. For example, the sine function is 2*pi-periodic and has the values [0 -1 0 1 0] at the sites (pi/2)*(-2:2). The difference, between the sine function and its periodic cubic spline interpolant at these sites, is only 2 percent. Not bad.

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

End Conditions Not Explicitly Covered by CSAPI or CSAPE

Any end condition not covered explicitly by csapi or csape can be handled by constructing the interpolant with the csape default side conditions, and then adding to it an appropriate scalar multiple of an interpolant to zero values and some side conditions. If there are two `nonstandard’ side conditions to be satisfied, you may have to solve a 2-by-2 linear system first.

For example, suppose that you want to compute the cubic spline interpolant s to the data

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

and enforce the condition

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

on the first and second derivatives of s at the point e.

The data were generated from a quartic polynomial that happens to satisfy this side condition with specific parameters

e = x(1);
a = 2; b = -3; c = 4;

To construct the interpolant that satisfies this specific condition, we first construct the interpolant with the default end conditions

and the first derivative of its first polynomial piece.

dp1 = fnder(fnbrk(pp1,1));

In addition, we construct the cubic spline interpolant to zero data values, specifying that it have a slope of 1 at e,

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

as well as constructing the first derivative of its first polynomial piece.

dp0 = fnder(fnbrk(pp0,1));

Then we compute lambda for both pp1 and pp0,

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

and construct the correct linear combination of pp1 and pp0 to get a cubic spline

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

that does satisfy the desired condition, as well as the default end condition at the right endpoint. We form this linear combination with the help of fncmb.

s = fncmb(pp0,(c-lam1)/lam0,pp1);

A plot of the interpolation error shows that s fits the quartic polynomial slightly better near e than the interpolant pp1 with the default conditions does.

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

If we want to enforce the condition

mu(s) := (D^3 s)(3) = 14.6

on the third derivative of the interpolant (the quartic satisfies this condition), then we construct an additional cubic spline interpolating to zero values, and with zero first derivative at the left endpoint, hence certain to be independent from pp0.

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

Then we find the coefficients d0 and d2 in the linear combination

s := pp1 + d0*pp0 + d2*pp2

that solves the linear system

lambda(s) = c

mu(s) = 14.6

Note that both pp0 and pp2 vanish at all interpolation sites, hence s will match the given data for any choice of d0 and d2.

For amusement, we use the MATLAB® encoding facility to write a loop to compute lambda(pp_j) and mu(pp_j), for j=0:2.

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

Given the values of lambda and mu for pp0, pp1, and pp2, we then solve for the coefficients that define the correct linear combination.

d = dd(:,[1,3])([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

For reassurance, we compare this error with the one obtained in complete cubic spline interpolation to this function:

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

The errors differ (and not by much) only near the end points, testifying to the fact that both pp0 and pp2 are sizable only near their respective end points.

As a final check, we verify that s satisfies the desired third derivative condition at 3.

This example shows how to use the csapi and csape commands from Curve Fitting Toolbox™ to construct cubic spline interpolants.

The CSAPI Command

The command

values = csapi(x,y,xx)

returns the values at xx of the cubic spline interpolant to the given data (x,y), using the not-a-knot end condition. This interpolant is a piecewise cubic function, with break sequence x, whose cubic pieces join together to form a function with two continuous derivatives. The «not-a-knot» end condition means that, at the first and last interior break, even the third derivative is continuous (up to round-off error).

Specifying only two data points results in a straight line interpolant.

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

Specifying three data points gives a parabola.

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

More generally, four or more data points give a cubic spline.

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

How to Check the Output of CSAPI

These look like nice interpolants, but how do we check that csapi performs as advertised?

We already saw that csapi interpolates, because we plotted the data points and the interpolant went right through those points. But to be sure that we get a cubic spline, it is best to start with data from a cubic spline of the expected sort and check whether csapi reproduces that cubic spline, i.e., gives back that cubic spline from which the data were taken.

Example: The Truncated Power Function

One simple example of a cubic spline function to check against is the truncated third power, i.e., the function

f(x)=((x-xi)+)3,

where xi is one of the breaks and the «+» subscript indicates the truncation function, provided by the command subplus:

 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

The truncated 3rd power is plotted below for the particular choice xi = 2. As expected, it is zero to the left of 2, and rises like (x-2)^3 to the right of 2.

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

Now we interpolate this particular cubic spline at the data sites 0:6, and plot the interpolant on top of the spline, in black.

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

The Interpolation Error

When comparing two functions, it is usually much more informative to plot their difference.

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

To put the size of their difference into context, you can also compute the maximum data value. This shows the error to be no worse than the inevitable round-off error.

A Truncated Power That Cannot be Reproduced

As a further test, we interpolate a truncated power whose csapi-produced interpolant at the sites 0:6 cannot coincide with it. For example, the first interior break of the interpolating spline is not really a knot since csapi uses the «not-a-knot» condition, hence the interpolant has three continuous derivatives at that site. This implies that we should not be able to reproduce the truncated 3rd power centered at that site since its third derivative is discontinuous across that site.

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

Since 1 is a first interior knot, it is not active for this interpolant.

The difference is as large as .18, but decays rapidly as we move away from 1. This illustrates that cubic spline interpolation is essentially local.

Using the ppform Instead of Values

It is possible to retain the interpolating cubic spline in a form suitable for subsequent evaluation, or for calculating its derivatives, or for other manipulations. This is done by calling csapi in the form

pp = csapi(x,y)

which returns the ppform of the interpolant. You can evaluate this form at some new points xx by the command

values = fnval(pp,xx)

You can differentiate the interpolant by the command

dpp = fnder(pp)

or integrate it by the command

ipp = fnint(pp)

which return the ppform of the derivative or the integral, respectively.

Example: Differentiating and Integrating the Interpolant

To show differentiation of an interpolant, we plot the derivative of this truncated power

f2′(x)=3((x-2)+)2,

(again in yellow) and then, on top of it, the derivative of our interpolant to the original truncated third power function (again in black).

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

Again, the more informative comparison is to plot their difference, and as before this is no bigger than the round-off error.

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

The second derivative of the truncated power is

f2′′(x)=6(x-2)+

A plot of the difference between this function and the second derivative of the interpolant to the original function shows that there are now jumps, but they are still within the round-off error.

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

The integral of the truncated power is

F2(x)=((x-2)+)4/4.

A plot of the difference between this function and the integral of the interpolant to the original function again shows that the errors are within the round-off error.

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

The CSAPE Command

Like csapi, the csape command provides a cubic spline interpolant to given data. However, it permits various additional end conditions. Its simplest version,

pp = csape(x,y)

uses the Lagrange end condition, which is a common alternative to the not-a-knot condition used by csapi. csape does not directly return values of the interpolant, but only its ppform.

For example, consider again interpolation to the function

f1(x)=((x-1)+)3,

which csapi fails to reproduce well. We plot the error of the not-a-knot interpolant returned by csapi (in black), along with the error of the interpolant obtained from csape (in red).

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

There is not much difference between the two interpolants in this case.

Other End Conditions: The ‘Natural’ Spline Interpolant

The csape command also provides ways to specify several other types of end conditions for an interpolating cubic spline. For example, the command

pp = csape(x,y,'variational')

uses the so-called ‘natural’ end conditions. This means that the second derivative is zero at the two extreme breaks.

This step shows how to apply ‘natural’ cubic spline interpolation to the function

f2(x)=((x-2)+)3,

and plot the error. The code below computes the ‘natural’ spline interpolant with an alternative argument syntax that is equivalent to the 'variational' argument: using 'second' specifies that csape should set the second derivative at the extreme data sites to the default value of 0.

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

Note the large error near the right end. This is due to the fact that the ‘natural’ end conditions implicitly insist on having a zero second derivative there.

Other End Conditions: Prescribing Second Derivatives

We can also explicitly use the correct second derivatives to get a small error. First, we compute the correct second derivative values of the truncated power at the endpoints.

endcond = 6*subplus(x([1 end])-2);

Then we create the interpolant, specifying that second derivatives at the endpoints are to be matched to the second derivative values we just computed. We do this by providing endcond(1) for the left endpoint condition, and endcond(2) for the right, along with the data values.

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

Other End Conditions: Prescribing Slopes

csape also permits specification of endpoint slopes. This is the clamped (or, complete) cubic spline interpolant. The statement

pp = csape(x,[sl,y,sr],'clamped')

creates the cubic spline interpolant to the data (x, y) that also has slope sl at the leftmost data site and slope sr at the rightmost data site.

Other End Conditions: Mixed End Conditions

It is even possible to mix these conditions. For example, our much-exercised truncated power function

f1(x)=((x-1)+)3

has slope 0 at x=0 and second derivative 30 at x=6 (the last data site).

Therefore, by matching the slope at the left end and the curvature at the right, we expect no error in the resulting interpolant.

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

Other End Conditions: Periodic Conditions

It is also possible to prescribe periodic end conditions. For example, the sine function is 2*pi-periodic and has the values [0 -1 0 1 0] at the sites (pi/2)*(-2:2). The difference, between the sine function and its periodic cubic spline interpolant at these sites, is only 2 percent. Not bad.

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

End Conditions Not Explicitly Covered by CSAPI or CSAPE

Any end condition not covered explicitly by csapi or csape can be handled by constructing the interpolant with the csape default side conditions, and then adding to it an appropriate scalar multiple of an interpolant to zero values and some side conditions. If there are two `nonstandard’ side conditions to be satisfied, you may have to solve a 2-by-2 linear system first.

For example, suppose that you want to compute the cubic spline interpolant s to the data

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

and enforce the condition

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

on the first and second derivatives of s at the point e.

The data were generated from a quartic polynomial that happens to satisfy this side condition with specific parameters

e = x(1);
a = 2; b = -3; c = 4;

To construct the interpolant that satisfies this specific condition, we first construct the interpolant with the default end conditions

and the first derivative of its first polynomial piece.

dp1 = fnder(fnbrk(pp1,1));

In addition, we construct the cubic spline interpolant to zero data values, specifying that it have a slope of 1 at e,

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

as well as constructing the first derivative of its first polynomial piece.

dp0 = fnder(fnbrk(pp0,1));

Then we compute lambda for both pp1 and pp0,

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

and construct the correct linear combination of pp1 and pp0 to get a cubic spline

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

that does satisfy the desired condition, as well as the default end condition at the right endpoint. We form this linear combination with the help of fncmb.

s = fncmb(pp0,(c-lam1)/lam0,pp1);

A plot of the interpolation error shows that s fits the quartic polynomial slightly better near e than the interpolant pp1 with the default conditions does.

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

If we want to enforce the condition

mu(s) := (D^3 s)(3) = 14.6

on the third derivative of the interpolant (the quartic satisfies this condition), then we construct an additional cubic spline interpolating to zero values, and with zero first derivative at the left endpoint, hence certain to be independent from pp0.

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

Then we find the coefficients d0 and d2 in the linear combination

s := pp1 + d0*pp0 + d2*pp2

that solves the linear system

lambda(s) = c

mu(s) = 14.6

Note that both pp0 and pp2 vanish at all interpolation sites, hence s will match the given data for any choice of d0 and d2.

For amusement, we use the MATLAB® encoding facility to write a loop to compute lambda(pp_j) and mu(pp_j), for j=0:2.

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

Given the values of lambda and mu for pp0, pp1, and pp2, we then solve for the coefficients that define the correct linear combination.

d = dd(:,[1,3])([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

For reassurance, we compare this error with the one obtained in complete cubic spline interpolation to this function:

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

The errors differ (and not by much) only near the end points, testifying to the fact that both pp0 and pp2 are sizable only near their respective end points.

As a final check, we verify that s satisfies the desired third derivative condition at 3.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325

#include <stdio.h> #include <cstring> #include <cmath> #include <iostream>   class knot { public:     double x, f, f2;     void Add(double arg, double func, double func2)     {         x = arg;         f = func;         f2 = func2;     }     knot() {} };   class vector { public:     double* x;     void Add(int m)     {         x = new double[m];     }     vector()     {     } };   knot* KnotArray; int n; // количество узлов интерполяции double** Coef; double* b;   //*************************************************************************/ //Решение системы уравнений с трехдиагональной матрицей //*************************************************************************/ void SolveTriDiag(double** TDM, double* F) {     double* alph = new double[n - 1];     double* beta = new double[n - 1];       int i;       alph[0] = -TDM[2][0] / TDM[1][0];     beta[0] = F[0] / TDM[1][0];       for (i = 1; i < n - 1; i++)     {         alph[i] = -TDM[2][i] / (TDM[1][i] + TDM[0][i] * alph[i - 1]);         beta[i] = (F[i] - TDM[0][i] * beta[i - 1]) / (TDM[1][i] + TDM[0][i] * alph[i - 1]);     }     b[n - 1] = (F[n - 1] - TDM[0][n - 1] * beta[n - 2]) / (TDM[1][n - 1] + TDM[0][n - 1] * alph[n - 2]);       for (i = n - 2; i > -1; i--)     {         b[i] = b[i + 1] * alph[i] + beta[i];     } }   //*************************************************************************/ //Построение таблицы коэффициентов кубического сплайна  y=f(x) //*************************************************************************/ int BuildSpline() {     double* a = new double[n - 1];     double* c = new double[n - 1];     double* d = new double[n - 1];     double* delta = new double[n - 1];     double* h = new double[n - 1];     double** TriDiagMatrix = new double* [3];       b = new double[n];       TriDiagMatrix[0] = new double[n];     TriDiagMatrix[1] = new double[n];     TriDiagMatrix[2] = new double[n];       double* f = new double[n];     double x3, xn;     int i;       if (n < 3)         return -1;       x3 = KnotArray[2].x - KnotArray[0].x;     xn = KnotArray[n - 1].x - KnotArray[n - 3].x;       for (i = 0; i < n - 1; i++)     {         a[i] = KnotArray[i].f;         h[i] = KnotArray[i + 1].x - KnotArray[i].x;         delta[i] = (KnotArray[i + 1].f - KnotArray[i].f) / h[i];         TriDiagMatrix[0][i] = i > 0 ? h[i] : x3;         f[i] = i > 0 ? 3 * (h[i] * delta[i - 1] + h[i - 1] * delta[i]) : 0;     }     TriDiagMatrix[1][0] = h[0];     TriDiagMatrix[2][0] = h[0];     for (int i = 1; i < n - 1; i++)     {         TriDiagMatrix[1][i] = 2 * (h[i] + h[i - 1]);         TriDiagMatrix[2][i] = h[i];     }     TriDiagMatrix[1][n - 1] = h[n - 2];     TriDiagMatrix[2][n - 1] = xn;     TriDiagMatrix[0][n - 1] = h[n - 2];         i = n - 1;     f[0] = ((h[0] + 2 * x3) * h[1] * delta[0] + powf(h[0], 2) * delta[1]) / x3;     f[n - 1] = (powf(h[i - 1], 2) * delta[i - 2] + (2 * xn + h[i - 1]) * h[i - 2] * delta[i - 1]) / xn;       SolveTriDiag(TriDiagMatrix, f);       /*Coef = new double* [4];     Coef[0] = new double[n - 1];     Coef[1] = new double[n - 1];     Coef[2] = new double[n - 1];     Coef[3] = new double[n - 1];     */     Coef = new double* [n - 1];     for (int count = 0; count < n - 1; count++)         Coef[count] = new double[4];     int j;       for (j = 0; j < n - 1; j++)     {         d[j] = (b[j + 1] + b[j] - 2 * delta[j]) / (h[j] * h[j]);         c[j] = 2 * (delta[j] - b[j]) / h[j] - (b[j + 1] - delta[j]) / h[j];           Coef[j][0] = a[j];         Coef[j][1] = b[j];         Coef[j][2] = c[j];         Coef[j][3] = d[j];     }     /*for (j = 0; j < n - 1; j++)//вывод значений коэффициентов полиномов     {         for (int i = 0; i < 4; i++)             printf("%lft", Coef[j][i]);         printf("n");     }*/     return 1; } /*************************************************************************/ //Построение таблицы коэффициентов кубического сплайна z=f(x) //*************************************************************************/ int BuildSpline_z() {     double* a = new double[n - 1];     double* c = new double[n - 1];     double* d = new double[n - 1];     double* delta = new double[n - 1];     double* h = new double[n - 1];     double** TriDiagMatrix = new double* [3];       b = new double[n];       TriDiagMatrix[0] = new double[n];     TriDiagMatrix[1] = new double[n];     TriDiagMatrix[2] = new double[n];       double* f = new double[n];     double x3, xn;     int i;       if (n < 3)         return -1;       x3 = KnotArray[2].x - KnotArray[0].x;     xn = KnotArray[n - 1].x - KnotArray[n - 3].x;       for (i = 0; i < n - 1; i++)     {         a[i] = KnotArray[i].f2;         h[i] = KnotArray[i + 1].x - KnotArray[i].x;         delta[i] = (KnotArray[i + 1].f2 - KnotArray[i].f2) / h[i];         TriDiagMatrix[0][i] = i > 0 ? h[i] : x3;         f[i] = i > 0 ? 3 * (h[i] * delta[i - 1] + h[i - 1] * delta[i]) : 0;     }     TriDiagMatrix[1][0] = h[0];     TriDiagMatrix[2][0] = h[0];     for (int i = 1; i < n - 1; i++)     {         TriDiagMatrix[1][i] = 2 * (h[i] + h[i - 1]);         TriDiagMatrix[2][i] = h[i];     }     TriDiagMatrix[1][n - 1] = h[n - 2];     TriDiagMatrix[2][n - 1] = xn;     TriDiagMatrix[0][n - 1] = h[n - 2];         i = n - 1;     f[0] = ((h[0] + 2 * x3) * h[1] * delta[0] + powf(h[0], 2) * delta[1]) / x3;     f[n - 1] = (powf(h[i - 1], 2) * delta[i - 2] + (2 * xn + h[i - 1]) * h[i - 2] * delta[i - 1]) / xn;       SolveTriDiag(TriDiagMatrix, f);       /*Coef = new double* [4];     Coef[0] = new double[n - 1];     Coef[1] = new double[n - 1];     Coef[2] = new double[n - 1];     Coef[3] = new double[n - 1];*/       Coef = new double* [n - 1];     for (int count = 0; count < n - 1; count++)         Coef[count] = new double[4];         int j;     for (j = 0; j < n - 1; j++)     {         d[j] = (b[j + 1] + b[j] - 2 * delta[j]) / (h[j] * h[j]);         c[j] = 2 * (delta[j] - b[j]) / h[j] - (b[j + 1] - delta[j]) / h[j];           Coef[j][0] = a[j];         Coef[j][1] = b[j];         Coef[j][2] = c[j];         Coef[j][3] = d[j];     }   } //*************************************************************************/ //Подсчет значения интерполянты в заданной точке //*************************************************************************/ double Interpolate(double x) {     //double result;     int i = 0;       while (KnotArray[i].x < x)         i++;     i--;     return Coef[i][0] + Coef[i][1] * (x - KnotArray[i].x) + Coef[i][2] * powf((x - KnotArray[i].x), 2) + Coef[i][3] * powf((x - KnotArray[i].x), 3);   }   //*************************************************************************/ //Загрузка данных //*************************************************************************/ int Load_Data() {     printf("Input filename with datan");     char FileName[20];     int i = 0;     FILE* File;     scanf_s("%s", &FileName, 20);     if (!fopen_s(&File, FileName, "r"))     {         printf("file is openedn");     }     else     {         printf("%s: file doesn't existn", FileName);         return -1;     }     double x, f, f2;     fscanf_s(File, "%d", &n);     KnotArray = new knot[n + 2];     while (!feof(File))     {         fscanf_s(File, "%lf%lf%lf", &x, &f, &f2);         KnotArray[i].Add(x, f, f2);         i++;         if (i == n)             return 1;     }     fclose(File);     return 1; }   int main() {     FILE* file_out;       int N;//количество значений функции-интерполянты     //double x = 0;     printf("Input the number of interpolant values:n");     scanf_s("%d", &N);     double* init = new double[N + 1];     double* fun = new double[N + 1];     double* fun_z = new double[N + 1];       if (Load_Data() != -1)     {         for (int i = 0; i <= N; i++)         {             init[i] = KnotArray[0].x + (KnotArray[n - 1].x - KnotArray[0].x) / N * i;           }         BuildSpline();         fun[0] = KnotArray[0].f;         fun_z[0] = KnotArray[0].f2;         fun[N] = KnotArray[n - 1].f;         fun_z[N] = KnotArray[n - 1].f2;         for (int i = 1; i < N; i++)         {             fun[i] = Interpolate(init[i]);             printf("i=%d f=%lfn", i, fun[i]);         }         BuildSpline_z();         for (int i = 1; i < N; i++)         {             fun_z[i] = Interpolate(init[i]);         }         if (fopen_s(&file_out, "D:interp.dat", "wb"))             printf("File could not be openedn");         else         {             for (int i = 0; i <= N; i++)             {                 printf("x:  %lfty_interp:   %lftz_interp:   %lfn", init[i], fun[i], fun_z[i]);                 fwrite(&init[i], sizeof(double), 1, file_out);                 fwrite(&fun[i], sizeof(double), 1, file_out);                 fwrite(&fun_z[i], sizeof(double), 1, file_out);             }         }     }     delete[] init;     delete[] fun;     delete[] fun_z;     for (int count = 0; count < n - 1; count++)         delete[] Coef[count];     system("pause");     return 0; }

Правила форума

В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте

его в существующую тему, а создайте новую в корневом разделе «Помогите решить/разобраться (М)».

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

Не ищите на этом форуме халяву

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

и указать конкретные затруднения.

Обязательно просмотрите тему

Правила данного раздела, иначе Ваша тема может быть удалена

или перемещена в Карантин, а Вы так и не узнаете, почему.

	Погрешность интерполяции сплайнами и оптимальная сетка 18.10.2007, 15:57
04/01/07 90	Всем привет! Вопрос к тем, кто изучал погрешность сплайн-интерполяции. Известно, что оценить ее сверху можно с использованием производных высших степеней интерполируемой функции. Используя эти сведения я попытался состряпать алгоритм формирования оптимальной сетки измерений. Но вскоре понял, что такой подход — тупик, поскольку производные (тем более высоких степеней) известны мне с погрешностью ~ +- 100% . Решил попробовать моделировать процедуру интерполяции для типовой (средне-обычной) функции и подбирать сетку таким образом, чтоб погрешность не превысила допустимой. Тут своии проблемы, связанные с взаимным влиянием интервалов. Т.е. при изменении шага измерения на одном из интервалов, меняется погрешность на другом и наоборот. В принципе, очень грубо, тут можно прикинуть в каких точках мерять, но слишком халтурно и «на глаз». Если кто-нибудь занимался подобной деятельностью, посоветуйте, пожалуйста, более-менее обоснованный способ расчета оптимальной сетки измерений. Заранее благодарен.

worm2	Re: Погрешность интерполяции сплайнами и оптимальная сетка 18.10.2007, 16:29
Заслуженный участник 01/08/06 2928 Уфа	oliva писал(а): Вопрос к тем, кто изучал погрешность сплайн-интерполяции. Известно, что оценить ее сверху можно с использованием производных высших степеней интерполируемой функции. Что значит «высших степеней»? 3-я это высшая или как? Для сплайн-интерполяции погрешность зависит от порядка сплайна, но не от количества точек. Например, для кубического сплайна на равномерной сетке: , где h — шаг сетки. Вас такая погрешность не устраивает? Может быть, вы имеете в виду полиномиальную интерполяцию (единым многочленом на всём отрезке) ?. Добавлено спустя 2 минуты 49 секунд: Или у Вас нет вообще никаких оценок для максимума модуля третьей производной?

oliva	Re: Погрешность интерполяции сплайнами и оптимальная сетка 18.10.2007, 16:46
04/01/07 90	worm2 писал(а): для кубического сплайна на равномерной сетке: , где h — шаг сетки. Вообще-то у меня другая формула для погрешности кубического сплайна. Но спорить не буду. Я так понял, что эта погрешность определяется свойствами (прерывностью, дифференциируемостью и т.п.) функции f. worm2 писал(а): Или у Вас нет вообще никаких оценок для максимума модуля третьей производной? Оценить максимум модуля я могу, но приблизительно. При этом получаю низкую точность определения шага измерения. Особенно это заметно в местах резкого изменения свойств f Добавлено спустя 3 минуты 28 секунд: worm2 писал(а): погрешность зависит от порядка сплайна, но не от количества точек. Например, для кубического сплайна на равномерной сетке: , где h — шаг Шаг и количество точек я рассматриваю как взаимоопределяющие факторы.

olga_helga	18.10.2007, 17:28
28/09/07 86	Стоп, я так понимаю что по сути это все дело сводится к нахождению новой таблицы значений, по которой будет построен более точный сплайн.Так? Добавлено спустя 4 минуты 17 секунд: А если так,то по сути это есть «задача зглаживания эксперементальных данных», основная идея которой построение новой таблицы данных, путем получения нового сглаженного значения сведением какой либо зависимости м/у и к линейной

oliva	18.10.2007, 17:29
04/01/07 90	olga_helga писал(а): Стоп, я так понимаю что по сути это все дело сводится к нахождению новой таблицы значений, по которой будет построен более точный сплайн.Так? Речь идет об интерполяции данных, полученных в процессе градуировки датчиков (в данном случае температуры). Характеристика этих датчиков в рабочем диапазоне меняется от почти линейной до сильно нелинейной. Нужно составить оптимальную сетку измерений, которая гарантирует нужную точность интерполяции при минимальном количестве измерений. Для этого я взял типовую характеристику и пытаюсь ее исследовать.

olga_helga	18.10.2007, 17:44
28/09/07 86	Вооо!!!! Тык вот, для зглаживания экспериментальных данных используют апроксимацию. Используя старую таблицу данных , получают новую , где -сглаженное значение. Так вот коим образом получают новую таблицу. 1)если зависимость линейная, то 2)если зависимость степенная,,то вводя замену переменных , получаем ,а -новый узел таблицы 3)если зависимость , то замена , -линейная зависимость, а -новая таблица. Ну вот что-то вроде етого.

oliva	18.10.2007, 17:52
04/01/07 90	olga_helga писал(а): для зглаживания экспериментальных данных используют апроксимацию. Да не нужно данные сглаживать Аппроксимацию со сглаживанием применяют когда данные получены с погрешностью и им не доверяют. Я решаю задачу интерполяции и измеренным точкам доверяю больше чем каким-либо расчетным. Интерполяция — нахождение кривой, которая ПРОХОДИТ через экспериментальные точки.

olga_helga	18.10.2007, 18:14
28/09/07 86

worm2	18.10.2007, 18:17
Заслуженный участник 01/08/06 2928 Уфа	Правильно ли я понял задачу: По известной функции f(x) на отрезке [a, b] найти минимальное число точек интерполяции на этом отрезке и построить сплайн на них, который гарантировал бы заданную погрешность определения f(x). Интересная задача. Можно попробовать решать тупым перебором на компьютере с каким-то шагом. Добавлено спустя 2 минуты 37 секунд: Я понимаю, что Ваша задача намного сложнее, но с этой упрощённой формулировки можно было бы начать.

olga_helga	18.10.2007, 18:21
28/09/07 86	Многочлен имеет n действительных корней,выражаемых формулой Добавлено спустя 3 минуты 14 секунд: В роде есть еще формулка для ,что то типа что узлы интерполяции выбираются как .Воть.Помогло?

Brukvalub	18.10.2007, 18:22
Заслуженный участник 01/03/06 13626 Москва	olga_helga писал(а): Помоему есть такая оценка погрешности интерполяции Это неверная оценка. Необходимо еще потребовать существование и умножить правую часть на

olga_helga	18.10.2007, 18:27
28/09/07 86	Извиняюсь!Я ошибочку сделала! В оригинале оценка для интерполирования n+1 раз непрерывно дифференцируемой ф-ции на отрезке [a,b] ф-ции f(x) многочленом Логранжа такая: .Вроде так.

Brukvalub	18.10.2007, 18:29
Заслуженный участник 01/03/06 13626 Москва

olga_helga	to worm2 18.10.2007, 18:31
28/09/07 86	Попробуйте последней из написанных формул воспользоваться.

oliva	Re: to worm2 18.10.2007, 19:36
04/01/07 90	olga_helga писал(а): Попробуйте последней из написанных формул воспользоваться. Я ж начал с того, что не могу с достаточной точностью знать производную выше первой степени, а в Вашей формуле присутствует производная степени (n+1). Что-же касается полиномов Чебышева — не все понял, но интересно. Сообщите пожалуйста источник, чтоб я вник. Добавлено спустя 5 минут 24 секунды: worm2 писал(а): Можно попробовать решать тупым перебором на компьютере с каким-то шагом. Над этим и «пыхтю» Проблема в том, что интервал измерения переменный и влияет на соседние. Т.е. при попытке подобрать один интервал (если подходить строго) нужно пересчитывать остальные, что усложняет, запутывает и ставит под сомнение результативность алгоритма.

Модераторы: Модераторы Математики, Супермодераторы