/* spline.c

*/

#include <float.h>

#include <stddef.h>

 

#include "spline.h"

 

#define TOL        (2.0*DBL_MIN/DBL_EPSILON)

#define square(x)  ((x)*(x))

 

 

double *spline (int n, double x[], double y[], double c[])

/*

 * Вычисление коэффициентов кубического интерполяционного сплайна,

 * записанного в виде:

 *

 * S(x) = w(i)*y(i+1) + w'(i)*y(i) +

 *        h(i)^2 * ((w(i)^3 - w(i))*c(i+1) + (w'(i)^3 - w'(i))*c(i)),

 *

 * где h(i)  = x(i+1) - x(i),

 *     w(i)  = (x - x(i)) / h(i),

 *     w'(i) = 1 - w(i).

 *

 * Аргументы:

 *

 * n - число заданных точек (узлов) (n >= 2),

 * x - массив абсцисс узлов в строго возрастающем порядке,

 * y - массив ординат узлов,

 * c - рабочий массив размером не менее 3*n элементов типа double.

 *     После возврата из функции spline в первых n элементах этого

 *     массива находятся коэффициенты сплайна c(i).

 *

 * Возвращаемое значение:

 *

 * Указатель на массив вычисленных коэффициентов c[] либо NULL, если

 * построить сплайн невозможно (n < 2, абсциссы узлов не упорядочены

 * либо среди них имеются совпадающие значения).

*/

{

   double *b = c + n;

   double *d = b + n;

 

   double t;

   int i;

 

   if (n < 2)

 

      return NULL;

 

   else if (n == 2)

 

      c[1] = c[0] = 0.0;

 

   else  {

      /* Формирование трехдиагональной системы линейных уравнений.

         b = диагональ, d = наддиагональ, c = правые части.

      */

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

         if ((d[i] = x[i+1] - x[i]) <= TOL) return (NULL);

         c[i+1] = (y[i+1] - y[i]) / d[i];

         if (i) {

            c[i] = c[i+1] - c[i];

            b[i] = 2.0 * (d[i-1] + d[i]);

         }

      }

      /* Граничные условия. Третьи производные в точках x[0] и x[n-1]

         вычисляются с помощью разделенных разностей.

      */

      b[0] = -d[0];

      b[n-1] = -d[n-2];

      if (n == 3)

         c[0] = c[n-1] = 0.0;

      else {

         c[0]   = c[2]/(x[3]-x[1]) - c[1]/(x[2]-x[0]);

         c[n-1] = c[n-2]/(x[n-1]-x[n-3]) - c[n-3]/(x[n-2]-x[n-4]);

         c[0]  *= square(d[0]) / (x[3]-x[0]);

         c[n-1] = -(c[n-1]*square(d[n-2]) / (x[n-1]-x[n-4]));

      }

 

      /* Прямой ход метода Гаусса */

 

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

         t = d[i] / b[i];

         b[i+1] -= t*d[i];

         c[i+1] -= t*c[i];

      }

 

      /* Обратная подстановка */

 

      c[n-1] /= b[n-1];

      for (i = n-2; i >= 0; i--)

         c[i] = (c[i] - d[i]*c[i+1]) / b[i];

   }

   return c;

}

 

/* splcvt.c

*/

#include <stddef.h>

 

#include "spline.h"

 

double *splcvt (int n, double x[], double y[], double c[])

/*

 * Преобразование сплайна, построенного с помощью функции spline,

 * к виду

 *

 * S(x) = y[i] + b[i]*(x-x[i]) + c[i]*(x-x[i])^2 + d[i]*(x-x[i])^3

 *

 * где x[i] <= x <= x[i+1].

 *

 * Аргументы:

 *

 * n - число заданных точек (узлов) (n >= 2),

 * x - массив абсцисс узлов в строго возрастающем порядке,

 * y - массив ординат узлов,

 * c - массив размером не менее 3*n элементов типа double.

 *     В первых n элементах этого массива должны находиться коэффициенты

 *     сплайна, рассчитанные функцией spline. После вызова splcvt массив

 *     содержит результат преобразования, т.е. коэффициенты b[i], c[i], d[i]

 *     в следующем порядке:

 *     - первые n элементов массива: коэффициенты c[i],

 *     - следующие n элементов:      коэффициенты b[i],

 *     - последние n элементов:      коэффициенты d[i].

 *

 * Возвращаемое значение:

 *

 * Указатель на массив вычисленных коэффициентов c[] либо NULL, если

 * построить сплайн невозможно (n < 2).

*/

{

   double *b = c + n;

   double *d = b + n;

   int i;

 

   if (n < 2)

 

      return NULL;

 

   else if (n == 2) {

 

      b[1] = b[0] = (y[1] - y[0]) / (x[1] - x[0]);

      d[1] = d[0] = 0.0;

 

   }

   else {

 

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

         d[i] = x[i+1] - x[i];

         b[i] = (y[i+1] - y[i]) / d[i];

      }

 

      b[n-1] = b[n-2] + d[n-2] * (c[n-2] + 2.0*c[n-1]);

 

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

         b[i] -= d[i] * (c[i+1] + 2.0*c[i]);

         d[i] = (c[i+1] - c[i]) / d[i];

         c[i] *= 3.0;

      }

      d[n-1] = d[n-2];

      c[n-1] *= 3.0;

 

   }

   return c;

}

 

/* seval1.c

*/

#include "spline.h"

 

double seval1 (int n, double z, double x[], double y[], double c[])

/*

 * Значение кубического сплайна в заданной точке.

 *

 * S(z) = y(i) + b(i)*(z-x(i)) + c(i)*(z-x(i))^2 + d(i)*(z-x(i))^3

 *

 * Номер интервала i выбирается таким образом, чтобы для значений z,

 * лежащих внутри диапазона табличных абсцисс, выполнялось условие:

 *

 *     x(i) <= z <= x(i+1).

 *

 * Если z < x(0), берется i = 0; если z >= x(n-1), берется i = n-1.

 *

 * Если по сравнению с предыдущим вызовом seval значение z не находится

 * в том же интервале [x(i), x(i+1)], тo для отыскания нужного интервала

 * применяется двоичный поиск.

 *

 * Аргументы:

 *

 * n - число заданных точек (узлов) (n >= 2),

 * z - абсцисса, для которой вычисляется значение сплайна,

 * x - массив абсцисс узлов в строго возрастающем порядке,

 * y - массив ординат узлов,

 * c - массив коэффициентов сплайна, рассчитанных функцией splcvt.

*/

{

   static int i = 0;

 

   double *b = c + n;

   double *d = b + n;

   double dx;

   int j, k;

 

   if (i >= n-1) i = 0;

 

   /* Двоичный поиск */

 

   if (z < x[i] || z > x[i+1]) {

      i = 0;

      j = n;

      while (j > i+1) {

         k = (i+j)/2;

         if (z < x[k])

            j = k;

         else

            i = k;

      }

   }

 

   /* Вычисление сплайна */

 

   dx = z - x[i];

   return (y[i] + dx*(b[i] + dx*(c[i] + dx*d[i])));

}

 

/* seval.c

*/

#include "spline.h"

 

#define square(x)  ((x)*(x))

 

double seval (int n, double z, double x[], double y[], double c[])

/*

 * Значение кубического сплайна в заданной точке.

 *

 * S(z) = w(i)*y(i+1) + w'(i)*y(i) +

 *        h(i)^2 * ((w(i)^3 - w(i))*c(i+1) + (w'(i)^3 - w'(i))*c(i)),

 *

 * где h(i)  = x(i+1) - x(i),

 *     w(i)  = (z - x(i)) / h(i),

 *     w'(i) = 1 - w(i).

 *

 * Номер интервала i выбирается таким образом, чтобы для значений z,

 * лежащих внутри диапазона табличных абсцисс, выполнялось условие:

 *

 *     x(i) <= z <= x(i+1).

 *

 * Если z < x(0), берется i = 0; если z > x(n-1), берется i = n-2.

 *

 * Если по сравнению с предыдущим вызовом seval значение z не находится

 * в том же интервале [x(i), x(i+1)], тo для отыскания нужного интервала

 * применяется двоичный поиск.

 *

 * Аргументы:

 *

 * n - число заданных точек (узлов) (n >= 2),

 * z - абсцисса, для которой вычисляется значение сплайна,

 * x - массив абсцисс узлов в строго возрастающем порядке,

 * y - массив ординат узлов,

 * c - массив коэффициентов сплайна, рассчитанных функцией spline.

*/

{

   static int i = 0;

 

   double h, w, w1;

   int j, k;

 

   if (i >= n-1) i = 0;

 

   /* Двоичный поиск */

 

   if (z < x[i] || z > x[i+1]) {

      i = 0;

      j = n-1;

      while (j > i+1) {

         k = (i+j)/2;

         if (z < x[k])

            j = k;

         else

            i = k;

      }

   }

 

   /* Вычисление сплайна */

 

   h = x[i+1] - x[i];

   w1 = 1.0 - (w = (z-x[i])/h);

   return (w*y[i+1] + w1*y[i] +

           square(h) * (w*(square(w)-1.0)*c[i+1] + w1*(square(w1)-1.0)*c[i]));

}

 

{====================}  UNIT SplinUNT; {======================}

 

{  В этот модуль включены подпрограммы SPLINE и SEVAL из книги

Форсайт Дж., Малькольм М., Моулер К., Машинные методы математи-

ческих вычислений: Пер. с анг.- М.: Мир, 1980. с. 91-94.

 

   Для повышения точности интерполяции используется арифметика

с плавающей точкой, определяемая типом

 

                    RealType = double

 

   Если Ваш компьютер не имеет арифметического сопроцессора, и

Вы хотите несколько увеличить  скорость вычислений ценой  ухуд-

шения точности, а также во всех других случаях, когда необходи-

мо сменить базовый вещественный тип (например, в целях экономии

памяти), Вы должны желаемым образом изменить определение вещес-

твенного типа, например: type RealType = real;

}

{-------------------------------------------------------------}

 

                         INTERFACE

{$N+,E+}

 

type

   RealType = double;              { базовый вещественный тип }

   RealTypeArray =

      array [1..(2*MaxInt) div Sizeof(RealType)] of RealType;

 

procedure Spline (N : word; var XI, YI, BI, CI, DI);

{

   Вычисляются коэффициенты B[i], C[i] и D[i], i = 1..N

   для кубического интерполяционного сплайна вида

 

   S(X) = Y[i]+B[i]*(X-X[i])+C[i]*(X-X[i])**2+D[i]*(X-X[i])**3

 

   для X[i] <= X <= X[i+1].

 

   ВХОДНАЯ ИНФОРМАЦИЯ:

 

   N - число заданных точек или узлов (N >= 2);

   X - абсциссы узлов в строго возрастающем порядке ;

   Y - ординаты узлов.

 

   ВЫХОДНАЯ ИНФОРМАЦИЯ:

 

   B, C, D, - массивы коэффициентов сплайна, записанного в указанной

              выше форме.

 

   Если обозначить дифференцирование символом ', то

   Y[i] = S(X[i]);

   B[i] = S'(X[i]);

   C[i] = S''(X[i])/2;

   D[i] = S'''(X[i])/6  (правосторонняя производная).

 

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

   сплайна.

}

{-------------------------------------------------------------}

 

function Seval (N : word; var U : RealType;

                var XI, YI, BI, CI, DI) : RealType;

{

   Эта функция вычисляет значение кубического сплайна:

 

   Seval = Y[i]+B[i]*(U-X[i])+C[i]*(U-X[i])**2+D[i]*(U-X[i])**3

 

   где X[i] <= U <= X[i+1]. Используется схема Горнера.

   Если U < X[1], то берется значение i = 1;

   если U > X[N], то i = N.

 

   ВХОДНАЯ ИНФОРМАЦИЯ:

 

   N - число заданных точек (N >= 2);

   U - абсцисса, для которой вычисляется значение сплайна;

   X, Y - массивы заданных абсцисс и ординат;

   B, C, D - массивы коэффициентов сплайна, вычисленные

             процедурой Spline.

 

   Если по сравнению с предыдущим вызовом Seval значение U не

   находится в том же интервале, то для отыскания нужного ин-

   тервала используется двоичный поиск.

}

{-------------------------------------------------------------}

 

                     IMPLEMENTATION

 

procedure Spline (N : word; var XI, YI, BI, CI, DI);

 

var

   X : RealTypeArray absolute XI;

   Y : RealTypeArray absolute YI;

   B : RealTypeArray absolute BI;

   C : RealTypeArray absolute CI;

   D : RealTypeArray absolute DI;

   T : RealType;

   i : integer;

 

begin { Spline }

 

   if N < 2 then Exit;

 

   if N = 2 then begin

 

      B[1] := (Y[2] - Y[1]) / (X[2] - X[1]);

      C[1] := 0;

      D[1] := 0;

      B[2] := B[1];

      C[2] := 0;

      D[2] := 0

 

   end

   else begin

 

    { Формирование трехдиагональной системы линейных уравнений:

      B = диагональ, D = наддиагональ, C = правые части. }

 

      D[1] := X[2] - X[1];

      C[2] := (Y[2] - Y[1]) / D[1];

      for i := 2 to N-1 do begin

         D[i] := X[i+1] - X[i];

         B[i] := 2 * (D[i-1] + D[i]);

         C[i+1] := (Y[i+1] - Y[i]) / D[i];

         C[i] := C[i+1] - C[i]

      end;

 

    { Граничные условия. Третьи производные в точках X[1] и X[N]

      оцениваются с помощью разделенных разностей. }

 

      B[1] := -D[1];

      B[N] := -D[N-1];

      if N = 3 then begin

         C[1] := 0;

         C[N] := 0

      end

      else begin

         C[1] := C[3]/(X[4]-X[2]) - C[2]/(X[3]-X[1]);

         C[N] := C[N-1]/(X[N]-X[N-2]) - C[N-2]/(X[N-1]-X[N-3]);

         C[1] :=  C[1] * sqr(D[1])/(X[4]-X[1]);

         C[N] := -C[N] * sqr(D[N-1])/(X[N]-X[N-3])

      end;

 

    { Прямой ход метода Гаусса }

 

      for i := 2 to N do begin

         T := D[i-1]/B[i-1];

         B[i] := B[i] - T*D[i-1];

         C[i] := C[i] - T*C[i-1]

      end;

 

    { Обратная подстановка }

 

      C[N] := C[N]/B[N];

      for i := N-1 downto 1 do C[i] := (C[i] - D[i]*C[i+1]) / B[i];

 

    { В C[i] теперь хранится значение Sigma(i), определяемое в книге:

      Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических

      вычислений. }

 

    { Вычисление коэффициентов сплайна B, C и D }

 

      B[N] := (Y[N]-Y[N-1])/D[N-1] + D[N-1]*(C[N-1]+2*C[N]);

      for i := 1 to N-1 do begin

         B[i] := (Y[i+1]-Y[i])/D[i] - D[i]*(C[i+1]+2*C[i]);

         D[i] := (C[i+1]-C[i])/D[i];

         C[i] := 3*C[i]

      end;

      C[N] := 3*C[N];

      D[N] := D[N-1]

 

   end { N <> 2 };

 

end { Spline };

 

{-------------------------------------------------------------}

 

function Seval (N : word; var U : RealType;

                var XI, YI, BI, CI, DI) : RealType;

 

var

   X : RealTypeArray absolute XI;

   Y : RealTypeArray absolute YI;

   B : RealTypeArray absolute BI;

   C : RealTypeArray absolute CI;

   D : RealTypeArray absolute DI;

   DX : RealType;

   j, k : integer;

 

const

   i : integer = 1;

 

begin { Seval }

 

   if i >= N then i := 1;

 

   if (U < X[i]) or (U > X[i+1]) then begin { Двоичный поиск }

      i := 1;

      j := N+1;

      repeat

         k := (i+j) div 2;

         if U < X[k] then j := k else i := k

      until j <= i+1

   end;

 

 { Вычисление сплайна }

 

   DX := U - X[i];

   Seval := Y[i] + DX*(B[i] + DX*(C[i] + DX*D[i]))

 

end { Seval };

 

{=================} END. { Unit SplinUNT } {==================}