3

我试图学习 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”并减少“订单”，它也会持续存在。

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

期待收到您关于如何解决上述数值微分问题的消息。此致

4

1 回答 1

0

您的代码存在一些奇怪的问题，表明您需要重新学习一些 python！我不知道你是如何在 python 中做出这些定义的，因为它们不是合法的语法。

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

其次，var未定义，尽管无论如何都没有必要（我认为您的意思是先导入 sympy 并使用sympy.var('y')）。sin也没有从数学（或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），它表示沿 y=0.5 的 g(x,y)（这在计算机科学中称为 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。您可以像我上面展示的那样使用 quad 进行集成。

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

于 2015-04-10T21:03:43.723 回答