计算科学 - 给出不正确符号的数值积分 - 吾爱随笔录

给出不正确符号的数值积分

计算科学一体化正交 C

2021-12-19 18:01:31

对于我的研究，我需要集成以下功能：

W (z) = \int_{0}^{\infty} d x w (x, z) = \int_{0}^{\infty} d x \frac{e^{x}}{(e^{x} + 1)^{2}} \log (\frac{e^{z^{2} / 4 x + x} + 1}{e^{z^{2} / 4 x + x} - e^{x}}) = \int_{0}^{\infty} d x \frac{e^{x}}{(e^{x} + 1)^{2}} \log (\frac{e^{z^{2} / 4 x} + e^{- x}}{e^{z^{2} / 4 x} - 1})

$W(z)=\int_0^{\infty}dx\ w(x,z)\\ =\int_0^{\infty}dx\frac{e^x}{(e^x+1)^2}\log{\left(\frac{e^{z^2/4x+x}+1}{e^{z^2/4x+x}-e^x}\right)}\\ =\int_0^{\infty}dx\frac{e^x}{(e^x+1)^2}\log{\left(\frac{e^{z^2/4x}+e^{-x}}{e^{z^2/4x}-1}\right)}$

其中。该函数的数值积分是较大代码的一部分，其中的结果用作其他函数的输入。到目前为止，我已将被积函数写为 $z > 0$ $W(z)$

double integrand__W(double x, double z){
    double arg = z*z/(4.0*x);
    
    double num = exp(arg+x)+1;
    double den1 = expm1(arg);
    double den2 = exp(x);

    num = isinf(num) ? arg+x : log(num);
    den1 = isinf(den1) ? arg : log(den1);
    den2 = x; //log(exp(x))=x
    double t1 = num-den1-den2;
    
    num = exp(x);
    double den = exp(x)+1;
    double t2 = isinf(den) ? exp(-x) : num/(den*den);
    
    return t1*t2;
}

其情节如下图所示：

对于数值积分，我使用的是Cubature，这是一个用于自适应多维积分的简单 C 包：

//integrator
struct fparams {
    double z;
};

int inf_W(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval){
    struct fparams * fp = (struct fparams *)fdata;
    double z = fp->z;
    double t = x[0]; 
    double aux = integrand__W(a_int+t*pow(1.0-t, -1.0), z)*pow(1.0-t, -2.0);
    if (!isnan(aux) && !isinf(aux))
    {
        fval[0] = aux;
    }
    else
    {
        fval[0] = 0.0;
    }
    return 0;
}

//range integration 1D
    size_t maxEval = 1e7;
    double xl[1] = { 0 };
    double xu[1] = { 1 };

    double W, W_ERR;
    struct fparams params = {z};
    hcubature(1, inf_W, &params, 1, xl, xu, maxEval, 0, 1e-5, ERROR_INDIVIDUAL, &W, &W_ERR);
    cout << "z: " << z << " | " << W << " , " << W_ERR << endl;

其中通过变量的变化可以在半无限区间上进行积分：

\int_{a}^{\infty} f (x) d x = \int_{0}^{1} f (a + \frac{t}{1 - t}) \frac{1}{(1 - t)^{2}} d t

$\int_a^{\infty}f(x)dx=\int_0^1 f\left(a+\frac{t}{1-t}\right)\frac{1}{(1-t)^2}dt$

分析上，我知道是非负的，所以积分本身应该是非负的。但是，由于缺乏准确性，我得到了一些不正确的结果： $w(x,z)$

z: 100 | -3.97632e-17 , 1.24182e-16

在Mathematica中，以高精度工作，我可以获得所需的结果：

w[x_, z_] := E^x/(E^x + 1)^2 Log[(E^(z^2/(4 x)) + E^-x)/(E^(z^2/(4 x)) - 1)]

W[z_?NumericQ] := NIntegrate[w[x, z], {x, 0, ∞},
  WorkingPrecision -> 40,
  Method -> "LocalAdaptive"]

W[100]

(* 4.679853458969239635780655689865016458810*10^-43 *)

到目前为止，我已经听到了关于这个问题的两种可能的诊断：

对于较大的，被积函数非常小，积分需要高精度。 $z$
积分方案允许负权重，有时可能会导致积分的准确性更高，但这里会导致假阴性。

我的问题：关于如何解决这个问题的任何想法，要么通过重写我的被积函数以达到所需的精度，要么切换到不同的集成方案？我需要在我的代码中解决这个积分，所以使用Mathematica不是一种选择。

3个回答

这可能是计算第二项时的准确性问题，因为当 $x \gg 1$ .

我会首先研究这个词：收集 $e^x$ 从分子和分母中取出并将其简化（留下一个和 $e^{-x}$ 在分子中）。

然而，与其猜测可能有必要更好地理解准确性问题的所在：你能画出被积函数吗？ $f(x)$ 对于不同的值 $x$ ? 并将其与 Mathematica 以更高精度计算的“精确”值进行比较？

如果您可以使用 Boost，您将获得exp-sinh quadrature的好处，它似乎可以很好地解决这个问题：

#include <iostream>
#include <cmath>
#include <boost/math/quadrature/exp_sinh.hpp>

using boost::math::quadrature::exp_sinh;
using std::exp;
using std::expm1;
using std::log;

int main() {
    exp_sinh<double> integrator;
    double z = 100.0;
    auto f = [z](double x) {
        double k1 = 1.0/(2 + exp(-x) +exp(x));
        double t = z*z/(4*x);
        double log_arg;
        if (t > 1) {
            log_arg = (1 + exp(-x)*exp(-t))/(1 - exp(-t));
        } else {
            log_arg = (exp(t) + exp(-x))/expm1(t);
        }
        return k1*log(log_arg);
    };
    double termination = sqrt(std::numeric_limits<double>::epsilon());
    double error;
    double L1;
    double Q = integrator.integrate(f, termination, &error, &L1);
    std::cout << "Q = " << Q << ", error estimate: " << error << "\n";
}

输出：

Q = 2.56597e-45, error estimate: 6.30636e-46

如果我double改为boost::multiprecision::float128，我得到

Q = 4.67986e-43, error estimate: 3.59871e-49

编辑：使其达到的差异boost::multiprecision::float128：

#include <iostream>
#include <cmath>
#include <boost/math/quadrature/exp_sinh.hpp>
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
using boost::math::quadrature::exp_sinh;
using std::exp;
using std::expm1;
using std::log;

int main() {
    exp_sinh<float128> integrator;
    float128 z = 100;
    auto f = [z](float128 x) {
        float128 k1 = 1/(2 + exp(-x) +exp(x));
        float128 t = z*z/(4*x);
        float128 log_arg;
        if (t > 1) {
            log_arg = (1 + exp(-x)*exp(-t))/(1 - exp(-t));
        } else {
            log_arg = (exp(t) + exp(-x))/expm1(t);
        }
        return k1*log(log_arg);
    };
    float128 termination = sqrt(std::numeric_limits<float128>::epsilon());
    float128 error;
    float128 L1;
    float128 Q = integrator.integrate(f, termination, &error, &L1);
    std::cout << "Q = " << Q << ", error estimate: " << error << "\n";
}
```

（我对 Cubature 没有任何经验，我不知道他们实现了哪些算法。）

如果您使用的积分规则具有负权重，例如 Newton-Cotes 规则，那么即使被积函数严格为正，您得到的结果也可能为负。这是我们在使用 MFEM 计算范数时经常会遇到的问题（可能其他 FEM 包也会有同样的问题）。但是，在这些情况下，您可以严格证明只有符号是错误的。所以对我们来说一个简单的解决方法是取结果的绝对值。

您的解决方案可能是使用具有正权重的正交规则，例如许多高斯正交规则之一。

其它你可能感兴趣的问题

上一篇如何在 MATLAB 中生成半正态变量？下一篇的正根Xq+ b x - bxq+bx−b