Таблично заданные функции и линейная интерполяция на C++

Время от времени в наших программах, при решении различных расчетных задач, возникает необходимость получать значения таблично заданных функций одной или двух переменных. Наиболее простым и надежным решением в данном случае является линейная интерполяция. В данной заметке я хочу привести пример определения значения таблично заданной функции двух переменных в виде реализации на C++.
В случаях, когда значения аргументов функции выходят за пределы таблицы, можно поступать по-разному: останавливать расчет и возвращать код ошибки или урезать значение аргумента до крайнего, представленного в таблице. В примере приведем первый вариант.
Стоит также учесть, что в зависимости от значений аргументов наш подход к работе с таблицей будет отличаться.

1. Значения обоих аргументов совпадают со значениями сетки. В этом случае, не прибегая ни к каким вычислениям, мы выбираем значение функции из таблицы.
2. Один из аргументов функции совпал со значением сетки. Тогда следует воспользоваться формулой линейной интерполяции для функции одной переменной.
3. Ни один из аргументов функции не совпал со значениями сетки. Используем в таком случае формулы линейной интерполяции для функции двух переменных, что, собственно, не сильно отличается от работы по п.2 - то же самое, только в два захода.

Хорошо, теперь, наконец, переходим к коду. Пусть наша таблица определена несколькими векторами, содержащими значения сетки по оси X, значения сетки по оси Y и сами значения функции, то есть Z.

const vector<double> X = { 1, 1.5, 2, 2.5 };
const vector<double> Y = { 25, 30, 40 };
const vector< vector<double> > Z =
{ {0.392, 0.476, 0.574, 0.658},
  {0.406, 0.588, 0.770, 1.050},
  {0.434, 0.700, 0.980, 1.119} };

Затем приведем функции линейной интерполяции. В случае 3D, для пущей наглядности, я привожу в коде схему обозначения точек.

double linInterp2D(
    double x1,  double y1,
    double x2,  double y2,
    double x
    ) {

    double y = y1 + (x - x1) * (y2 - y1) / (x2 - x1);

    return y;
}

double linInterp3D(
    double x1,  double y1,
    double x2,  double y2,
    double z11, double z21, double z12, double z22,
    double x,   double y
    ) {

    /*

         |
      y2-| z12       z22
         |
       y-|      z
         |
      y1-| z11       z21
         |_______________
           |    |    |
           x1   x    x2

    */

    double z1 = z11 + (z12 - z11) * (y - y1) / (y2 - y1);
    double z2 = z21 + (z22 - z21) * (y - y1) / (y2 - y1);
    double z = z1 + (z2 - z1) * (x - x1) / (x2 - x1);

    return z;
}

И, наконец, нахождение значения таблично заданной функции по заданным x и y.

double findZ(double x, double y) {

    double x1 = 0; double x1i = 0;
    double y1 = 0; double y1i = 0;
    double x2 = 0; double x2i = 0;
    double y2 = 0; double y2i = 0;
    double z11 = 0;
    double z21 = 0;
    double z12 = 0;
    double z22 = 0;

    bool x_inside_range = false;
    bool y_inside_range = false;

    for (size_t i=0; i<X.size()-1; i++) {

        if ( x == X[i] ) {
            x1 = x2 = X[i];
            x1i = x2i = i;
            x_inside_range = true;
            break;
        }
        else if ( (x > X[i]) && (x < X[i+1]) ) {
            x1 = X[i];
            x1i = i;
            x2 = X[i+1];
            x2i = i+1;
            x_inside_range = true;
            break;
        }
        else if ( x == X[i+1] ) {
            x1 = x2 = X[i+1];
            x1i = x2i = i+1;
            x_inside_range = true;
            break;
        }
    }

    if ( !x_inside_range ) {
        cout << "x is outside of X axis range!\n";
        return -666;
    }

    for (size_t i=0; i<Y.size()-1; i++) {

        if ( y == Y[i] ) {
            y1 = y2 = Y[i];
            y1i = y2i = i;
            y_inside_range = true;
            break;
        }
        else if ( (y > Y[i]) && (y < Y[i+1]) ) {
            y1 = Y[i];
            y1i = i;
            y2 = Y[i+1];
            y2i = i+1;
            y_inside_range = true;
            break;
        }
        else if ( y == Y[i+1] ) {
            y1 = y2 = Y[i+1];
            y1i = y2i = i+1;
            y_inside_range = true;
            break;
        }
    }

    if ( !y_inside_range ) {
        cout << "y is outside of Y axis range!\n";
        return -666;
    }

    z11 = Z[y1i][x1i];
    z21 = Z[y1i][x2i];
    z12 = Z[y2i][x1i];
    z22 = Z[y2i][x2i];

    if ( (x1 == x2) && (y1 == y2) ) {
        return Z[y1i][x1i];
    }
    else if ( (x1 == x2) && (y1 != y2) ) {
        return linInterp2D(y1, z11, y2, z12, y);
    }
    else if ( (x1 != x2) && (y1 == y2) ) {
        return linInterp2D(x1, z11, x2, z21, x);
    }
    else {
        return linInterp3D(x1, y1, x2, y2, z11, z21, z12, z22, x, y);
    }
}

Вот и все. Вроде бы совершенно ничего сложного, но такие вот заготовочки весьма полезны, когда сроки поджимают, а некоторые функции уже написаны. Это здорово экономит время.