Время от времени в наших программах, при решении различных расчетных задач, возникает необходимость получать значения таблично заданных функций одной или двух переменных. Наиболее простым и надежным решением в данном случае является линейная интерполяция. В данной заметке я хочу привести пример определения значения таблично заданной функции двух переменных в виде реализации на C++.
В случаях, когда значения аргументов функции выходят за пределы таблицы, можно поступать по-разному: останавливать расчет и возвращать код ошибки или урезать значение аргумента до крайнего, представленного в таблице. В примере приведем первый вариант.
Стоит также учесть, что в зависимости от значений аргументов наш подход к работе с таблицей будет отличаться.
В случаях, когда значения аргументов функции выходят за пределы таблицы, можно поступать по-разному: останавливать расчет и возвращать код ошибки или урезать значение аргумента до крайнего, представленного в таблице. В примере приведем первый вариант.
Стоит также учесть, что в зависимости от значений аргументов наш подход к работе с таблицей будет отличаться.
- Значения обоих аргументов совпадают со значениями сетки. В этом случае, не прибегая ни к каким вычислениям, мы выбираем значение функции из таблицы.
- Один из аргументов функции совпал со значением сетки. Тогда следует воспользоваться формулой линейной интерполяции для функции одной переменной.
- Ни один из аргументов функции не совпал со значениями сетки. Используем в таком случае формулы линейной интерполяции для функции двух переменных, что, собственно, не сильно отличается от работы по п.2 - то же самое, только в два захода.
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} };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;
}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);
}
}