Calcul de l'intégrale numérique à l'aide de l'intégrale ou du quadgk

Question

J'utilise MATLAB pour calculer l'intégrale numérique d'une fonction complexe incluant l'exposant naturel.

Je reçois un avertissement:

valeur infinie ou non un nombre rencontré

si j'utilise la fonction integral, tandis qu'une autre erreur est levée:

sortie de la fonction doit avoir la même taille que l'entrée

si j'utilise la fonction quadgk.

Je pense que la raison pourrait être que l'intégrande est infinie lorsque la variable ep est proche de zéro.

Code illustré ci-dessous. J'espère que vous pouvez m'aider à comprendre.

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

J'ai réécrite vos équations afin qu'un humain puisse comprendre réellement:

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 

Maintenant, pour la première itération, ep = 0.01, la valeur de l'argument de la fonction exp() à l'intérieur f est énorme. En fait, si je retravaille la fonction de retourner l'argument à l'exposant (pas la valeur):

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 

et imprimer sa valeur à quelques exemples de noeuds comme ceci:

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 

Je reçois ceci:

% 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 

Ainsi, integral() doit faire face à des choses comme exp(10^7). Cela peut ne pas être un problème en soi si l'argument tombe assez rapidement, mais comme indiqué ci-dessus, ce n'est pas le cas.

Donc, fondamentalement, vous demandez l'intégrale d'une fonction dont la valeur est comprise entre exp(~10^7) et exp(~10^3). Inutile de dire que le d(ef) de l'ordre de 0.01 ne va pas compenser cela, et ce sera non fini en arithmétique en virgule flottante.

Je suppose que vous avez un problème de mise à l'échelle. A en juger par les noms de vos variables ainsi que les équations, je pense que cela a quelque chose à voir avec la thermodynamique; une forme retravaillée de la loi de Planck? Dans ce cas, je vérifierais si vous travaillez en nanomètres; quelques facteurs de 10^(-9) s'insinueront, en réajustant votre integrand à la gamme computableablement computable.

Dans tous les cas, il serait judicieux de vérifier toutes vos unités, car c'est quelque chose comme ça qui gâche les chiffres.

NB: le exp() maximum que vous pouvez calculer est d'environ exp(709.7827128933840)