使用离散方法计算导数

12

我正在寻找一种使用离散快速方法计算导数的方法。因为目前我不知道我所拥有的方程类型,所以我正在寻找类似于我们可以在积分中找到的欧拉方法等的离散方法。

10个回答

11
我认为您正在寻找在某一点计算的导数。 如果是这种情况,那么有一种简单的方法。您需要知道在一个点,比如说a,导数的值。它由差商极限给出,当h->0时:

difference quotient

您实际上需要实现极限函数。所以你需要:

  • 定义一个epsilon,设置它更小来提高精度,更大来提高速度
  • 在起始h中计算差商,假设h=0.01,将其存储在f1
  • 现在在DO-WHILE循环中:

    1- 将h除以2(或10),重要的是要使它更小
    2- 使用新的h值重新计算差商,将其存储在f2
    3- 设置diff = abs(f2-f1)
    4- 赋值f1 = f2
    5- 重复从第1点开始 while (diff>epsilon)

  • 最后您可以将f1(或f2)作为f'(a)的值返回

请记住: 您假设函数在a处是可微的。 由于您的计算机只能处理有限的十进制数字,因此您得到的每个结果都会出现错误。

Python示例:

def derive(f, a, h=0.01, epsilon = 1e-7):
    f1 = (f(a+h)-f(a))/h
    while True: # DO-WHILE
        h /= 2.
        f2 = (f(a+h)-f(a))/h
        diff = abs(f2-f1)
        f1 = f2
        if diff<epsilon: break
    return f2

print "derivatives in x=0"
print "x^2: \t\t %.6f" % derive(lambda x: x**2,0)
print "x:\t\t %.6f" % derive(lambda x: x,0)
print "(x-1)^2:\t %.6f" % derive(lambda x: (x-1)**2,0)

print "\n\nReal values:"
print derive(lambda x: x**2,0)
print derive(lambda x: x,0)
print derive(lambda x: (x-1)**2,0)

输出:

derivatives in x=0
x^2:         0.000000
x:       1.000000
(x-1)^2:     -2.000000


Real values:
7.62939453125e-08
1.0
-1.99999992328

由于只使用了结果的前6位数字,我第一次得到了“精确”值。请注意,我使用了1e-7作为epsilon。之后打印了真实计算出的值,它们显然在数学上是错误的。epsilon的大小选择取决于您希望结果有多精确。


7

在计算数值(“有限”)导数方面,有相当多的理论(和已经建立的实践)。要正确地获取所有细节,以使您相信结果,这并不是简单的。如果您可以通过笔和纸或计算机代数系统(如MapleMathematicaSageSymPy)获得函数的解析导数,那么这是远远最好的选择。

如果您无法获得解析形式,或者您不知道函数(只知道它的输出),那么数值估计是您唯一的选择。Numerical Recipies in C中的本章是一个很好的开始。


请记住,Numerical Recipes 中的代码不是公有领域或开源的。复制需自负风险。 - Joseph Holsten
1
是的,绝对不会计划从NR in C(或任何其他受版权保护的来源)窃取代码。但是,您可以通过阅读该章节学到足够的知识,以便编写自己的代码 ;) - Barry Wark

4

一个简单的方法是计算每个导数点上f值的微小变化。例如,要计算∂f/∂x,您可以使用以下公式:

epsilon = 1e-8
∂f/∂x(x, y, z) = (f(x+epsilon,y,z) - f(x-epsilon, y, z))/(epsilon * 2);

其他的偏导数在y和z方向上也类似。

epsilon的选择取决于f的内容、所需精度、浮点类型以及可能的其他因素。我建议您使用您感兴趣的函数尝试不同的值。


除非你真的想深入了解你正在使用的语言的浮点实现,否则很难选择一个好的 epsilon。但请发布您发现的奇怪浮点错误! - Joseph Holsten

4
要进行数值微分,它总是一个近似值,有两种常见情况:
1. 您有一种方法(算法、方程)可以计算任何给定x的f(x)的值;或者
2. 您有f(x)在一组等间距的x值上的值(f(1),f(1.5),f(2)等)。



  1. 如果您已经有了算法,那么最好参考Andrea Ambu的答案*:
    1. 从f(a+h)和f(a-h)的值开始,并执行
      f'(a) = { f(a+h)-f(a-h) } / 2h
      这被称为中心差分。
    2. 缩小偏移量h的大小,直到f'(a)停止变化,如Andrea的答案所述。
    • 请注意,在没有检查它是否大于0的情况下除以2h,您将获得除以零错误。
  2. 如果您有一组f(x)=f(xn)=fn的函数值,则我们可以使用与上述方法类似的方法,但我们的精度将受到xn的值和它们之间的间隔的限制。
  3. 第一个近似只是使用与上面相同的中心差分算子;
    f'n = (fn+1 - fn-1)/2h,
    其中h是xn值之间的等间距。
  4. 接下来使用五点差分, 尽管它只有4个系数,如维基百科页面上详细说明的那样。这使用fn-2到fn+2的值:f'n = (-fn+2 + 8fn+1 - 8fn-1 +fn-2)/12h。
  5. 有更多点的高阶近似,但对于递减回报而言,它们变得越来越难以计算,并在数值上变得更加不稳定。
  6. **注意**:这些有限差分公式依赖于f在xn-1≤ x ≤n+1范围内大致与多项式相同的形状。对于正弦波之类的函数,它们在计算导数时可能会非常糟糕。

* Andrea的回答使用前向差分算子{f(a+h) - f(a)}/h,而不是中心差分算子{f(a+h) - f(a-h)}/2h,但在数值解中前向差分算子的精度较低。


2
除了使用类似Maple的符号数学语言,你最好能够在各个点上近似计算导数。(如果需要一个函数,那么就插值。)
如果你已经有了想要使用的函数,那么你应该使用向后差商公式Richardson外推法来改善你的误差。
同时请注意,这些方法适用于单变量函数。然而,每个变量的偏导数将其他变量视为常数。

理查逊外推法很有用,但为什么要使用向后的分裂差分而不是中心的呢? - Jesse Rusak
这些似乎是两个变量方法,不像问题中的多元方法。我是否漏掉了什么?如果没有,我不能建议将一个双变量方法改为更多变量的方法。这真的容易出错。 - Joseph Holsten

2

自动微分是最准确、概念最棒的方式来完成这种任务,只是稍微复杂一些。


1

从正式意义上讲,不行。你要么在描述离散函数的(部分)导数,要么是在寻求一种近似连续函数(部分)导数的数值方法。

离散函数没有导数。如果你回顾一下导数的 epsilon-delta 定义,你会发现你需要能够在你想要求导的点附近计算函数的值。如果函数只在 x、y 和 z 的整数值处有值,那么这是没有意义的。因此,无论快速的取什么值,都无法找到离散函数的导数。

如果你想要一个数值方法来精确计算连续函数的导数,那么你也会失望。导数的数值方法是启发式的,而不是算法性的。没有一种数值方法可以保证得到精确的解决方案。幸运的是,存在许多好的启发式方法。Mathematica默认使用Brent's principle axis method的专门版本。我建议你使用GNU Scientific Library,它有一个非常好的Brent方法实现。我的一个数学课程的整个成绩都归功于GSL。如果需要,大多数数值微分库都有几种不同的方法可用。

如果你真的想要,我可以给你一些示例代码。让我知道。


非常感谢您提供的所有信息,我希望能找到一些简单的东西,但似乎对于我的问题并没有简单的解决方案。我将采用任何不使用导数的其他方法。再次感谢您。 - Ryan

0
我会假设你的函数比你发布的简单函数更复杂,因为闭式解法太过简单。
当你使用“离散”这个词时,让我想到你需要“有限差分”。你需要一些离散化来计算近似值。
Df/Dx ~ (f2-f1)/(x2-x1),等等。

0

-1
如果您所指的函数是线性的,那么导数就很简单了。对于'x'的导数是'a';对于'y'的导数是'b';对于'z'的导数是'c'。如果方程式更为复杂,您需要一个代表解决方案的公式而不是经验解决方案,请提交更复杂的方程式。
此致敬礼

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接