文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
首页
学习
活动
专区
圈层
工具
MCP广场
返回腾讯云官网
社区首页 >问答首页 >为什么我的导数函数会产生奇怪的结果？

问为什么我的导数函数会产生奇怪的结果？
EN

Stack Overflow用户
提问于 2012-10-09 05:22:56
回答 3查看 465关注 0票数 3

我正在写一个小型微积分库，供个人使用。它是标准的微积分工具--取一阶导数、n阶导数、Riemann和等。我所面临的一个问题是，第n阶导数函数对h(有限差分)的某些值返回明显的假结果。

代码在这里及以下：

代码语言：javascript
运行
复制
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)的五阶导数。我们都知道，答案是-1。当我试图改变NDEPS (“N个导数epsilon")的值时，问题就出现了。

  • 当NDEPS=1.5625e-2 (1/64.0)：sin()在x=pi:-1.0003e+00处的五阶导数时(看起来不错)。
  • 当NDEPS=1e-5 (1/100000.0)：sin()在x=pi:2.4302e+11的五阶导数时(我在这里称之为扯淡)。

这一切为什么要发生？它与sin()函数的性质有关吗？还是因为浮点精度问题？

c
floating-point
calculus
EN

回答 3

Stack Overflow用户

回答已采纳

发布于 2012-10-09 06:53:03

下面是您代码的改编本。除了使用更多的空白之外，主要的修改是将NDEPS变成nthDerivative()函数的一个参数，这样就可以用不同的值调用它，并添加大量的打印。我还必须对普通的derivative()函数有创造性；代码编译正确(但我真的不想用断言assert(sin == fun);来做哲学或神学上的陈述，但它确实意味着代码编译时没有警告，它承认这个派生函数的局限性)。

代码

代码语言：javascript
运行
复制
#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);
}

输出

在MacOSX10.7.5上使用GCC 4.7.1，输出如下：

代码语言：javascript
运行
复制
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

疯狂的迹象..。

代码语言：javascript
运行
复制
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)

注意，在上一次迭代中，结果是如何被数十亿个值所困扰的。你至少有一个数值稳定性问题，或者你需要复习一下导数公式。请注意，即使是使用较大epsilon的运行，在以后的迭代中也会出现大值的倾向。

使用“正式”derivative()函数

将问题中现在出现的导数函数插入上面的代码，在第4次迭代中得到更不稳定的答案：

代码语言：javascript
运行
复制
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上的零的出现在多大程度上加剧了这个问题。

票数 1
EN

Stack Overflow用户

发布于 2012-10-09 07:30:40

数值微分是一个难题。关键问题是有限差分逼近

代码语言：javascript
运行
复制
f'(x) =(approx) (f(x+h)-f(x-h)) / (2*h)

是取消的配方。

这意味着您必须在两个错误之间进行仔细的权衡:如果h太大，那么在有限差分中会产生一个很大的近似误差。如果h太小，您将产生一个很大(可能是灾难性的)数值错误。维基百科关于数值微分的文章很好地说明了这一点。这些问题随着衍生产品数量的增加而加剧。

这就是为什么自动微分如此重要的原因--本质上，它允许你计算一个精确的导数，当给出的只是一个计算基函数的算法时。如果函数足够简单，可以象征性地确定导数，那么很可能应该这样做。

如果你确实需要与数值微分，有很多的数字技巧，你可以应用。数值累加有一个很好的待遇--看看5.7节，它从第186页开始。

票数 6
EN

Stack Overflow用户

发布于 2012-10-09 07:28:09

这是个棘手的话题。数值微分是病态的，会产生抵消和舍入误差.在你的情况下，这个问题被大大放大了，在这里，你要做重复的区分。

通过使用m点来确定特定微分阶N的系数，而不是重复微分，你将得到更好的高阶导数估计。

例如，就像您在一阶导数中使用系数一样：

代码语言：javascript
运行
复制
{1, -8, 8, -1} / (12*h)

用5点逼近二阶导数的系数是：

代码语言：javascript
运行
复制
{-1, 16, -30, 16, -1} / (12*h²)

你可以通过在样本点周围做泰勒级数展开来解决一般问题的系数，找出这些展开的线性组合，给出想要的导数。

对于m点{x + n1*h，x+n2*h，x+n3*h，…，x+nm*h}，估计N阶导数的系数(k)可由下列方程组计算：

代码语言：javascript
运行
复制
M*k=b

其中M是m*m矩阵：

代码语言：javascript
运行
复制
M[i,j] = n[j]^(i-1), i,j = 1..m

和

代码语言：javascript
运行
复制
b = [0, 0, 0, …, N!/h^N, …, 0, 0, 0],

哪里N！表示N的阶乘。

票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/12793280

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2025 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档