#include <stdio.h>
#include <math.h>
void spline( int n, double *x, double *y, double *b, double *c, double *d);
double seval( int n, double u, double *x, double *y, double *b, double *c,
double *d);
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;
register int ni;
double *px, *py, *pa, *pb, *pc, *pd;
/*-------------------------------------------------------------------------*/
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[1];
py = &y[1];
pd = &d[1];
pb = &b[1];
pc = &c[1];
ni = nm1 - 1;
while (ni--) {
*pd = *(px + 1) - *px;
*pb = 2. * (*(pd - 1) + *pd );
*(pc + 1) = ( *(py + 1) - *py ) / *pd;
*pc = *(pc + 1) - *pc;
px++; py++; pd++; pb++; pc++;
}
/*--------------------------------------------------------------------------
* Граничные условия. Третьи производные в точках 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] );
}
/*-------------*
* Прямой ход. *
*-------------*
*/
pb = &b[1];
pd = d;
pc = &c[1];
ni = nm1;
while (ni--) {
t = *pd / *(pb - 1);
*pb -= t * *pd;
*pc -= t * *(pc - 1);
pb++; pd++; pc++;
}
/*------------------------*
* Oбратная подстановка. *
*------------------------*
*/
c[nm1] = c[nm1] / b[nm1];
pc = &c[nm1-1];
pb = &b[nm1-1];
pd = &d[nm1-1];
ni = nm1;
while (ni--) {
*pc = ( *pc - *pd * *(pc + 1) ) / *pb;
pc--; pd--; pb--;
}
/*----------------------------------*
* В С(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;
while (ni--) {
*pb = ( *(py + 1) - *py ) / *pd - *pd * ( *(pc + 1) + 2. * *pc );
*pd = ( *(pc + 1) - *pc ) / *pd;
*pc = 3. * *pc;
pb++; pc++; pd++; py++;
}
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;
}
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;
static 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);
}
main() {
/*********************************
* Иллюстрирующая программа для *
* программ SPLINE & SEVAL. *
*********************************
*/
double x[10], y[10], b[10], c[10], d[10];
double s, u;
int i, n = 10;
for ( i = 0; i < n; i++ ) {
x[i] = i;
y[i] = x[i] * x[i] * x[i];
}
spline ( n, x, y, b, c, d );
u = 2.5;
s = seval( n, u, x, y, b, c, d );
printf(" %f %f \n ", s);
}