public static void
spline( int n, double[] x, double[] y, double[] b, double[] c, double[] d)
{
/*--------------------------------------------------------------------------
* ********************** Подпрограмма SPLINE ***************************
*__________________________________________________________________________
* Вычисляются коэффициенты B(i), C(i) и D(i), i = 1,2,...,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))
* D(i) = S'''(X(i)) ( правосторонняя производная )
*__________________________________________________________________________
* C помощью сопровождающей подпрограммы-функции SEVAL
* можно вычислить значения сплайна.
*--------------------------------------------------------------------------
*/
double t;
int nm1, ib, i;
int ni;
double[] px = new double[n+1];
double[] py = new double[n+1];
double[] pb = new double[n+1];
double[] pc = new double[n+1];
double[] pd = new double[n+1];
/*-------------------------------------------------------------------------*/
nm1 = n -1;
if ( n < 2 ) return;
if ( n >= 3 ) {
/*------------------------------*
* B = диагональ *
* D = наддиагональ *
* С = правые части *
*------------------------------*
*/
d[0] = x[1] - x[0];
c[1] = ( y[1] - y[0] ) / d[0];
px = x;
py = y;
pd = d;
pb = b;
pc = c;
for( i = 0; i < n-2; i++) {
pd[1+i] = px[2+i] - px[1+i];
pb[1+i] = 2. * (pd[i] + pd[1+i]);
pc[2+i] = (py[2+i] -py[1+i])/pd[1+i];
pc[1+i] = pc[2+i] - pc[1+i];
}
int j=1;
/*--------------------------------------------------------------------------
* Граничные условия. Третьи производные в точках X(i) и X(N)
* вычисляются с помощью разделенных разностей.
*--------------------------------------------------------------------------
*/
b[0] = -d[0];
b[nm1] = -d[nm1 - 1];
c[0] = 0.;
c[nm1] = 0.;
if ( n != 3 ) {
c[0] = c[2] / ( x[3] - x[1] ) - c[1] / ( x[2] - x[0] );
c[nm1] = c[n - 2] / ( x[nm1] - x[n - 3] ) - c[n - 3] / ( x[n - 2]
- x[n - 4]);
c[0] = c[0] * d[0] * d[0] / ( x[3] - x[0] );
c[nm1] = -c[nm1] * d[nm1-1] * d[nm1-1] / ( x[n - 1] - x[n - 4] );
}
/*-------------*
* Прямой ход. *
*-------------*
*/
pd = d; // внимание, здесь другой сдвиг
ni = n-1;
for( i = 0; i < n-1; i++) {
t = pd[i] / pb[i];
pb[i+1] -= t * pd[i];
pc[i+1] -= t * pc[i];
}
/*------------------------*
* Oбратная подстановка. *
*------------------------*
*/
c[nm1] = c[nm1] / b[nm1];
ni = nm1;
for( i = n-1; i > 0; i--) {
pc[i-1] = ( pc[i-1] - pd[i-1] * pc[i] ) / pb[i-1];
}
/*----------------------------------*
* В С(i) теперь храниться величина *
* SIGMA(i), определенная в $ 4.4 . *
*__________________________________*
* Вычислить коэффициенты полиномов.*
*----------------------------------*
*/
b[nm1] = ( y[nm1] - y[nm1 - 1] ) / d[nm1 - 1] + d[nm1 - 1] *
( c[nm1 - 1] + 2. * c[n - 1] );
pb = b; pc = c; pd = d; py = y;
ni = nm1;
for( i = 0; i < n-1; i++) {
pb[i] = (py[i+1] - py[i])/pd[i] - pd[i] * (pc[i+1] + 2.0* pc[i]);
pd[i] = (pc[i+1] - pc[i])/pd[i];
pc[i] = 3.0*pc[i];
}
c[nm1] = 3. * c[nm1];
d[nm1] = d[nm1 - 1];
return;
}
b[0] = ( y[1] - y[0] ) / ( x[1] - x[0] );
c[0] = 0.;
d[0] = 0.;
b[1] = b[0];
c[1] = 0.;
d[1] = 0.;
return;
}
static double seval( int n, double u, double[] x, double[] y, double[] b, double[] c, double[] d)
{
/***********************----------------------------------------------------
* Подпрограмма SEVAL * Эта подпрограмма вычисляет значение
*********************** кубического сплайна.
* 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(i), то берется значение i = 1.
* Если U >= X(i), то берется значение i = N.
*__________________________________________________________________________
* ВХОДНАЯ ИНФОРМАЦИЯ:
* N - число заданных точек.
* U - абцисса, для которой вычисляется значение сплайна.
* Х, Y - массивы заданных абцисс и ординат.
* B, C, D - массивы коэффициентов сплайна, вычисленные
* подпрограммой SPLINE.
*__________________________________________________________________________
* Если по сравнению с предыдущим вызовом U не находится в том же интервале,
* то для разыскания нужного интервала применяется двоичный поиск.
*--------------------------------------------------------------------------
*/
int j, k;
double dx, seval;
int i=1;
if ( i >= n ) i = 1;
if( (u<x[i-1])||(u>x[i]) )
{
/******************
* Двоичный поиск *
*/
i = 0;
j = n;
while ( j > i+1 )
{
k = ( i + j ) / 2;
if ( u < x[k-1] ) j = k;
if ( u >= x[k-1] ) i = k;
}
}
/********************
* Вычислить сплайн *
********************
*/
dx = u - x[i-1];
seval = y[i-1] + dx * ( b[i-1] + dx * ( c[i-1] + dx * d[i-1] ) );
return( seval);
}