3

私は個人的な使用のために小さな微積分ライブラリを書いています。その中には標準的な微積分ツールがあります-一次導関数、n次導関数、リーマン和などを取ります。私が直面した問題の1つは、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;
}

この関数をテストするために私が与えるサンプル問題は、x = piの場合のsin（x）の5次導関数です。ご存知のとおり、答えは-1です。問題は、NDEPS（「N次微分イプシロン」）の値を変更しようとすると発生します。

  • NDEPS = 1.5625e-2（1 / 64.0）の場合：x = pi：-1.0003e + 00でのsin（）の5次導関数（ただし、問題ないように見えます）。
  • NDEPS = 1e-5（1 / 100000.0）の場合：x = pi：2.4302e + 11でのsin（）の5次導関数（ここではでたらめと呼んでいます）。

なぜこれが起こるのですか？それはsin（）関数の性質に関連していますか？それとも、浮動小数点の精度の問題が原因ですか？

4

3 に答える 3

6

数値微分は難しい問題です。重要な問題は、有限差分近似です

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

キャンセルのレシピです。

これは、2つのエラーの間で慎重なトレードオフを選択する必要があることを意味します。hが大きすぎると、有限差分で大きな近似エラーが発生します。hが小さすぎると、大きな（おそらく壊滅的な）数値エラーが発生します。これは、数値微分に関するウィキペディアの記事で非常にうまく説明されています。導関数の次数が増えると、問題は悪化します。

これが、自動微分が非常に重要である理由です。本質的に、与えられたすべてが基本関数を計算するアルゴリズムである場合、正確な導関数を計算することができます。関数が十分に単純で、導関数を記号的に決定できる場合は、おそらくそれを行う必要があります。

数値微分を行う必要がある場合は、適用できる数値のフィネスがたくさんあります。数値レシピには良い扱いがあります。186ページから始まるセクション5.7を参照してください。

于 2012-10-09T07:30:40.993 に答える
1

これがあなたのコードの適応です。より多くの空白を使用することを除いて、主な変更は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);
}

出力

Mac OS X10.7.5でGCC4.7.1を使用すると、出力は次のようになります。

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でのゼロの出現は、問題をどの程度悪化させるのだろうか。

于 2012-10-09T06:53:03.120 に答える
1

これはトリッキーなテーマです。数値微分は悪条件であり、キャンセルと丸め誤差を引き起こします。あなたが繰り返し分化を行うあなたの場合、問題は著しく拡大します。

繰り返し微分を行う代わりに、m点を使用して特定の微分次数Nの係数を決定することにより、はるかに優れた高階微分推定が得られます。

たとえば、一次導関数の係数を使用しているように見えるのと同じように、次のようになります。

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

5点を使用して2階導関数を近似するための係数は次のとおりです。

{-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はam *m行列です。

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

と

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

ここでN！Nの階乗を示します。

于 2012-10-09T07:28:09.240 に答える