|
Текст подпрограммы и версий ids1r_c.zip |
Тексты тестовых примеров tids1r_c.zip |
Вычисление значений вторых производных вещественной табличной функции от одного переменного в узлах неравномерной сетки методом интерполирующих кубических сплайнов при известных условиях на вторые производные в концевых точках.
Пусть в узлах сетки xk : x1 <x2 < ... < xN заданы значения yk = y (xk), а также значения u и v, определяющие линейные условия y''1 = u*y''2, y''N = v*y''N - 1, где y''i = y''(xi). Интерполирующая функция y (x) на каждом из отрезков [xk, xk + 1] представляется кубическим полиномом Эрмита [1], построенным по значениям yi, y''i и yi + 1, y''i + 1, где y''i = y''(xi) определяются из условия непрерывности первой производной функции y (x) в узловых точках [2]. Это обеспечивает гладкость второго порядка y (x) на всем отрезке [x1, xN].
| 1. |
И.С.Березин, Н.П.Жидков, Методы вычислений, т. 1, Физматгиз, 1962. |
| 2. | Дж.Алберг, Э.Нильсон, Дж.Уолш, Теория сплайнов и ее приложения, M., "Мир", 1972. |
int ids1r_c (integer *n, real *x, real *y, real *y2, real *f,
real *g)
Параметры
| n - | число узлов сетки, n ≥ 3 (тип: целый); |
| x - | одномерный массив значений узлов сетки размерности n (тип: вещественный); |
| y - | одномерный массив значений табличной функции размерности n (тип: вещественный); |
| y2 - | одномерный массив вычисляемых значений вторых производных в узлах сетки размерности n; на входе y2 (1) = u, y2 (n) = v; (тип: вещественный); |
| f, g - | одномерные рабочие массивы размерности n (тип: вещественный). |
Версии: нет
Вызываемые подпрограммы: нет
Замечания по использованию: нет
На сетке x = (k-1)*π/5 , k = 1, ..., 6 брались значения yk = cos (xk), u = y1/y2, v = y6/y5
int main(void)
{
/* Builtin functions */
double cos(double);
/* Local variables */
extern int ids1r_c(int *, float *, float *, float *, float *, float *);
static float f[6], g[6];
static int k, n, i;
static float x[6], y[6], y2[6], pi;
pi = .62831853059999998f;
n = 6;
for (k = 1; k <= 6; ++k) {
x[k - 1] = (k - 1) * pi;
y[k - 1] = (float)cos(x[k - 1]);
/* l1: */
}
y2[0] = y[0] / y[1];
y2[n - 1] = y[n - 1] / y[n - 2];
ids1r_c(&n, x, y, y2, f, g);
for (i = 1; i <= 6; ++i) {
printf("\n %16.7e \n",y2[i-1]);
}
return 0;
} /* main */
Результат: y2 = ( -1.03, -0.84, -0.32, 0.32, 0.84, 1.03 )