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

Содержание 

1. Теоретические  предпосылки

2. Кубически  интерполяционный сплайн 

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

4. Построение  кубического сплайна функции

5. Тестирование  программы

6. Заключение

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

  1. Теоретические предпосылки

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

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

Теперь о математической части. Пусть заданы точки x1 , x2 , ..., xn  и соответствующие им значения y1 , y2 , ..., yn  функции f(x). На каждом из отрезков [xi , xi+1 ], i=1, 2, ..., n-1 функцию приближаем при помощи полинома третьей степени:

S(x) = yi  + c1i (x-xi ) + c2i (x-xi ) 2 + c3i (x-xi ) 3,

   xi  < x < xi+1

Для вычисления коэффициентов c1i , c2i , c3i , i = 1, 2, ..., n-1 решается система линейных уравнений, построенная  из условия непрерывности производной S'(x) в узлах сетки и дополнительных краевых условий на вторую производную, которые имеют вид:

2*S''1  + b1 *S''2  = b2

b3 *S''N-1  + 2*S''N  = b4

здесь возможны два случая. Случай первый, когда  известны значения первой производной  в краевых точках (y'1  = y'(x1 ), y'n  = y'(xn )), тогда следует положить:

b1  = 1, b2  = (6/(x2 -x1 )) * ((y2 -y1 ) / (x2 -x1 )-y'1 ),

b3  = 1, b4  = (6/(xn -xn-1 )) * (y'n  - (yn -yn-1 )/(xn -xn-1 ))

Случай второй, когда известны значения второй производной (y''1  = y''(x1 ), y''n  = y''(xn )), тогда полагаем:

b1  = 0, b2  = 2*y''1

b3  = 0, b4  = 2*y''N

Отметим, что  точки xi  должны быть отсортированы  по возрастанию.

На вход алгоритма (блок-схема см Приложение 1) подаются массивы x, y, на выходе получаем матрицу c:array[1..3,1..n]. Чтобы получить значение интерполируемой функции f в точке x необходимо в начале определить какому из отрезков разбиения она принадлежит (т.е. найти такое i, что xi  < x < xi+1 ) и затем воспользоваться формулой S(x) с сответствующими c1i , c2i , c3i , yi .  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

2. Кубически интерполяционный  сплайн

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

1) 

2) Для любого  частичного промежутка     -многочлен третьей степени

3)   

Для задания  надо определить 4 коэффициента для каждого промежутка  , т.е.  параметров.

Условия 1) требуют  чтобы во внутренних узлах сплайн и его производные до 2-го порядка были непрерывны.

Это дает  условия для определения параметров, еще  условие содержится в 3).

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

Возьмем в качестве граничных условия

4) 

Для построения кубического интерполяционного  сплайна могут быть использованы различные подходы. Проведем построение сплайна, исходя из условий 1) - 4). Из 1) и 2) следует, что ??? непрерывная функция, линейная на каждом  т.е. ??? - линейный сплайн.

Обозначив  , получаем

 (33)

для  .

Интегрируя (5), получаем

 (34)

 (35)

и  - постоянные интегрирования.

Условия 3) дают:

 (36)

Из (36) получаем:

Подставляя  и  в (7), получаем:

 (37)

После преобразования 

из (37) получаем

 (38)

Из (34) получаем

 (39)

Из (39) находим  односторонние пределы производной для узла  , 

 (40)

 (41)

Подставляя (40) и (41) в условие непрерывности  в узле  получаем :

 (42)

Дополняя (42) равенствами  из условия 4) :  , получаем систему уравнений относительно  вида :

 (43) 
с квадратной матрицей  .

 
и квадратной матрицей 

Координатами  вектора  являются значения  .

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

Матрица с диагональным преобладанием невырождена. Следовательно, система (42) однозначно разрешима, т.е. существует единственный кубический интерполяционный сплайн. Кроме условий 4) - условий "свободного провисания" интерполяционной кривой в точках  и  , могут быть известны наклоны интерполяционной кривой в граничных точках. Тогда условия на границах имеют вид:

 (44) 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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

   Пусть задана таблица значений функции f(x_i) = y_i ( ), в которой они расположены по возрастанию значений аргумента: x₀ < x₁ < … < x_n. Чтобы построить кубический сплайн, требуется определить коэффициенты a_i0, a_i1, a_i2, a_i3, которые задают интерполяционный кубический многочлен

   

   на  каждом интервале интерполирования [x_i-1, x_i],  .

   Таким образом, необходимо определить 4n коэффициентов a_ij ( ,  ), для чего требуется 4n уравнений. Необходимые уравнения определяются следующими условиями.

   1. Условия непрерывности функции:

   

   2. Условия непрерывности 1-х и  2-х производных функции:

   

   3. Граничные условия:

   

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

   Задача  определения кубического сплайна  существенно упрощается при использовании  многочлена Эрмита. Кубический многочлен  Эрмита на интервале [x_i-1, x_i] определяется с помощью значений функции y_i-1, y_i и ее производных y¢_i-1, y¢_i. Так как значения производных в общем случае могут быть неизвестны, обозначим их как y¢_i-1 = S_i-1; y¢_i = S_i. При построении сплайна переменные S_i называются наклонами сплайна в соответствующих точках x_i.

   Запишем многочлен Эрмита для интервала [x_i-1, x_i], где h_i =  x_i - x_i-1:

   

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

     

   Чтобы определить сплайн, нужно задать условия  непрерывности второй производной:

   

   Для записи этих условий в развернутом  виде определим кубический многочлен Эрмита на интервале [x_i, x_i+1], где h_i+1 =  x_i+1 - x_i:

   

   Определим вторые производные многочленов Q_i(x) и Q_i+1(x) в точке x = x_i:

                                 (4)

                             (5)

   Отсюда  условие непрерывности вторых производных имеет вид:

        (6)

   Это условие порождает систему линейных уравнений относительно наклонов сплайна S_i, которая содержит n - 1 уравнение и n + 1 переменную. Чтобы определить два недостающих уравнения используются граничные условия. Например, для естественного кубического сплайна:

     

   Указанные граничные условия могут быть получены из уравнения (5) для i = 0 и из уравнения (4) для i = n соответственно. В развернутом виде:

 

                               
 
 
 
 
 
 
 
 
 
 
 
 
 

4. Построение

кубического сплайна функции 

План: 

1) вывод расчётных  формул;

2) текст программы;

3) тестирование. 

Текст программы.

#include <iostream.h>

#include <fstream.h>

#include <conio.h>

#include <math.h>

#include <dos.h>

#include "mat_vec.h" // классы для работы с матрицами  и векторами

#include "progonka.h" // решение системы ур-ний (для  3-х диагональных матриц)

#include "funct.h" // второстепеннные функции программы  (рисование и т.д.) 

// "корень" программы

void spline (float step, int dop, int n, double* &x,double* &y,double* &x1,double* &y1) {

int k = 0;

matrica Sp(n, n-1);

for (int i = 1; i <= (n-1); i++) {

Sp(i,n) = 3*(y[i-1] - 2*y[i] + y[i+1])/pow(step,2);

Sp(i,i) = 4;

if (i < (n-1)) Sp(i,i+1) = 1;

if (i > 1) Sp(i,i-1) = 1;

}

float *tmp; 

progonka(Sp, tmp); // решение  системы уравнений методом прогонки

// (см. файл "progonka.h") 

vector a(n),b(n+1),c(n),d(n); // вычисление коэф-тов многочленов

b(1) = 0;

b(n+1) = 0;

for(int index = 0; index < n-1; index++)

b(index+2) = tmp[index];

delete [] tmp;

for (i = 1; i <= n; i++)

{

d(i) = y[i-1];

a(i) = (b(i+1)- b(i))/(-3*step);

c(i) = (y[i] - d(i) - pow(step,2)*b(i) + pow(step,3)*a(i) )/(-step);

}

i=0; 

//построение  графика сплайна при помощи  полученный коэф-тов (см. выше)

for (i=0; i < n; i++)