Частные производные

Ответ 1

Если вы правильно поняли, я думаю, что следующее может работать:

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).

Ответ 2

Эта проблема полностью решена с помощью функционального программирования. Действительно, \partial_ {xy} f является частной производной вдоль x of\partial_y f.

Я полагаю, у вас есть функция черного ящика (или объект функции) 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));

Ответ 3

Кажется, что вы можете просто создать цикл, основанный на измерении (количество двоичных мест), а затем вернуться к следующей двоичной цифре.

Грубый (не С++) псевдокод:

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

Как рекурсия работает, так это то, что у вас есть проблема, где ее можно разбить на подзадачи, эквивалентные большей проблеме, только меньше. Поэтому, поскольку вы используете все двоичные цифры в размерном месте, вы просто делаете вычисления в верхнем измерении, рекурсивным на две подзадачи на основе значений 0 и 1 в значении цифры. Нижняя часть рекурсии - это уровень размерности = 1. Поскольку вы подчеркиваете, что вам нужно только выяснить, как структурировать рекурсию цикла, и уже определили математику, эта структура должна работать на вас.

Ответ 4

Хорошо, чтобы начать с ответов, моя первая мысль заключалась в интерполяции с полиномами Чебышева. Тогда аппроксимация может быть легко дифференцирована (или интегрирована).

Научная библиотека Gnu имеет реализацию.

Заметьте, я не специалист по числовым методам, поэтому я не могу объяснить проблемы с этим методом. Однако он должен работать достаточно хорошо, если вы хотите локальное приближение. Если кто-то знает лучше, не стесняйтесь.