АЛГОРИТМЫ

Новости

Рассылка новостей

Форум

AlgoPascal

Редактор блок-схем

Статьи

О сайте

Контакты
Содержание - Интерполяция - Интерполяция кубическими сплайнами

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

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

Теперь о математической части. Пусть заданы точки 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  должны быть отсортированы по возрастанию.

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

Для удобства работы эти операции выделены в две отдельные функции. Функция Spline3BuildTable по граничным условиям и точкам строит таблицу коэффициентов, в которую также помещает информацию об опорных точках (в элементы с первыми индексами 0 и 4), а функция Spline3Interpolate по готовой таблице коэффициентов позволяет интерполировать функцию.

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

Если нашли ошибку в алгоритме - сообщите!


Реализации алгоритма на различных языках:

Реализация алгоритма на C++
Реализация алгоритма на Delphi
Реализация алгоритма на Visual Basic 6

Блоксхемы:

Посмотреть блок-схему алгоритма
Скачать блок-схему алгоритма

Реализация алгоритма на AlgoPascal:

unit Spline3Unit;

interface Spline3BuildTable, Spline3Interpolate;

implementation
(*************************************************************************
Этот модуль содержит две функции - Spline3BuildTable и Spline3Interpolate.

Процедура Spline3BuildTable служит  для  постройки  таблицы  коэффициентов
кубического сплайна по заданным точкам и граничным условиям, накладываемым
на производные. Функция принимает параметры:
    N       - число отрезков разбиения, не меньше одного.
    DiffN   - тип граничного условия. "1" соответствует граничным условиям
              накладываемым на первые производные, "2" - на вторые.
    xs      - массив абсцисс опорных точек с номерами от 0 до n.
    ys      - массив ординат опорных точек с номерами от 0 до n.
    BoundL  - левое граничное условие. Если DiffN равно 1, то первая произ-
              водная  на  левой  границе  равна BoundL, иначе BoundL равна
              вторая.
    BoundR  - аналогично BoundL
Функция возвращает:
    ctbl    - в этот массив помещается таблица коэффициентов сплайна.


Функция Spline3Interpolate по построенной таблице коэффициентов  вычисляет
значение интерполируемой функции в заданной точке. Параметры:
    N       - число отрезков разбиения
    c       - таблица коэффициентов, построенная функцией
    X       - точка, в которой ведем расчет

*************************************************************************)

procedure Spline3BuildTable(
            N       :   Integer;
    const   DiffN   :   Integer;
            x       :   array of Real;
            y       :   array of Real;
    const   BoundL  :   Real;
    const   BoundR  :   Real;
    out     ctbl    :   array of array of Real);
var
    C       :   Boolean;
    E       :   Integer;
    G       :   Integer;
    Tmp     :   Real;
    nxm1    :   Integer;
    I       :   Integer;
    J       :   Integer;
    DX      :   Real;
    DXJ     :   Real;
    DYJ     :   Real;
    DXJP1   :   Real;
    DYJP1   :   Real;
    DXP     :   Real;
    DYP     :   Real;
    YPPA    :   Real;
    YPPB    :   Real;
    PJ      :   Real;
    b1      :   Real;
    b2      :   Real;
    b3      :   Real;
    b4      :   Real;
begin
    //сортируем массив точек по возрастанию
    g:=((n+1) div 2);
    repeat
        i:=g;
        repeat
            j:=i-g;
            c:=True;
            repeat
                if x[j]<=x[j+g] then
                begin
                    c:=False;
                end
                else
                begin
                    Tmp:=x[j];
                    x[j]:=x[j+g];
                    x[j+g]:=Tmp;
                    Tmp:=y[j];
                    y[j]:=y[j+g];
                    y[j+g]:=Tmp;
                end;
                j:=j-1
            until not((j>=0)and(C));
            i:=i+1
        until not(i<=n);
        g:=g div 2;
    until not(g>0);

    //подготовка таблиц/коэффициентов
    SetBounds(ctbl, [0, 4], [0, N]);
    N:=N+1;
    
    if DiffN=1 then
    begin
        b1:=1;
        b2:=(6/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-BoundL);
        b3:=1;
        b4:=(6/(x[n-1]-x[n-2]))*(BoundR-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
    end
    else
    begin
        b1:=0;
        b2:=2*BoundL;
        b3:=0;
        b4:=2*BoundR;
    end;

    nxm1:=n-1;
    if n>=2 then
    begin
        if n>2 then
        begin
            dxj:=x[1]-x[0];
            dyj:=y[1]-y[0];
            j:=2;
            while j<=nxm1 do
            begin
                dxjp1:=x[j]-x[j-1];
                dyjp1:=y[j]-y[j-1];
                dxp:=dxj+dxjp1;
                ctbl[1,j-1]:=dxjp1/dxp;
                ctbl[2,j-1]:=1-ctbl[1,j-1];
                ctbl[3,j-1]:=6*(dyjp1/dxjp1-dyj/dxj)/dxp;
                dxj:=dxjp1;
                dyj:=dyjp1;
                j:=j+1
            end;
        end;
        ctbl[1,0]:=-b1/2;
        ctbl[2,0]:=b2/2;
        if n<>2 then
        begin
            j:=2;
            while j<=nxm1 do
            begin
                pj:=ctbl[2,j-1]*ctbl[1,j-2]+2;
                ctbl[1,j-1]:=-ctbl[1,j-1]/pj;
                ctbl[2,j-1]:=(ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,J-2])/pj;
                j:=j+1
            end;
        end;
        yppb:=(b4-b3*ctbl[2,nxm1-1])/(b3*ctbl[1,nxm1-1]+2);
        i:=1;
        while i<=nxm1 do
        begin
            j:=n-i;
            yppa:=ctbl[1,j-1]*yppb+ctbl[2,j-1];
            dx:=x[j]-x[j-1];
            ctbl[3,j-1]:=(yppb-yppa)/dx/6;
            ctbl[2,j-1]:=yppa/2;
            ctbl[1,j-1]:=(y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]*dx)*dx;
            yppb:=yppa;
            i:=i+1
        end;
        
        //дополняем таблицу коэффициентов
        for i:=1 to n do
        begin
            ctbl[0,i-1]:=y[i-1];
            ctbl[4,i-1]:=x[i-1];
        end;
    end;
end;

function Spline3Interpolate(
    const   N       :   Integer;
    const   c       :   array of array of Real;
    const   X       :   Real):Real;
var
    I   :   Integer;
    L:      Integer;
    Half:   Integer;
    First:  Integer;
    Middle: Integer;
begin
    //двоичный поиск отрезка для интерполяции
    L:=N;
    First:=0;
    while L>0 do
    begin
        Half:=L div 2;
        Middle:=First+Half;
        if C[4,Middle]<X then
        begin
            First:=Middle+1;
            L:=L-Half-1;
        end
        else
            L:=Half;
    end;
    I:=First-1;
    if I<0 then
        I:=0;

    Result:=c[0,I] + (X-c[4,i])*
        (C[1,I] + (X-c[4,i])*(C[2,I] + C[3,i]*(X-c[4,i])));
end;

end.

 


Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU