我正在寻找一种使用离散快速方法计算导数的方法。因为目前我不知道我所拥有的方程类型,所以我正在寻找类似于我们可以在积分中找到的欧拉方法等的离散方法。
我正在寻找一种使用离散快速方法计算导数的方法。因为目前我不知道我所拥有的方程类型,所以我正在寻找类似于我们可以在积分中找到的欧拉方法等的离散方法。
您实际上需要实现极限函数。所以你需要:
现在在DO-WHILE循环中:
1- 将h除以2(或10),重要的是要使它更小
2- 使用新的h值重新计算差商,将其存储在f2中
3- 设置diff = abs(f2-f1)
4- 赋值f1 = f2
5- 重复从第1点开始 while (diff>epsilon)
请记住: 您假设函数在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的大小选择取决于您希望结果有多精确。
在计算数值(“有限”)导数方面,有相当多的理论(和已经建立的实践)。要正确地获取所有细节,以使您相信结果,这并不是简单的。如果您可以通过笔和纸或计算机代数系统(如Maple、Mathematica、Sage或SymPy)获得函数的解析导数,那么这是远远最好的选择。
如果您无法获得解析形式,或者您不知道函数(只知道它的输出),那么数值估计是您唯一的选择。Numerical Recipies in C中的本章是一个很好的开始。
一个简单的方法是计算每个导数点上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的内容、所需精度、浮点类型以及可能的其他因素。我建议您使用您感兴趣的函数尝试不同的值。
* Andrea的回答使用前向差分算子{f(a+h) - f(a)}/h,而不是中心差分算子{f(a+h) - f(a-h)}/2h,但在数值解中前向差分算子的精度较低。
从正式意义上讲,不行。你要么在描述离散函数的(部分)导数,要么是在寻求一种近似连续函数(部分)导数的数值方法。
离散函数没有导数。如果你回顾一下导数的 epsilon-delta 定义,你会发现你需要能够在你想要求导的点附近计算函数的值。如果函数只在 x、y 和 z 的整数值处有值,那么这是没有意义的。因此,无论快速的取什么值,都无法找到离散函数的导数。
如果你想要一个数值方法来精确计算连续函数的导数,那么你也会失望。导数的数值方法是启发式的,而不是算法性的。没有一种数值方法可以保证得到精确的解决方案。幸运的是,存在许多好的启发式方法。Mathematica默认使用Brent's principle axis method的专门版本。我建议你使用GNU Scientific Library,它有一个非常好的Brent方法实现。我的一个数学课程的整个成绩都归功于GSL。如果需要,大多数数值微分库都有几种不同的方法可用。
如果你真的想要,我可以给你一些示例代码。让我知道。