使用積分或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)