Текст подпрограммы и версий
ip03r_c.zip , ip03d_c.zip
Тексты тестовых примеров
tip03r_c.zip , tip03d_c.zip

Подпрограмма:  ip03r_c

Назначение

Полиномиальная интерполяция функции двух переменных, заданной на прямоугольной сетке, по модифицированной схеме Невилла.

Математическое описание

Пусть заданы узлы неравномерной прямоугольной сетки:

   X1A  =  { x1(1), x2(1), ..., xM(1) } ;    X2A  =  { x1(2), x2(2), ..., xN(2) }

 и значения  функции     f (x(1), x(2))

 в узлах этой сетки:     YA (I, J)  =  f (X1A(I), X2A(J)) ,

   I  =  1, 2, ..., M ;     J  =  1, 2, ..., N .

 Предполагается, что
    x1(1) < x2(1) < ... < xM(1)   и    x1(2) < x2(2) < ... < xN(2). 

Подпрограмма ip03r_c вычисляет значение Y = f (X1, X2) в заданной точке плоскости с координатами (X1, X2) и оценку ошибки DY полученного значения. При этом по горизонтальному направлению используется интерполяционный полином степени MP, а по вертикальному направлению - интерполяционный полином степени NP. Вычисление соответствующих значений интерполяционных полиномов осуществляется по модифицированной схеме Невилла.

Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.

Использование

    int ip03r_c (real *x1a, real *x2a, real *ya, integer *m,
            integer *n, integer *mp, integer *np, real *x1, real *x2, real *y, 
            real *dy, real *ym, real *yn, real *c, real *d, integer *ierr)

Параметры

x1a -
x2a  
вещественные векторы длины m и n, содержащие узлы  { x1 (1), x2 (1), ..., xm (1) } и { x1 (2), x2 (2), ..., xn (2) } соответственно;
ya - вещественный двумерный массив размеров m на n, содержащий значения интерполируемой функции двух переменных в узлах заданной сетки;
m, n - длины векторов x1a и x2a соответственно, m ≥ mp + 1, n ≥ np + 1 (тип: целый);
mp, np - заданные степени интерполяционных полиномов, mp ≥ 1, np ≥ 1 (тип: целый);
x1, x2 - координаты заданной точки на плоскости, в которой ищется значение интерполируемой функции (тип: вещественный);
y - вещественная переменная, содержащая значение интерполируемой функции в заданной точке плоскости;
dy - вещественная переменная, содержащая оценку ошибки вычисленного значения интерполируемой функции;
ym, yn - вещественные векторы длины m и n соответственно, используемые в подпрограмме в качестве рабочих;
c, d - вещественные векторы длины  max (mp + 1, np + 1), используемые в подпрограмме в качестве рабочих;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом
ierr=65 - когда по крайней мере два узла сетки по какому - либо направлению совпадают;
ierr=66 - когда число узлов сетки по какому - либо направлению меньше или равно степени соответствующего интерполяционного полинома.

Версии

ip03d_c - полиномиальная интерполяция функции двух переменных, заданной на прямоугольной неравномерной сетке, по модифицированной схеме Невилла в режиме удвоенной точности; при этом параметры x1a, x2a, ya, x1, x2, y, dy, ym, yn, c, d должны иметь тип double.

Вызываемые подпрограммы

ip02r_c - вычисление значения интерполяционного полинома в заданной точке по модифицированной схеме Невилла; используется в подпрограмме ip03r_c;
ip02d_c - вычисление значения интерполяционного полинома в заданной точке по модифицированной схеме Невилла в режиме удвоенной точности; используется в подпрограмме ip03d_c.

Замечания по использованию

  Заданная точка на плоскости не обязательно должна лежать внутри заданной прямоугольной сетки. В этом случае подпрограммы ip03r_c и ip03d_c выполняют полиномиальную экстраполяцию двух переменных.

Пример использования

int main(void)
{
    /* Builtin functions */
    double sin(double);

    /* Local variables */
    extern int ip03r_c(float *, float *, float *, int *, int *, int *,
                       int *, float *, float *, float *, float *,
                       float *, float *, float *, float *, int *);
    static int ierr;
    static float c__[6], d__[6];
    static int i__, j, m, n;
    static float r__, y, x1, x2, ya[150] /* was [15][10] */;
    static int mp, np;
    static float dy, ym[15], yn[10], x1a[15], x2a[10];
    int i__1, i__2;

#define ya_ref(a_1,a_2) ya[(a_2)*15 + a_1 - 16]

    m = 15;
    n = 10;
    r__ = 0.f;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        x1a[i__ - 1] = r__;
/* l1: */
        r__ += .1f;
    }
    r__ = 0.f;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        x2a[i__ - 1] = r__;
/* l2: */
        r__ += .1f;
    }
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
/* l3: */
            ya_ref(i__, j) = (float)sin(x1a[i__ - 1]) * (float)sin(x2a[j - 1]);
        }
    }
    mp = 5;
    np = 4;
    x1 = .55f;
    x2 = .65f;
    ip03r_c(x1a, x2a, ya, &m, &n, &mp, &np, &x1, &x2, &y, &dy, ym, yn, c__,
            d__, &ierr);

    printf("\n %16.7e %16.7e \n",y, dy);
    printf("\n %5i \n",ierr);
    return 0;
} /* main */


 Результаты:
     
             y = 0.316323 
          dy = 0.616967e-07 
       ierr = 0