高斯函数的数值积分 - 不定积分

Question

我的方法

fun = @(y) (1/sqrt(pi))*exp(-(y-1).^2).*log(1 + exp(-4*y))
integral(fun,-Inf,Inf)

这给出了 NaN。

所以我尝试绘制它。

y= -10:0.1:10;
plot(y,exp(-(y-1).^2).*log(1 + exp(-4*y)))

然后明白域（重要部分）是从-4到+4。

因此将限制更改为

integral(fun,-10,10)

但是我不想总是绘制图表然后知道它的限制。那么有没有办法直接从-Inf到Inf知道积分。

Answer 1

讨论

如果你的积分总是这样的形式

我会使用高阶Gauss–Hermite quadrature rule。 它类似于构成 quadgk 基础的 Gauss-Legendre-Kronrod 规则，但它是专门为使用标准高斯乘数​​在实线上进行积分而定制的。

用替换x = y-1重写你的方程，我们得到

.

然后可以使用任意阶的 Gauss-Hermite 规则（在合理范围内）计算积分：

>> order           = 10;
>> [nodes,weights] = GaussHermiteRule(order);
>> f               = @(x) log(1 + exp(-4*(x+1)))/sqrt(pi);
>> sum(f(nodes).*weights)
ans =
    0.1933

我注意到下面的函数构建了一个完整的order x order 矩阵来计算nodes，所以它不应该太大。 有一种方法可以通过显式计算权重来避免这种情况，但我决定偷懒。 此外，100 阶事件，高斯乘数约为2E-98，因此被积函数的贡献极小。 虽然这本身不是自适应的，但在大多数情况下，高阶规则就足够了……我希望。

代码

function [nodes,weights] = GaussHermiteRule(n)
    % ------------------------------------------------------------------------------
    %  Find the nodes and weights for a Gauss-Hermite Quadrature integration.
    %

    if (n < 1)
        error('There is no Gauss-Hermite rule of order 0.');
    elseif  (n < 0) || (abs(n - round(n)) > eps())
        error('Given order ''n'' must be a strictly positive integer.');
    else
        n = round(n);
    end

    %   Get the nodes and weights from the Golub-Welsch function
    n = (0:n)'  ;
    b = n*0     ;
    a = b + 0.5 ;
    c = n       ;
    [nodes,weights] = GolubWelsch(a,b,c,sqrt(pi));

end


function [xk,wk] = GolubWelsch(ak,bk,ck,mu0)
    %GolubWelsch
    %   Calculate the approximate* nodes and weights (normalized to 1) of an orthogonal 
    %   polynomial family defined by a three-term reccurence relation of the form
    %       x pk(x) = ak pkp1(x) + bk pk(x) + ck pkm1(x)
    %   
    %   The weight scale factor mu0 is the integral of the weight function over the
    %   orthogonal domain.
    %

    %   Calculate the terms for the orthonormal version of the polynomials
    alpha = sqrt(ak(1:end-1) .* ck(2:end));

    %   Build the symmetric tridiagonal matrix
    T = full(spdiags([[alpha;0],bk,[0;alpha]],[-1,0,+1],length(alpha),length(alpha)));

    %   Calculate the eigenvectors and values of the matrix
    [V,xk] = eig(T,'vector');

    %   Calculate the weights from the eigenvectors - technically, Golub-Welsch requires
    %   a normalization, but since MATLAB returns unit eigenvectors, it is omitted.
    wk = mu0*(V(1,:).^2)';

end

Answer 2

我已经成功地使用数值变量变换转换了这种无限有界积分，如数值配方 3e，第 4.5.3 节中所述。基本上，您替换 y=c*tan(t)+b，然后在 (-pi/2,pi/2) 中对 t 进行数值积分，这会将 y 从 -infinity 扫描到 infinity。您可以调整 c 和 b 的值来优化过程。这种方法在很大程度上避免了尝试确定域中的截止点的问题，但是为了使用正交可靠地工作，您必须知道被积函数没有远离 y=b 的特征。

Answer 3

一个快速而肮脏的解决方案是寻找一个位置，你的函数足够小，然后将其作为限制。这假设对于x>0，函数fun 单调递减，对于所有x，fun(x) 的大小与fun(-x) 大致相同。

%// A small number
epsilon = eps;
%// Stepsize for searching bound
stepTest = 1;
%// Starting position for searching bound
position = 0;
%// Not yet small enough
smallEnough = false;

%// Search bound
while ~smallEnough
   smallEnough = (fun(position) < eps);
   position = position + stepTest;
end

%// Calculate integral
integral(fun, -position, position)

如果您对绘制函数感到满意，并通过眼睛决定可以剪切的位置，那么我想这段代码就足够了。