使用积分或quadgk计算数值积分

Question

我使用MATLAB来计算包括自然指数在内的复数函数的数值积分。

我得到一个警告：

无限或没有非数字值遇到

，如果我使用的功能integral，而另一个则会引发错误：

输出的功能必须与输入尺寸相同

如果我使用功能quadgk。

我认为其原因可能是当变量ep接近零时，被积函数是无限的。

代码如下所示。希望你们能帮我弄明白。

close all 
clear 
clc 
%% 
N = 10^5; 
edot = 10^8; 
yita = N/edot; 
kB = 8.6173324*10^(-5); 
T = 300; 
gamainf = 0.115; 
dTol = 3; 
K0 = 180; 
K = K0/160.21766208; 
nu = 3*10^12; 
i = 1; 
data = []; 
%% lambda = ec/ef < 1 
for ef = 0.01:0.01:0.1 
    for lambda = 0.01:0.01:0.08 
     ec = lambda*ef; 
     f = @(ep) exp(-((32/3)*pi*gamainf^3*(0.5+0.5*sqrt(1+2*dTol*K*(ep-ec)/gamainf)-dTol*K*(ep-ec)/gamainf).^3/(K*(ep-ec)).^2-16*pi*gamainf^3*(0.5+0.5*sqrt(1+2*dTol*K*(ep-ec)/gamainf)-dTol*K*(ep-ec)/gamainf).^2/((1+dTol*K*(ep-ec)/(gamainf*(0.5+0.5*sqrt(1+2*dTol*K*(ep-ec)/gamainf)-dTol*K*(ep-ec)/gamainf)))*(K*(ep-ec)).^2))/(kB*T)); 
     q = integral(f,0,ef,'ArrayValued',true); 
     % q = quadgk(f,0,ef); 
     prob = 1-exp(-yita*nu*q); 
     data(i,1) = ef; 
     data(i,2) = lambda; 
     data(i,3) = q; 
     i = i+1; 
    end 
end 

Answer 1

我已经重写你的方程，使人类能够真正理解它：现在

function integration 

    N  = 1e5; 
    edot = 1e8; 
    yita = N/edot; 
    kB  = 8.6173324e-5; 
    T  = 300; 
    gamainf = 0.115; 
    dTol = 3; 
    K0  = 180; 
    K  = K0/160.21766208; 
    nu  = 3e12; 
    i  = 1; 
    data = []; 

    %% lambda = ec/ef < 1 
    for ef = 0.01:0.01:0.1 
     for lambda = 0.01:0.01:0.08 
      ec = lambda*ef; 

      q = integral(@f,0,ef,'ArrayValued',true); 
      % q = quadgk(f,0,ef); 

      prob = 1 - exp(-yita*nu*q); 
      data(i,:) = [ef lambda q]; 

      i = i+1; 
     end 
    end 


    function y = f(ep) 

     G = K*(ep - ec); 
     r = dTol*G/gamainf; 
     S = sqrt(1 + 2*r); 
     x = (1 + S)/2 - r; 
     Z = 16*pi*gamainf^3; 

     y = exp(-Z*x.^2.*(2*x/(3*G.^2) - 1/(G.^2*(1 + r/x)))) /... 
            (kB*T)); 

    end 

end 

，对于第一次迭代，ep = 0.01的exp()函数的自变量的内部f值巨大的。事实上，如果我返工功能的参数返回指数（不是值）：

function y = f(ep) 

    % ... all of the above 

    % NOTE: removed the exp() to return the argument 
    y = -Z*x.^2.*(2*x/(3*G.^2) - 1/(G.^2*(1 + r/x)))) /... 
          (kB*T); 

end 

，并打印出它的值，在某些例子节点像这样：

for ef = 0.01 : 0.01 : 0.1 
    for lambda = 0.01 : 0.01 : 0.08 

     ec = lambda*ef; 

     zzz(i,:) = [f(0) f(ef/4) f(ef)]; 

     i = i+1; 
    end 
end 

zzz 

我得到这样的：

% f(0)      f(ef/4)     f(ef) 
zzz = 
    7.878426438111721e+07  1.093627454284284e+05  3.091140080273912e+03 
    1.986962280947140e+07  1.201698288371587e+05  3.187767404903769e+03 
    8.908646053687230e+06  1.325435523124976e+05  3.288027743119838e+03 
    5.055141696747510e+06  1.467952125661714e+05  3.392088351112798e+03 
    ... 
    3.601790797707676e+04  2.897200140791236e+02  2.577170427480841e+01 
    2.869829209254144e+04  3.673888685004256e+02  2.404148067956737e+01 
    2.381082059148755e+04  4.671147785149462e+02  2.238181495716831e+01 

所以，integral()必须处理像exp(10^7)。如果论证很快就会下降，这本身可能不是问题，但如上所示，它不会。

所以基本上你要求的是一个函数的积分，其范围在exp(~10^7)和exp(~10^3)之间。不用说，0.01的d(ef)不会补偿这一点，并且在浮点运算中它将是非有限的。

我怀疑你有缩放问题。从变量的名称和方程来看，我认为这与热力学有关。普朗克法则的修订形式？在这种情况下，我会检查你是否在纳米工作; 10^(-9)的一些因素将会蔓延，将您的积分重新调整到可计算的范围。

在任何情况下，检查所有单位是明智的，因为这就像是搞乱了数字。

注：您可以计算最大exp()大约是exp(709.7827128933840)