Текст подпрограммы и версий ip04r_c.zip , ip04d_c.zip |
Тексты тестовых примеров tip04r_c.zip , tip04d_c.zip |
Бикубическая интерполяция функции двух переменных, заданной вместе с частными производными первого порядка и смешанной производной второго порядка в вершинах элементарного сеточного прямоугольника.
Пусть X1L, X1U, X2L и X2U следующим образом определяют вершины элементарного сеточного прямоугольника:
(X1L,X2U) | | (X1U,X2U) ______________________________________ | 4 3 | | (X1,X2) | | 1 2 | ______________________________________ (X1L,X2L) | | (X1U,X2L)
Эти вершины пронумерованы на рисунке цифрами, начиная с нижней левой вершины против часовой стрелки. Пусть в каждой из вершин известны значения интерполируемой функции
Y = f(x1, x2) , обеих ее частных производных первого порядка: Y1 = ∂f / ∂x1 , Y2 = ∂f / ∂x2 и смешанной производной Y12 = ∂2 f / ∂x1 ∂x2 .
Подпрограмма ip04r_c при помощи бикубической интерполяции вычисляет в точке (X1,X2), заданной внутри исходного элементарного сеточного прямоугольника, значение функции
AY = f (X1, X2) , обеих ее частных производных первого порядка: AY1 = ∂ f (X1, X2) / ∂x1 , AY2 = ∂ f (X1, X2) / ∂x2 и смешанной производной AY12 = ∂2 f (X1, X2) / ∂x1 ∂x2.
Используются следующие формулы:
4 4 AY = ∑ ∑ ci j t i-1 u j-1 ; i=1 j=1 4 4 AY1 = ∑ ∑ (i - 1) ci j t i-2 u j-1 ; i=1 j=1 4 4 AY2 = ∑ ∑ (j - 1) ci j t i-1 u j-2 ; i=1 j=1 4 4 AY12 = ∑ ∑ (i - 1) (j - 1) ci j t i-2 u j-2 ; i=1 j=1
Здесь t = (X1 - X1L) / (X1U - X1L), U = (X2 - X2L) / (X2U - X2L). Коэффициенты ci j ( i, j = 1, 2, 3, 4 ) вычисляются в подпрограмме ip04r_c.
Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.
int ip04r_c (real *y, real *y1, real *y2, real *y12, real *x1l, real *x1u, real *x2l, real *x2u, real *x1, real *x2, real *ay, real *ay1, real *ay2, real *ay12, real *c)
Параметры
y, y1 - y2, y12 | вещественные векторы длины 4, содержащие соответственно значения интерполируемой функции, обеих ее частных производных первого порядка и смешанной производной в вершинах заданного элементарного сеточного прямоугольника, пронумерованных начиная с нижней левой вершины против часовой стрелки; |
x1l - x1u | левая и правая координаты вершин элементарного сеточного пpямoyгольника по первой переменной (тип: вещественный); |
x2l - x2u | нижняя и верхняя координаты вершин элементарного сеточного прямоугольника по второй переменной (тип: вещественный); |
x1, x2 - | координаты точки, заданной внутри исходного элементарного сеточного прямоугольника, в которой выполняется интерполяция (тип: вещественный); |
ay, ay1 - ay2, ay12 | вещественные переменные, значения которых полагаются равными соответственно значениям в точке (x1, x2) интерполируемой функции, ее частных производных первого порядка и смешанной производной; |
c - | вещественный двумерный массив размеров 4 на 4, в котором помещаются коэффициенты, используемые при бикубической интерполяции. |
Версии
ip04d_c - | бикубическая интерполяция функции двух переменных в режиме удвоенной точности; при этом все формальные параметры должны иметь тип double. |
Вызываемые подпрограммы: нет
Замечания по использованию: нет
int main(void) { /* Builtin functions */ double sin(double), cos(double); /* Local variables */ extern int ip04r_c(float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *, float *); static float c__[16] /* was [4][4] */, y[4], x1, y1[4], y2[4], x2, ay, y12[4], ty, ay1, ay2, x1l, x2l, x1u, x2u, ty1, ty2, ay12, ty12; int i; x1l = 0.f; x1u = .1f; x2l = 0.f; x2u = .15f; y[0] = (float)sin(x1l) * (float)sin(x2l); y[1] = (float)sin(x1u) * (float)sin(x2l); y[2] = (float)sin(x1u) * (float)sin(x2u); y[3] = (float)sin(x1l) * (float)sin(x2u); y1[0] = (float)cos(x1l) * (float)sin(x2l); y1[1] = (float)cos(x1u) * (float)sin(x2l); y1[2] = (float)cos(x1u) * (float)sin(x2u); y1[3] = (float)cos(x1l) * (float)sin(x2u); y2[0] = (float)sin(x1l) * (float)cos(x2l); y2[1] = (float)sin(x1u) * (float)cos(x2l); y2[2] = (float)sin(x1u) * (float)cos(x2u); y2[3] = (float)sin(x1l) * (float)cos(x2u); y12[0] = (float)cos(x1l) * (float)cos(x2l); y12[1] = (float)cos(x1u) * (float)cos(x2l); y12[2] = (float)cos(x1u) * (float)cos(x2u); y12[3] = (float)cos(x1l) * (float)cos(x2u); x1 = .06f; x2 = .07f; ip04r_c(y, y1, y2, y12, &x1l, &x1u, &x2l, &x2u, &x1, &x2, &ay, &ay1, &ay2, &ay12, c__); for (i = 0; i <= 12; i+= 4) { printf("\n %16.7e %16.7e %16.7e %16.7e \n", c__[i], c__[i+1], c__[i+2], c__[i+3]); } ty = (float)sin(x1) * (float)sin(x2); ty1 = (float)cos(x1) * (float)sin(x2); ty2 = (float)sin(x1) * (float)cos(x2); ty12 = (float)cos(x1) * (float)cos(x2); printf("\n %16.7e %16.7e \n",ay, ay1); printf("\n %16.7e %16.7e \n",ay2, ay12); printf("\n %16.7e %16.7e \n",ty, ty1); printf("\n %16.7e %16.7e \n",ty2, ty12); return 0; } /* main */ Результаты: ay = 0.4194o5e-02 ay1 = 0.698169e-01 ay2 = 0.598171e-01 ay12 = 0.995755