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

Question

Я пытаюсь написать алгоритм, который будет выполнять 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 - это значение dx смежных точек в измерении y. Это относится ко всем смешанным частям. Так что dzyx - это dyx смежных точек в измерении z. Я не уверен, что этот абзац является актуальной информацией для вопроса, просто подумал, что я бы поставил здесь для полноты.

Любые предложения указателей помощи приветствуются. Часть, выделенная жирным шрифтом, - это часть, которую мне нужно осознать.

:: EDIT ::

Я собираюсь попытаться быть более точным, приведя пример того, что мне нужно. Это всего лишь двухмерный случай, но я думаю, что он иллюстрирует весь процесс.

Мне нужна помощь в разработке алгоритма, который будет генерировать значения в столбцах dx, dy, dyx, et. и др.

|  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) неизвестно, известны только его значения; поэтому аналитическая дифференциация бесполезна, она должна быть только числовой.

Любые предложения по указателям помощи приветствуются. Часть, выделенная жирным шрифтом, - это часть, которую мне нужно осознать.

:: РЕДАКТИРОВАТЬ - СНОВА ::

Начал Gist здесь: https://gist.github.com/1195522

Answer 1

Эта проблема четко решена с помощью функционального программирования.В самом деле, \ part_ {xy} f является частной производной по x от \ part_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);
}

Теперь для вычисления \ part__{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));

Answer 2

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

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, а селектор значения -0b01).

Answer 3

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

Грубый (не 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

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

Answer 4

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

В Научной библиотеке Gnu есть реализация.

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