Почему моя производная функция дает странные результаты?

Я пишу небольшую библиотеку исчисления для личного использования. В нем есть стандартные инструменты исчисления - получение первой производной, n-й производной, сумм Римана и т. д. Одна проблема, с которой я столкнулся, заключается в том, что функция n-й производной возвращает явно фиктивные результаты для определенных значений h (конечная разность).

Код здесь и ниже:

typedef double(*math_func)(double x);
inline double max ( double a, double b ) { return a > b ? a : b; }

//Uses the five-point-stencil algorithm.
double derivative(math_func f,double x){
    double h=max(sqrt(DBL_EPSILON)*x,1e-8);
    return ((-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h));
}

#define NDEPS (value)
double nthDerivative(math_func f, double x, int N){
    if(N<0) return NAN; //bogus value of N
    if(N==0) return f(x);
    if(N==1) return derivative(f,x);
    double* vals=calloc(2*N+9,sizeof(double)); //buffer region around the real values
    if(vals==NULL){ //oops! no memory
        return NAN;
    }
    int i,j;
    //don't take too small a finite difference
    double h=max(sqrt(DBL_EPSILON)*x,NDEPS);
    for(i=-(N+4);i<=N+4;i++){
        vals[i+N+4]=derivative(f,x+h*i);
    }
    //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    for(j=1;j<N;j++){
        double *vals2=calloc(2*N+9,sizeof(double));
        for(i=2;i<2*N+7;i++){
            vals2[i]=(-vals[i+2]+8*vals[i+1]-8*vals[i-1]+vals[i-2])/(12*h);
        }
        free(vals);
        vals=vals2;
        //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    }
    double result=vals[N+4];
    free(vals);
    return result;
}

Пример задачи, которую я даю для проверки этой функции, — это 5-я производная от sin(x), когда x=pi. Как мы все знаем, ответ равен -1. Проблема возникает, когда я пытаюсь изменить значение NDEPS («N-я производная эпсилон»).

Когда NDEPS=1,5625e-2 (1/64,0): 5-я производная от sin() при x=pi:-1,0003e+00 (хотя выглядит нормально).
Когда NDEPS=1e-5 (1/100000,0): 5-я производная от sin() при x=pi:2,4302e+11 (здесь я называю ерундой).

Почему это происходит? Связано ли это с природой функции sin()? Или это из-за проблем с точностью с плавающей запятой?

c floating-point calculus

o_o 09.10.2012 источник

comment

Вероятно, проблема с точностью с плавающей запятой. - Basile Starynkevitch 09.10.2012

comment

Хотя я не могу открыть вашу линию (заблокирована моим FW), я сильно подозреваю, что это из-за точности. - FamZheng 09.10.2012

comment

Поскольку мы не можем видеть вашу «идеально точную производную» функцию, мы не можем легко сказать, что происходит не так. Мы должны предположить, что это как-то связано с точностью с плавающей запятой. Я предполагаю, что мы можем использовать double derivative(math_func fun, double x) { return cos(x); } для работы с этим конкретным тестовым примером, но нам не нужно гадать, как вы написали код. Ознакомьтесь с тем, как написать SSCCE (краткий, автономный, правильный (компилируемый) пример). - Jonathan Leffler 09.10.2012

comment

@JonathanLeffler: я отредактировал вопрос, включив в него определение derivative(). - o_o 09.10.2012

Ответы (3)

arrow_upward
1
arrow_downward

Вот адаптация вашего кода. Основная модификация, помимо использования большего количества пробелов, состоит в том, чтобы сделать NDEPS аргументом функции nthDerivative(), чтобы ее можно было вызывать с разными значениями, и добавить обильную печать. Мне также пришлось проявить изобретательность в отношении простой функции derivative(); код компилируется правильно (но я действительно не пытаюсь делать философские или богословские утверждения с утверждением assert(sin == fun);, но это означает, что код компилируется без предупреждений и признает ограничения этой производной функции).

Код

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define max(a, b)    (((a) > (b)) ? (a) : (b))

#define PRIe_double    "%21.15e"

typedef double(*math_func)(double x);

static double derivative(math_func fun, double x) { assert(sin == fun); return cos(x); }

static double nthDerivative(math_func f, double x, int N, double NDEPS)
{
    if (N < 0) return NAN; //bogus value of N
    if (N == 0) return f(x);
    if (N == 1) return derivative(f, x);

    double* vals = calloc(2*N+9, sizeof(double)); //buffer region around the real values
    if (vals == NULL) //oops! no memory
        return NAN;

    int i, j;
    //don't take too small a finite difference

    double h = max(sqrt(DBL_EPSILON)*x, NDEPS);
    printf("h = " PRIe_double "\n", h);

    for (i = -(N+4); i <= N+4; i++)
    {
        vals[i+N+4] = derivative(f, x+h*i);
        printf("%2d: deriv(" PRIe_double ") = " PRIe_double "\n", i, x+h*i, vals[i+N+4]);
    }

    for (j = 1; j < N; j++)
    {
        printf("Iteration %d\n", j);
        double *vals2 = calloc(2*N+9, sizeof(double));
        for (i = 2; i < 2*N+7; i++){
            vals2[i] = (-vals[i+2] + 8*vals[i+1] - 8*vals[i-1] + vals[i-2]) / (12*h);
        }
        free(vals);
        vals = vals2;
        for (i = 0; i < 2*N+9; i++)
            printf("%2d: " PRIe_double "\n", i, vals[i]);
    }
    double result = vals[N+4];
    free(vals);
    return result;
}

int main(void)
{
    double val = M_PI;
    double eps;
    double r;

    eps = 1.0 / 64.0;
    r = nthDerivative(sin, val, 5, eps);
    printf("5th Derivative of sin(x) at x = " PRIe_double " = " PRIe_double " (eps = %f)\n", val, r, eps);

    eps = 1.0 / 100000.0;
    r = nthDerivative(sin, val, 5, eps);
    printf("5th Derivative of sin(x) at x = " PRIe_double " = " PRIe_double " (eps = %f)\n", val, r, eps);

    return(0);
}

Выход

При использовании GCC 4.7.1 в Mac OS X 10.7.5 вывод:

h = 1.562500000000000e-02
-9: deriv(3.000967653589793e+00) = -9.901285883701071e-01
-8: deriv(3.016592653589793e+00) = -9.921976672293290e-01
-7: deriv(3.032217653589793e+00) = -9.940245152582091e-01
-6: deriv(3.047842653589793e+00) = -9.956086864580017e-01
-5: deriv(3.063467653589793e+00) = -9.969497940760287e-01
-4: deriv(3.079092653589793e+00) = -9.980475107000991e-01
-3: deriv(3.094717653589793e+00) = -9.989015683384429e-01
-2: deriv(3.110342653589793e+00) = -9.995117584851364e-01
-1: deriv(3.125967653589793e+00) = -9.998779321710066e-01
 0: deriv(3.141592653589793e+00) = -1.000000000000000e+00
 1: deriv(3.157217653589793e+00) = -9.998779321710066e-01
 2: deriv(3.172842653589793e+00) = -9.995117584851364e-01
 3: deriv(3.188467653589793e+00) = -9.989015683384429e-01
 4: deriv(3.204092653589793e+00) = -9.980475107000991e-01
 5: deriv(3.219717653589793e+00) = -9.969497940760287e-01
 6: deriv(3.235342653589793e+00) = -9.956086864580017e-01
 7: deriv(3.250967653589793e+00) = -9.940245152582091e-01
 8: deriv(3.266592653589793e+00) = -9.921976672293291e-01
 9: deriv(3.282217653589793e+00) = -9.901285883701071e-01
Iteration 1
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -1.091570566584531e-01
 3: -9.361273104952932e-02
 4: -7.804555123490846e-02
 5: -6.245931771829009e-02
 6: -4.685783565504131e-02
 7: -3.124491392324735e-02
 8: -1.562436419383969e-02
 9: -1.776356839400250e-15
10: 1.562436419384087e-02
11: 3.124491392324558e-02
12: 4.685783565504250e-02
13: 6.245931771828653e-02
14: 7.804555123490846e-02
15: 9.361273104953050e-02
16: 1.091570566584471e-01
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 2
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -3.577900251527073e+00
 3: 1.660540592568783e+00
 4: 9.969497901146779e-01
 5: 9.980475067341610e-01
 6: 9.989015643694564e-01
 7: 9.995117545137316e-01
 8: 9.998779281977730e-01
 9: 9.999999960264082e-01
10: 9.998779281978486e-01
11: 9.995117545137319e-01
12: 9.989015643693869e-01
13: 9.980475067340947e-01
14: 9.969497901149178e-01
15: 1.660540592568512e+00
16: -3.577900251527123e+00
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 3
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 6.553266640232313e+01
 3: 1.898706817407991e+02
 4: -5.267598134705870e+01
 5: 3.608762837830827e+00
 6: 4.685783548517186e-02
 7: 3.124491378285654e-02
 8: 1.562436412277357e-02
 9: 3.226456139297321e-12
10: -1.562436412279785e-02
11: -3.124491378869069e-02
12: -4.685783548888563e-02
13: -3.608762837816180e+00
14: 5.267598134704988e+01
15: -1.898706817408119e+02
16: -6.553266640231030e+01
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 8.382087654791743e+03
 3: -5.062815705775390e+03
 4: -7.597917560836845e+03
 5: 3.261984801532626e+03
 6: -4.336626618856812e+02
 7: 1.791410702361821e+01
 8: -9.998779233550453e-01
 9: -9.999999914294619e-01
10: -9.998779238596159e-01
11: 1.791410702341709e+01
12: -4.336626618847606e+02
13: 3.261984801532444e+03
14: -7.597917560838102e+03
15: -5.062815705774387e+03
16: 8.382087654792240e+03
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = -9.999999914294619e-01 (eps = 0.015625)
h = 1.000000000000000e-05
-9: deriv(3.141502653589793e+00) = -9.999999959500000e-01
-8: deriv(3.141512653589793e+00) = -9.999999968000000e-01
-7: deriv(3.141522653589793e+00) = -9.999999975500000e-01
-6: deriv(3.141532653589793e+00) = -9.999999982000000e-01
-5: deriv(3.141542653589793e+00) = -9.999999987500000e-01
-4: deriv(3.141552653589793e+00) = -9.999999992000000e-01
-3: deriv(3.141562653589793e+00) = -9.999999995500000e-01
-2: deriv(3.141572653589793e+00) = -9.999999998000000e-01
-1: deriv(3.141582653589793e+00) = -9.999999999500000e-01
 0: deriv(3.141592653589793e+00) = -1.000000000000000e+00
 1: deriv(3.141602653589793e+00) = -9.999999999500000e-01
 2: deriv(3.141612653589793e+00) = -9.999999998000000e-01
 3: deriv(3.141622653589793e+00) = -9.999999995500000e-01
 4: deriv(3.141632653589793e+00) = -9.999999992000000e-01
 5: deriv(3.141642653589793e+00) = -9.999999987500000e-01
 6: deriv(3.141652653589793e+00) = -9.999999982000000e-01
 7: deriv(3.141662653589793e+00) = -9.999999975500000e-01
 8: deriv(3.141672653589793e+00) = -9.999999968000000e-01
 9: deriv(3.141682653589793e+00) = -9.999999959500000e-01
Iteration 1
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -7.000000116589669e-05
 3: -5.999999941330713e-05
 4: -5.000000598739025e-05
 5: -3.999999683331386e-05
 6: -2.999999600591015e-05
 7: -2.000000257999327e-05
 8: -1.000000082740371e-05
 9: 0.000000000000000e+00
10: 1.000000082740371e-05
11: 2.000000165480742e-05
12: 2.999999508072429e-05
13: 3.999999590812801e-05
14: 5.000000413701854e-05
15: 5.999999663774957e-05
16: 6.999999746515327e-05
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 2
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -3.583333244325556e+00
 3: 1.666666318844711e+00
 4: 1.000000128999663e+00
 5: 1.000000691821058e+00
 6: 9.999995738881511e-01
 7: 9.999997049561470e-01
 8: 1.000000198388602e+00
 9: 1.000000075030488e+00
10: 1.000000144419428e+00
11: 9.999996509869722e-01
12: 9.999995893079155e-01
13: 1.000000645561765e+00
14: 1.000000028771196e+00
15: 1.666666187776715e+00
16: -3.583333074708150e+00
17: 0.000000000000000e+00
18: 0.000000000000000e+00
Iteration 3
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 1.027777535146502e+05
 3: 2.972222191231725e+05
 4: -8.263881528669108e+04
 5: 5.555518108303881e+03
 6: -6.636923522614542e-02
 7: 4.677328483600658e-02
 8: 1.991719546327412e-02
 9: -3.148201866790915e-03
10: -2.319389535987426e-02
11: -4.176186143567406e-02
12: 6.726872146719150e-02
13: -5.555525175695851e+03
14: 8.263880834779718e+04
15: -2.972222015189417e+05
16: -1.027777456120210e+05
17: 0.000000000000000e+00
18: 0.000000000000000e+00

Признаки безумия...

Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: 2.050347140226725e+10
 3: -1.240740057099195e+10
 4: -1.858796490195886e+10
 5: 7.986101364079453e+09
 6: -1.059021715700324e+09
 7: 4.630176289959386e+07
 8: -3.687893612405425e+03
 9: -2.136279835945886e+03
10: -2.968840021291520e+03
11: 4.630204773690500e+07
12: -1.059022490436399e+09
13: 7.986100736580715e+09
14: -1.858796331554353e+10
15: -1.240739964045201e+10
16: 2.050347017082775e+10
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = -2.136279835945886e+03 (eps = 0.000010)

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

Использование «официальной» функции `derivative()`

Включение производной функции, которая теперь присутствует в вопросе, в приведенный выше код дает еще менее стабильные ответы на 4-й итерации:

Iteration 4
 0: 0.000000000000000e+00
 1: 0.000000000000000e+00
 2: -6.248925477935563e+09
 3: 1.729549900845405e+11
 4: -2.600544559219368e+11
 5: -2.755286326619338e+11
 6: 8.100546069731433e+11
 7: -2.961111189495415e+09
 8: -7.936480686806423e+11
 9: 2.430177384467434e+11
10: 2.389084910067162e+11
11: -6.461168564124718e+10
12: -4.574822745530297e+10
13: -8.923883451146609e+10
14: 7.042030613792160e+10
15: 4.988386306820556e+10
16: -1.793262395787471e+10
17: 0.000000000000000e+00
18: 0.000000000000000e+00
5th Derivative of sin(x) at x = 3.141592653589793e+00 = 2.430177384467434e+11 (eps = 0.000010)

Интересно, насколько появление нулей в индексах массива 0, 1, 17 и 18 усугубляет проблему.

Jonathan Leffler 09.10.2012

comment

Намерение выделить 2N+9 пробелов и получить реальное значение (N+4) состояло в том, чтобы значения медленно разрушались по направлению к центру, т.е. диапазон допустимых значений, окруженный двумя секциями бессмысленных значений, медленно смыкающихся... затем в в конце функции действительное значение уходит. (Или мне бы хотелось, чтобы это было именно так).‹br›Каждая итерация значения в позициях 0, 1, 17 и 18 остаются нетронутыми и, следовательно, равны 0. - o_o; 09.10.2012

arrow_upward
6
arrow_downward

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

f'(x) =(approx) (f(x+h)-f(x-h)) / (2*h)

это рецепт для отмены.

Это означает, что вы должны выбрать осторожный компромисс между двумя ошибками: если h слишком велико, вы получите большую ошибку аппроксимации в конечной разности. Если h слишком мало, вы получите большую (возможно, катастрофическую) числовую ошибку. Это хорошо показано в статье Википедии о численном дифференцировании. Проблемы усугубляются по мере увеличения порядка производной.

Вот почему автоматическое дифференцирование так важно — по сути, оно позволяет вам вычислить точную производную, когда все это задан алгоритм, вычисляющий базовую функцию. Если функция достаточно проста, чтобы вы могли символически определить производную, то вам, вероятно, следует это сделать.

Если вам нужно использовать численное дифференцирование, вы можете применить множество числовых хитростей. Числовые рецепты имеет хорошую трактовку — взгляните на раздел 5.7, который начинается на странице 186.

Martin B 09.10.2012

comment

Понятно... Спасибо за объяснение! - o_o; 09.10.2012

arrow_upward
1
arrow_downward

Это сложная тема. Численное дифференцирование плохо обусловлено и приводит к ошибкам сокращения и округления. Проблема значительно усугубляется в вашем случае, когда вы делаете повторную дифференциацию.

Вы получите гораздо лучшие оценки производных высокого порядка, определив коэффициенты для вашего конкретного порядка дифференцирования N с использованием m точек вместо повторного дифференцирования.

Например, так же, как вы, кажется, используете коэффициенты для производной первого порядка:

{1, -8, 8, -1} / (12*h)

коэффициенты для аппроксимации производной второго порядка с использованием 5 точек:

{-1, 16, -30, 16, -1} / (12*h²)

Вы можете найти коэффициенты для общей задачи, разложив ряд Тейлора по точкам выборки и найдя линейную комбинацию этих разложений, которая дает искомую производную.

для m точек {x + n[1]*h, x+n[2]*h, x+n[3]*h, …, x+n[m]*h} коэффициенты (k), оценивающие N-ю производную можно вычислить по следующей системе уравнений:

M*k=b

Где M - матрица m * m:

M[i,j] = n[j]^(i-1), i,j = 1..m

а также

b = [0, 0, 0, …, N!/h^N, …, 0, 0, 0],

где Н! обозначает факториал N.

digitalvision 09.10.2012

Почему моя производная функция дает странные результаты?

Ответы (3)

Код

Выход

Использование «официальной» функции derivative()

Похожие вопросы

Использование «официальной» функции `derivative()`