python中的双反导数计算

Question

我有以下问题。我有一个函数f在python中使用Numpy函数定义。函数平滑，可积于正实。我想构造函数的双反导数（假设0处的反导数的值和斜率都是0），这样我就可以在任何小于100的正实数上计算它。

f在x的反导数定义：

integrate f(s) with s from 0 to x

fatx的双重反导数的定义：

integrate (integrate f(t) with t from 0 to s) with s from 0 to x

f的实际形式并不重要，所以为了方便起见，我将使用一个简单的。但是请注意，即使我的例子有一个已知的封闭形式，我的实际函数却没有。

import numpy as np

f = lambda x: np.exp(-x)*x

我的解决方案是使用简单的数值积分将反导数构造为数组：

N = 10000
delta = 100/N

xs = np.linspace(0,100,N+1)

vs = f(xs)
avs = np.cumsum(vs)*delta
aavs = np.cumsum(avs)*delta

这当然有效，但它给了我数组而不是函数。但这不是一个大问题，因为我可以使用样条插值aavs，得到一个函数并去掉数组。

from scipy.interpolate import UnivariateSpline

aaf = UnivariateSpline(xs, aavs)

函数aaf近似为f的双反导数。

问题是，即使它工作，有相当多的开销之前，我可以得到我的函数和精度是昂贵的。

我的另一个想法是用样条插值f，并取其反导数，但是这会引入数值误差，对于我想要使用的函数来说，数值误差太大了。

有没有更好的办法？我所说的更好是指在不牺牲准确性的情况下更快。

编辑：我希望可以使用某种傅里叶变换来避免两次积分。我希望有一些方便的vs转换，允许将组件的值与xs相乘，并转换回双反导数。我玩了一会儿，但我迷路了。

编辑：我发现通过使用梯形规则而不是简单的求和，可以大大提高准确性。使用Simpson的规则可以进一步提高精度，但是对于Numpy数组来说有点麻烦。

编辑：正如@user202729正确地抱怨的那样，这似乎是错误的。它看起来不对劲的原因是我跳过了一些细节。我在这里解释为什么我说的有意义，但它并不影响我的问题。

我的实际目标不是找到f的双重反导数，而是找到它的一个变换。我跳过了这一点，因为我认为这只会混淆问题。

函数f在x接近0或无穷大时呈指数衰减。我通过从0开始求和并上升到大约f的峰值来最小化积分中的数值误差。这可确保相对误差近似恒定。然后我从一个非常大的x的相反方向开始，回到峰值。然后我对反导数值做同样的处理。

然后我用另一个对数字错误敏感的函数变换aavs。然后我找到误差较大的区域（值剧烈振荡）并删除这些值。最后，我用样条近似我认为是好的值。

现在，如果我使用样条曲线来近似f，它引入了一个绝对误差，这是一个相当大的区间内的主导项。这两次被“集成”，最终在aavs中成为一个相当大的相对错误。然后，一旦我转换了AAV，我发现“好的区域”已经大大缩小了。

Answer 1

首先，我创建了一个更直观的代码版本。这里，我将累积和值乘以箱子宽度。我相信与箱子宽度问题相关的原始版本代码中有一个小错误。

import numpy as np
f = lambda x: np.exp(-x)*x
N = 1000
xs = np.linspace(0,100,N+1)

domainwidth = ( np.max(xs) - np.min(xs) )
binwidth = domainwidth / N
vs = f(xs)
avs = np.cumsum(vs)*binwidth
aavs = np.cumsum(avs)*binwidth

接下来，为了可视化，这里有一些非常简单的绘图代码：

import matplotlib
import matplotlib.pyplot as plt
plt.figure()
plt.scatter( xs, vs )
plt.figure()
plt.scatter( xs, avs )
plt.figure()
plt.scatter( xs, aavs )
plt.show()

第一个积分与示例表达式的已知结果相匹配，可以在Wolfram上看到

下面是一个从二阶导数中提取元素的简单函数。请注意，int是一个错误的舍入函数。我想这就是你已经实现的。

def extract_double_antideriv_value(x):
    return aavs[int(x/binwidth)]
singleresult = extract_double_antideriv_value(50.24)
print('singleresult', singleresult)

无论需要什么完整的计算步骤，我们都需要在开始优化之前了解它们。你有一百万种不同的功能需要集成吗？如果您只需要多次查询单个双反导数，那么您的原始解决方案应该是相当理想的。