積分または四重積分を使用した数値積分の計算

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ため、f内部exp()関数の引数の値があります巨大。実際に、私は指数（ない値）に引数を返すために機能を手直し場合：

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)のいくつかの要素が入り込み、被積分関数を計算可能な範囲にリスケールします。

いずれにしても、数字を台無しにするようなものなので、すべてのユニットをチェックすることをお勧めします。

NB：あなたは計算することができ、最大exp()は約exp(709.7827128933840)