Pers.narod.ru. Алгоритмы. Построение интерполяционного кубического сплайна

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

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

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

(1)

Функцию φ(x) и будем использовать для аппроксимации зависимости f(x), заданной в узлах х0, x1,...,xn значениями f0,f1,...,fn.

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

(2)

где ai, bi, ci и di - коэффициенты сплайна, определяемые из дополнительных условий; i = 1, 2,..., n - номер сплайна.

В отличие от полиномиальной интерполяции, когда вся аппроксимируемая зависимость описывается одним полиномом, при сплайновой интерполяции на каждом интервале [хi-1, хi] строится отдельный полином третьей степени вида (2) со своими коэффициентами.

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

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

(3)

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

(4, 5)

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

(6, 7)

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

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

(8, 9),

где

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

(10,11)

Из условий непрерывности производных (4, 5) при переходе в точке xi от i-го к (i+1)-му сплайну с учетом выражений (10, 11) получим следующие соотношения:

(12, 13)

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

(14, 15)

Соотношения (8-9), (12-15) представляют собой полную систему линейных алгебраических уравнений относительно коэффициентов сплайнов ai, bi, ci и di. Но, прежде чем решать эту систему, выгодно преобразовать её так, чтобы неизвестными была только одна группа коэффициентов сi.

Из уравнения (13) коэффициенты di выразим через коэффициенты сi:

(16)

Объединяя уравнения (8-9) с соотношением (16), представим коэффициенты bi также через коэффициенты ci:

(17)

После подстановки выражений (16-17) в соотношение (12) получим уравнение, в которое входят только неизвестные коэффициенты ci. Для симметричности записи в полученном уравнении уменьшим значение индекса i на единицу

(18), где 2<=i<=n.

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

. (19)

Таким образом, n-1 уравнение вида (18) вместе с условиями (14) и (19) образует систему линейных алгебраических уравнений для определения коэффициентов ci. Коэффициенты di и bi вычисляются после нахождения ci по формулам (16) и (17), коэффициенты ai равны значениям аппроксимируемой функции в узлах в соответствии с формулой (8).

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

Рассмотрим алгоритм метода прогонки применительно к решаемой задаче. Для сокращения записи уравнение (18) представим в виде

(20), где

(21)

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

Из уравнения (20) при i = 2 с учетом граничного условия (14) установим связь коэффициента c2 с коэффициентом c3

(22), где к2 и l2- прогоночные коэффициенты,

.

Затем, подставляя выражение (22) в уравнение (20) при i=3, получим линейную связь коэффициента с3 с коэффициентом с4:

.

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

(23)

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

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

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

(24)

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

а также полагая c2 <> 0, задаем начальные коэффициенты k1=0 и l1=0. Затем по формуле (24) вычислим все n пар прогоночных коэффициентов ki и li.

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

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

(25)

Далее последовательно применим формулу (23) при i = n-1, n-2, ... ,2 и вычислим значения искомых величин cn-1, cn-2, ... ,c2. Эта процедура называется обратным ходом метода прогонки.

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

После определения всех коэффициентов ci другие коэффициенты сплайнов вычисляются по формулам (8), (16) и (17), после чего аппроксимирующую функцию φ(x) можно рассчитать с помощью соотношения (2) в любой точке х на интервале [х0, хn].

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

Программа запрашивает количество точек n, а затем известные значения (xi, yi), i =1,2,...,n.

{Кубический интерполяционный сплайн}
program Spline;
uses crt; {модуль управления экраном}
type vector=array [0..100] of real; {Нумеруем точки с нуля}
var x,y,c:vector;
    x0,x9,h,x1,p,p1,p2,e:real;
    n,i:integer;
 
procedure InputData (var n:integer; var x,y:vector); {Ввод исходных данных}
var i:integer;
begin
 write ('Введите количество точек n (3<n<100): ');
 read (n);
 for i:=1 to n do begin
  write ('Введите точку x[',i,'], f(x[',i,']): ');
  read (x[i-1], y[i-1]);
 end;
end;
 
procedure Coeff(n:integer; var x,f,c:vector);
{Вычисление коээфициентов сплайна}
var i,j,m:integer;
    a,b,r:real;
    k:vector;
begin
 {Прямой ход прогонки}
 k[1]:=0; c[1]:=0;
 for i:=2 to n do begin
  j:=i-1;
  m:=j-1;
  a:=x[i]-x[j];
  b:=x[j]-x[m];
  r:=2*(a+b)-b*c[j];
  c[i]:=a/r;
  k[i]:=(3.0*((f[i]-f[j])/a-(f[j]-f[m])/b)-b*k[j])/r;
 end;
 {Обратный ход прогонки}
 c[n]:=k[n];
 for i:=n-1 downto 2 do c[i]:=k[i]-c[i]*c[i+1];
end;
 
procedure Spl (n:integer; var x,f,c:vector; x1:real; var p,p1,p2:real);
{Построение сплайна. x,f - исходные данные, c - вектор коэффициентов,
наденный процедурой Coeff, x1 - значение x, для которого строим сплайн,
p - значение сплайна в точке, p1,p2 - 1-я и 2-я производные}
var i,j:integer;
    a,b,d,q,r:real;
begin
 i:=1;
 while (x1>x[i]) and (i<>n) do i:=i+1; {Ищем номер соседнего узла}
 {Промежуточные переменные и коэффициенты}
 j:=i-1; a:=f[j]; b:=x[j]; q:=x[i]-b;
 r:=x1-b; p:=c[i]; d:=c[i+1];
 b:=(f[i]-a)/q - (d+2*p)*q/3.0;
 d:=(d-p)/q*r;
 {Считаем значения сплайна и его производных:}
 p1:=b+r*(2*p+d);
 p2:=2*(p+d);
 p:=a+r*(b+r*(p+d/3.0))
end;
 
begin
 clrscr; {очистить экран}
 writeln ('Построение кубического интерполяционного сплайна');
 InputData (n,x,y);
 
 Coeff (n,x,y,c); {Нашли коэффициенты C с помощью прогонки}
 {Строим таблицу значений сплайна}
 writeln ('Значение X':19,'Значение F(X)':19,'Сплайн':19,'Ошибка':19);
 for i:=0 to n do begin
  Spl (n,x,y,c,x[i],p,p1,p2);
  e:=abs(y[i]-p);
  writeln (x[i]:19:8,y[i]:19:8,p:19:8,e:19:8);
 end;
 reset (input); readln;
end.

Вот пример вывода этой программы.

Построение кубического интерполяционного сплайна
Введите количество точек n (3<n<100): 4
Введите точку x[1], f(x[1]): 0 0
Введите точку x[2], f(x[2]): 1 1
Введите точку x[3], f(x[3]): 2 2
Введите точку x[4], f(x[4]): 3 3
         Значение X      Значение F(X)             Сплайн             Ошибка
         0.00000000         0.00000000         0.00000000         0.00000000
         1.00000000         1.00000000         1.00000000         0.00000000
         2.00000000         2.00000000         2.00000000         0.00000000
         3.00000000         3.00000000         3.00000000         0.00000000
         0.00000000         0.00000000         0.00000000         0.00000000

Приведённая ниже версия программы отличается только тем, что сплайн строится по точкам, полученным вычислением функции, заданной подпрограммой f(x). Вводится границы и шаг по оси X для построения сплайна. Погрешность считается относительной.

{Кубический интерполяционный сплайн}
program Spline;
uses crt; {модуль управления экраном}
type vector=array [0..100] of real; {Нумеруем точки с нуля}
var x,y,c:vector;
    x0,x9,h,x1,p,p1,p2,e:real;
    n,i:integer;

function f(x:real):real; {функция варианта}
begin
 f:=sqr(sqr(x))-5*sqr(x)-x+1;
end;

procedure Input (var n:integer; var x,y:vector); {Ввод исходных данных}
var i:integer;
begin
 writeln ('Введите границы по оси X для построения сплайна:');
 read (x0,x9);
 writeln ('Введите шаг по X для построения значений сплайна:');
 read (h);
 n:=round((x9-x0)/h+1);
 x[0]:=x0; y[0]:=f(x0);
 for i:=1 to n do begin {Вычислить значения функции}
  x[i]:=x[i-1]+h;
  y[i]:=f(x[i]);
 end;
end;

procedure Coeff(n:integer; var x,f,c:vector);
{Вычисление коээфициентов сплайна}
var i,j,m:integer;
    a,b,r:real;
    k:vector;
begin
 {Прямой ход прогонки}
 k[1]:=0; c[1]:=0;
 for i:=2 to n do begin
  j:=i-1;
  m:=j-1;
  a:=x[i]-x[j];
  b:=x[j]-x[m];
  r:=2*(a+b)-b*c[j];
  c[i]:=a/r;
  k[i]:=(3.0*((f[i]-f[j])/a-(f[j]-f[m])/b)-b*k[j])/r;
 end;
 {Обратный ход прогонки}
 c[n]:=k[n];
 for i:=n-1 downto 2 do c[i]:=k[i]-c[i]*c[i+1];
end;

procedure Spl (n:integer; var x,f,c:vector; x1:real; var p,p1,p2:real);
{Построение сплайна. x,f - исходные данные, c - вектор коэффициентов,
наденный процедурой Coeff, x1 - значение x, для которого строим сплайн,
p - значение сплайна в точке, p1,p2 - 1-я и 2-я производные}
var i,j:integer;
    a,b,d,q,r:real;
begin
 i:=1;
 while (x1>x[i]) and (i<>n) do i:=i+1; {Ищем номер соседнего узла}
 {Промежуточные переменные и коэффициенты}
 j:=i-1; a:=f[j]; b:=x[j]; q:=x[i]-b;
 r:=x1-b; p:=c[i]; d:=c[i+1];
 b:=(f[i]-a)/q - (d+2*p)*q/3.0;
 d:=(d-p)/q*r;
 {Считаем значения сплайна и его производных:}
 p1:=b+r*(2*p+d);
 p2:=2*(p+d);
 p:=a+r*(b+r*(p+d/3.0))
end;

begin
 clrscr; {очистить экран}
 writeln ('Построение кубического интерполяционного сплайна');
 writeln ('Функция Y(x)=x^4+5x^2-x+1');
 Input (n,x,y);

 Coeff (n,x,y,c); {Нашли коэффициенты C с помощью прогонки}
 {Строим таблицу значений сплайна}
 writeln ('Значение X':19,'Функция':19,'Сплайн':19,'Отн.ошибка':19);
 for i:=0 to n do begin
  Spl (n,x,y,c,x[i],p,p1,p2);
  e:=abs(y[i]-p)/abs(y[i]);
  writeln (x[i]:19:8,y[i]:19:8,p:19:8,e:19:8);
 end;
end.

Рейтинг@Mail.ru

вверх гостевая; E-mail
Hosted by uCoz