我正在尝试编写一个算法来执行 N 维混合偏导数。我知道我需要能够实现什么,但我似乎无法提出实现 N 维案例所需的正确循环/递归。
这是前 4 个维度的模式:
| 1D wzyx | 2D | 3D | 4D |
----------------------------------------------------------
| dx (0001) | dx (0001) | dx (0001) | dx (0001) |
| | dy (0010) | dy (0010) | dy (0010) |
| | dyx (0011) | dyx (0011) | dyx (0011) |
| | | dz (0100) | dz (0100) |
| | | dzx (0101) | dzx (0101) |
| | | dzy (0110) | dzy (0110) |
| | | dzyx (0111) | dzyx (0111) |
| | | | dw (1000) |
| | | | dwx (1001) |
| | | | dwy (1010) |
| | | | dwyx (1011) |
| | | | dwz (1100) |
| | | | dwzx (1101) |
| | | | dwzy (1110) |
| | | | dxyzw (1111) |
每个维度的导数个数(因为它遵循二进制模式)是(2^dim)-1;例如,2^3 = 8 - 1 = 7。
dyx的导数是y维度上相邻点的dx值。这适用于所有混合部分。所以 dzyx 是 z 维中相邻点的 dyx。我不确定这一段是否是该问题的相关信息,只是为了完整起见,我想把它放在这里。
欢迎任何帮助指针建议。粗体部分是我需要实现的部分。
::EDIT::
我将通过提供我需要的示例来尝试更明确一点。这只是一个 2D 案例,但它在某种程度上体现了我认为的整个过程。
我需要帮助想出将在 dx、dy、dyx 等列中生成值的算法。人。
| X | Y | f(x, y) | dx | dy | dyx |
-------------------------------------------------------------------------
| 0 | 0 | 4 | (3-4)/2 = -0.5 | (3-4)/2 | (-0.5 - (-2.0))/2 |
| 1 | 0 | 3 | (0-4)/2 = -2.0 | (2-3)/2 | (-2.0 - (-2.0))/2 |
| 2 | 0 | 0 | (0-3)/2 = -1.5 | (-1-0)/2 | (-1.5 - (-1.5))/2 |
| 0 | 1 | 3 | (2-3)/2 = -0.5 | (0-4)/2 | (-0.5 - (-0.5))/2 |
| 1 | 1 | 2 | (-1-3)/2 = -2.0 | (-1-3)/2 | (-1.5 - (-2.0))/2 |
| 2 | 1 | -1 | (-1-2)/2 = -1.5 | (-4-0)/2 | (-1.5 - (-1.5))/2 |
| 0 | 2 | 0 | (-1-0)/2 = -0.5 | (0-3)/2 | (-0.5 - (-0.5))/2 |
| 1 | 2 | -1 | (-4-0)/2 = -2.0 | (-1-2)/2 | (-2.0 - (-2.0))/2 |
| 2 | 2 | -4 |(-4--1)/2 = -1.5 |(-4--1)/2 | (-1.5 - (-1.5))/2 |
f(x, y) 是未知的,只有它的值是已知的;所以解析微分没用,只能是数字。
欢迎任何帮助指针建议。粗体部分是我需要实现的部分。
::EDIT - AGAIN::
在这里开始一个要点:https://gist.github.com/1195522
最佳答案
这个问题被函数式编程彻底解决了。实际上,\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));
关于c++ - 偏导数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/7304511/