Scipy 的数值微分

Question

我试图学习 Scipy，将其用于混合集成和微分，但在最初的步骤中我遇到了以下问题。

对于数值微分，似乎唯一适用于可调用函数的 Scipy 函数是 scipy.derivative() 如果我是对的!？但是，我无法使用它:

1st) 当我不打算指定要进行微分的点时，例如当微分在积分下时，积分应该将数值分配给其被积函数的变量，而不是我。作为一个简单的例子，我在 Sage 的笔记本中尝试了这段代码:

import scipy as sp
from scipy import integrate, derivative
var('y')
f=lambda x: 10^10*sin(x)
g=lambda x,y: f(x+y^2)
I=integrate.quad( sp.derivative(f(y),y, dx=0.00001, n=1, order=7) , 0, pi)[0]; show(I)
show( integral(diff(f(y),y),y,0,1).n() )

它还发出警告“警告:检测到舍入误差的发生，这会阻止达到要求的公差。误差可能被低估了。”而且我不知道这个警告代表什么，因为即使增加“dx”和减少“order”它仍然存在。

2nd) 当我想在上面的例子中找到像 g(x,y) 这样的多变量函数的导数和像 sp.derivative(g(x,y),(x ,0.5), dx=0.01, n=1, order=3) 给出了错误，这是很容易预料到的。

期待收到您的来信，了解如何使用数值微分解决上述问题。最好的问候

Answer 1

您的代码存在一些奇怪的问题，这表明您需要温习一些 python!我不知道你是怎么在 python 中定义这些定义的，因为它们不是合法的语法。

首先，我认为您使用的是旧版本的 scipy。在最近的版本中(至少从 0.12+ 开始)你需要 from scipy.misc import derivative。 derivative 不在 scipy 全局命名空间中。

其次，var 未定义，尽管无论如何都没有必要(我认为您打算先导入 sympy 并使用 sympy.var('y')) . sin 也没有从 math(或 numpy，如果您愿意)中导入。 show 在 sympy 或 scipy 中不是有效函数。

^ 不是 python 中的幂运算符。你是说**

您在这里似乎混淆了符号和数值微积分运算的概念。 scipy 不会在数字上区分涉及符号对象的表达式——导数的第二个参数应该是你希望求导的点(即一个数字)。正如您所说，您正在尝试进行数值微分，我会为此目的解决问题。

from scipy import integrate
from scipy.misc import derivative
from math import *

f = lambda x: 10**10*sin(x)
df = lambda x: derivative(f, x, dx=0.00001, n=1, order=7)
I = integrate.quad( df, 0, pi)[0]

现在，最后一个表达式生成了您提到的警告，并且返回的值在绝对值 -0.0731642869874073 处不是非常接近于零，尽管相对于 f 的规模来说这还不错。您必须了解有限差分中的舍入误差问题。您的函数 f 在 0 到 10^10 之间的区间内变化!这可能看起来很矛盾，但是使微分的 dx 值太小实际上会放大舍入误差并导致数值不稳定。请参阅此处的第二张图(“显示由于舍入误差和公式误差而难以选择 h 的示例”)以获取解释:http://en.wikipedia.org/wiki/Numerical_differentiation

事实上，在这种情况下，您需要增加它，比如 0.001:df = lambda x: derivative(f, x, dx=0.001, n=1, order =7)

然后，您可以安全地集成，而不会出现可怕的舍入。

I=integrate.quad( df, 0, pi)[0]

我不建议丢弃 quad 的第二个返回值。这是对发生的事情的重要验证，因为它是“对结果绝对误差的估计”。在这种情况下，I == 0.0012846582250212652 并且绝对误差为 ~ 0.00022，这还不错(暗示的区间仍然不包括零)。也许对四边形的 dx 和绝对公差进行更多调整会给您带来更好的解决方案，但希望您能理解。

对于你的第二个问题，你只需要创建一个合适的标量函数(称之为 gx)来代表 g(x,y) 沿 y=0.5(这在计算机科学中被称为 Currying) .

g = lambda x, y: f(x+y**2)
gx = lambda x: g(x, 0.5)

derivative(gx, 0.2, dx=0.01, n=1, order=3)

为您提供 x=0.2 处的导数值。自然地，鉴于 f 的规模，其值(value)是巨大的。您可以像我上面展示的那样使用 quad 进行集成。

如果你想对 g 本身进行微分，你需要一个不同的数值微分函数。我认为 scipy 或 numpy 不支持这一点，尽管您可以通过制作 2D 精细网格(大小 dx)并使用 numpy.gradient 来破解中心差异计算。可能还有其他我不知道的库解决方案，但我知道我的 PyDSTool 软件包含一个函数 diff 可以做到这一点(如果你重写 g 以取一个数组参数)。它使用 Ridder 的方法，并受到 Numerical Recipes 伪代码的启发。