|
Текст подпрограммы и версий is01r_c.zip , is01d_c.zip |
Тексты тестовых примеров tis01r_c.zip , tis01d_c.zip |
Построение одномерного сглаживающего кубического сплайна.
Пусть на сетке x1 < x2 < ... < xN заданы измеренные значения f1, f2, ..., fN некоторой гладкой функции. Требуется найти дважды непрерывно дифференцируемую функцию s (x), которая минимизирует функционал
xN
∫ ( g''(x) )2 dx
xi
среди всех функций g (x) таких, что
N
∑ [ ( g(xi) - fi ) / ( dfi ) ]2 ≤ SM ,
i=1
где dfi > 0 и SM ≥ 0 - заданные числа.
Числа dfi характеризуют точность задания дискретных данных fi, а SM - требуемую меpу их сглаживания.
Решением задачи является сглаживающий кубический сплайн, который на каждом интервале [xi, xi + 1] представляется в виде
s(x) = c3 i h3 + c2 i h2 + c1 i h + yi , h = x - xi ,
где c3 i , c2 i , c1 i , i = 1, 2, ..., N - 1 суть коэффициенты, а yi - значения сглаживающего сплайна в узлах xi , i = 1, 2, ..., N.
Reinsch C.H. Smoothing by Spline Functions. Numerishe Mathematik, 1967, v.10, 177 - 183.
Морозов B.A. O задаче дифференцирования и некоторых алгоритмах приближения экспериментальной информации, в сб. "Вычислительные методы и программирование", 1970, вып.ХIV, Изд - во МГУ, 46 - 62.
int is01r_c (real *x, real *f, real *df, integer *nx, real *sm,
real *y, real *c, real *rab, integer *ierr)
Параметры
| x - | вещественный массив длины nx, в котоpом задаются значения узлов сетки такие, что x (i) < x (i + 1), i = 1, ..., nx - 1; |
| f - | вещественный массив длины nx, в котоpом задаются значения сглаживаемой дискретной функции; |
| df - | вещественный массив длины nx, в котоpом задаются значения весов df (i) > 0 , i = 1, ..., nx; |
| nx - | заданное число узлов, nx ≥ 2 (тип: целый); |
| sm - | неотрицательное число, задающее меpу сглаживания (тип: вещественный); |
| y - | вещественный массив длины nx значений сглаживающего сплайна в узлах сетки; |
| c - | двумерный вещественный массив размера 3* (nx - 1) значений коэффициентов сглаживающего сплайна; |
| rab - | двумерный рабочий массив размера 7* (nx + 2) (тип: вещественный); |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы: |
| ierr=65 - | если nx < 2; |
| ierr=66 - | если массив x не упорядочен строго по возрастанию. |
Версии
| is01d_с - | построение одномерного сглаживающего кубического сплайна с повышенной точностью. |
Вызываемые подпрограммы
| utis10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы is01r_c. |
| utis11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы is01d_c. |
Замечания по использованию: нет
int main(void)
{
/* Initialized data */
static float x[5] = { 1.f,2.f,3.f,4.f,5.f };
static float f[5] = { 1.f,4.f,20.f,16.f,25.f };
static float df[5] = { .2f,.2f,10.f,.2f,.2f };
/* Builtin functions */
double sqrt(double);
/* Local variables */
extern int is01r_c(float *, float *, float *, int *, float *, float *,
float *, float *, int *);
static int ierr;
static float c__[12] /* was [3][4] */;
static int i__;
static float y[5], sm;
static int nx;
static float rab[49] /* was [7][7] */;
#define c___ref(a_1,a_2) c__[(a_2)*3 + a_1 - 4]
nx = 5;
sm = nx - 1.f - (float)sqrt((float)((float) (nx + nx)));
is01r_c(x, f, df, &nx, &sm, y, c__, rab, &ierr);
printf("\n %5i \n",ierr);
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %16.7e \n",y[i__-1]);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %16.7e %16.7e %16.7e %16.7e \n",
c___ref(i__, 1), c___ref(i__, 2), c___ref(i__, 3), c___ref(i__, 4));
}
return 0;
} /* main */
Результаты:
ierr = 0
y(1) = .999 ; y(2) = 4.003 ; y(3) = 10.85 ; y(4) = 16.003 ; y(5) = 24.999
c__(1, 1) = 1.785 c__(1, 2) = 5.442 c__(1, 3) = 6.000 c__(1, 4) = 6.558
c__(2, 1) = 0.000 c__(2, 2) = 3.657 c__(2, 3) = -3.099 c__(2, 4) = 3.657
c__(3, 1) = 1.219 c__(3, 2) = -2.252 c__(3, 3) = 2.252 c__(3, 4) = -1.219