将振荡(插值)函数与GSL和python集成在一起
我试图将一个高度振荡的数据与python中的GSL科学图书馆和吡咯烷酮中的qawo函数集成起来。由于我正在处理数据，所以我认为插值函数可以工作，但是GSL给出的结果是不正确的-- !!。让我以函数Sin(x)/(1+x²)为例进行解释。
以下代码工作正常：
import pygsl
from pygsl import integrate
def f1(x,y):
return 1./(1 + x**2)
sys = integrate.gsl_function(f1, None)
w = integrate.workspace(100)
cyclew = integrate.workspace(1000000)
table = integrate.qawo_table(1, 10000, integrate.SINE, 100)
flag, result, error = integrate.qawo(sys, 0, 1e-8, 1e-8, 100, w, table)为0.626761提供应有的。但是如果我们用上面的函数模拟数据点..。
xarr = np.linspace(0,1e15,1e30)
yarr = np.sin(xarr)/(1.+xarr**2)
interp = interpol.interp1d(xarr,yarr)
def fgsl(x,y):
return interp(x)
syst = integrate.gsl_function(fgsl, None)
w = integrate.workspace(1000)
cyclew = integrate.workspace(100000000)
table = integrate.qawo_table(1, 1e10, integrate.SINE, 100)
flag, result, error = integrate.qawo(syst, 0, 1e-15, 1e-15, 100, w, table)给出了一个完全错误的结果：4.2426e-21
此外，如果我们将yarr与simps功能集成在一起：
import scipy.integrate as ints
res = ints.simps(yarr,xarr)给出了一个很好的近似：0.64676099。
假设我不能用辛普森法则。有人知道如何在gsl中使用插值函数吗？或者如何更改代码以完成集成？
提前谢谢。
回答 1
Stack Overflow用户
发布于 2015-01-25 19:14:01
你的例子中的数字不太合理，它们会破坏任何自适应方案。让我解释一下原因。
您正在尝试集成一个周期为2*Pi从0到10^10的振荡函数！没有任何自适应方案能够“看到”该区间上的振荡行为，并且它们将收敛到w/错误的结果(错误收敛)！请记住，自适应方案采用自顶向下的方法。他们在整个区间上应用了一些规则，然后将该区间分成两部分，然后在每个细分中应用相同的规则。经过几个循环(通常是4或5)后，该方案通过比较连续步骤的部分结果开始检查收敛性。在您的示例中，该方案将需要很多细分才能最终看到振荡行为，这是一个典型的情况，可能发生错误收敛！
如何在开放区间(a,\infinity)上集成振荡函数？对qawf积分集成方案的解释是相当完整的。在只包含几个振荡的子区间上集成函数，并检查结果是如何收敛的，然后进行推断！
还有其他的数字不太合理。为什么您需要在每个dx=1e-15处取样sin(x)/(1+x^2)？任何合理的自适应方案都可以将sin(x)从0到2piw/ ~10-20样本点进行集成。
辛普森的规则并没有失败，因为它不是一个适应性的方案。python代码将根据您提供的x数组确定' dx‘，它将一直使用该dx到1e10！但是，我非常肯定，在代码中舍入错误是非常糟糕的，因为您选择了dx~1e-15。
编辑1第一部分:实际上这个问题不仅仅是由被积体的振荡行为引起的。假设包络1/x^2收敛很快--如果是x>>1的话，你的函数实际上是零。所以，由于你在巨大的间隔0,1e10中积分这个包络，自适应积分认为结果很小，因为它不能看到函数不可忽略的小(子)区间。(您可能认为，积分例程将在接近区间0,1e10中均匀地分配评价点-这对于高斯积分不是完全正确的，但它是接近的--因此，在积分不可忽略的区间~0,1e3范围内，其中一个点的概率很小。经过几次迭代之后，集成例程将得到积分接近于零的结果)。
编辑1第二部分:我仍然认为(在阅读了您的评论之后)问题在于您插入的数字(或者python包装器有一个bug。)请使用合理的数字尝试您的示例，就像我在下面的C++代码中所做的那样
int main()
{
const double omega = 1;
auto g = [](double x)->double {return 1.0/(1.+x*x);};
auto f = [&](double x)->double {return std::sin( omega * x) * g(x);};
const int points_per_cycle = 20;
const int n_cycles = 10;
const int size = n_cycles * points_per_cycle + 1;
const double xmin = 0.0;
const double xmax = n_cycles * (2 * M_PI);
const double L = xmax-xmin;
std::vector<double> x(size);
std::vector<double> y(size);
for (int i = 0; i <size; ++i) {
x[i] = i * L/(size-1);
y[i] = f(x[i]);
}
std::cout.precision(8);
// interpolation
InterpolationGSL<std::vector<double>> GSLinterpol(x, y, GSLIT::cspline, false);
// Integral of the interpolation
std::cout << GSLinterpol.If((1+1e-12)*xmin, (1-1e-12)*xmax) << std::endl;
// SECOND GSL INTEGRATION
gsl_integration_workspace* w = gsl_integration_workspace_alloc (1000);
gsl_integration_qawo_table* wf = gsl_integration_qawo_table_alloc
(omega, L, GSL_INTEG_SINE, 1000);
int status = gsl_integration_qawo_table_set (wf, omega, L, GSL_INTEG_SINE);
if(status) std::cerr<< "error: " << std::string(gsl_strerror (status)) << std::endl;
double result;
double abserr;
std::function<double(double)> gg( std::cref(g) );
GslFunction Fp(gg);
gsl_function *Fgsl = static_cast<gsl_function*>(&Fp);
status = gsl_integration_qawo (Fgsl, xmin, 0.0, 1e-5, 1000, w, wf, &result, &abserr);
if(status) std::cerr<< "error: " << std::string(gsl_strerror (status)) << std::endl;
std::cout << result << std::endl;
return 0;
}这段代码使用了我的函数和内插包装器--所以您可能觉得代码有点奇怪--但重要的一点是，它在合理的间隔内计算您提到的相同的积分，结果如下
Interpolation integral: 0.64631754
GSL integral: 0.64650827https://stackoverflow.com/questions/28131328
复制相似问题
腾讯云开发者
