Partial Derivatives

函数式编程很好地解决了这个问题。的确， partial__ {xy}f是partial_y f沿x的偏导数。

我假设你有一个黑盒函数(或函数对象)f，将其值作为指向内存缓冲区的指针。它的签名被假定为

double f(double* x);

现在，这里是一个代码来获得f的(二阶有限差分)偏导数:

template <typename F>
struct partial_derivative
{
    partial_derivative(F f, size_t i) : f(f), index(i) {}
    double operator()(double* x)
    {
        // Please read a book on numerical analysis to tweak this one
        static const double eps = 1e-4;
        double old_xi = x[index];
        x[index] = old_xi + eps;
        double f_plus = f(x);
        // circumvent the fact that a + b - b != a in floating point arithmetic
        volatile actual_eps = x[index];
        x[index] = old_xi - eps;
        actual_2eps -= x[index]
        double f_minus = f(x);
        return (f_plus - f_minus) / actual_2eps;
    }
private:
    F f;
    size_t index;
};
template <typename F>
partial_derivative<F> partial(F f, index i)
{
    return partial_derivative<F>(f, i);
}

现在，要计算partial_{123}f，可以这样做:

boost::function<double(double*)> f_123 = partial(partial(partial(f, 0), 1), 2);

如果需要全部计算:

template <typename F>
boost::function<double(double*)> mixed_derivatives(F f, size_t* i, size_t n_indices)
{
    if (n_indices == 0) return f;
    else return partial(mixed_derivatives(f, i + 1, n_indices - 1), i[0]);
}

所以你可以这样做:

size_t indices[] = { 0, 1, 2 };
boost::function<double(double*)> df_123 
    = mixed_derivatives(f, indices, sizeof(indices) / sizeof(size_t));

如果我理解正确的话，我认为以下方法可以奏效:

function partial_dev(point, dimension):
    neighbor_selector = top_bit(dimension)
    value_selector = dimension XOR neighbor_selector
    prev_point = point_before(point,neighbor_selector)
    next_point = pointafter(point,neighbor_selector)
    if value_selector == 0:
        return (f[prev_point] - f[next_point])/2
    else:
        return ( partial_dev(prev_point, value_selector) -
                 partial_dev(next_point, value_selector) )/2

这个想法是:你的维度值的顶部位是在之前和之后点被选择的坐标。如果维度值的其余部分为0，则使用该点的f值进行偏导数计算。如果不是，则得到由其余位表示的偏导数来计算值。

如果需要计算所有维度值的所有值，则根本不需要递归:只需使用维度选择器作为数组索引，其中每个数组元素包含该维度的完整值集。初始化为vals[0][coords] = f(coords)。然后计算vals[1], vals[2]，计算vals[3]时，使用vals[1]作为值表而不是vals[0](因为3 = 0b11，其中邻居选择器为0b10, value_selector为0b01)。

看起来似乎可以根据维度(二进制位的数量)进行循环，然后递归到下一个二进制位。

粗略(非c++)伪代码:

Function partialcalc(leadingdigit, dimension)
  If dimension > 1 {
    For i = 1 to dimension {
      //do stuff with these two calls
      partialcalc(0, i - 1)
      partialcalc(1, i - 1)
    }
  }
  Else {
    //partialcalc = 1D case
  }
return partialcalc
End Function

递归的工作方式是，你有一个问题，它可以被分解成子问题，这些子问题与大问题等效，只是小一些。既然你用的是所有的二进制数字的第5位，那么你就可以在最上面的维度上做计算通过递归到两个子问题基于维度上的0和1的值。递归的底部是dimension = 1级别。既然您强调了您只需要弄清楚如何构造循环递归，并且已经弄清楚了数学运算，那么这个结构应该适合您。

好吧，从答案开始，我的第一个想法是用切比雪夫多项式进行插值。这个近似可以很容易地微分(或积分)。

Gnu Scientific Library有一个实现。

注意，我不是数值方法的专家，所以我不能解释这个方法的问题。然而，如果你想要一个局部近似值，它应该可以很好地工作。