4
$\begingroup$

I have a simulation code I developed in Mathematica 10.2. I use Nintegrate to calculate some values. It works fine and each run takes about 170s.

However When I run it in my university's computer (same specs ) that has the Mathematica 11.3. I get a bunch of errors in the Nintegrate part of the code.
(General::munfl: Exp[-21236.2] is too small to represent as a normalized machine number; precision may be lost.
General::stop: Further output of General::munfl will be suppressed during this calculation.
Infinite expression 1/(0. +0. I) encountered.)

I suspected some problems with numerical precision, So I manually changed the precision of all declared Machine numbers to 24 and the errors stopped. But the code now takes 370s per run...
I can't afford that difference in timing.

My question: Does anyone know how to run the Nintegrate function as it was before ?

Here is a simplified version of the code I run:

a = 25;   AreaFactor = 2.45 10^9;  cM = 0.262468;     
mL = 1.0;   mW = 0.8;  mR = 1.0; meff1 = 0.4;      meff2 = 0.4;     
e1 = 10.0;       e2 = 10.0;     
Lt1 = 10.0;     Lt2 = 10.0;    
U1EF = 1.2;      U2EF = 1.2;      
pFLMajor = 1.09;       pFLMinor = 0.421;      
pFWellMajor = 0.115 ; pFWellMinor = 0.115 ;
pFRMajor = 1.09;       pFRMinor = 0.421;        
cML = mL*cM; cMW = mW*cM; cMR = mR*cM;
cI1 = meff1*cM;  cI2 = meff2*cM;  
Tem = 2.5 ;(* Kelvin *) kB = 8.61733*10^(-5); (* Kelvin^-1 eV*)
Xx1 = -N[Sqrt[2.0*cMW*kB*Tem]]; Xx2 = +N[Sqrt[2.0*cMW*kB*Tem]];
nni = 110;   
ddx = (Xx2 - Xx1)/nni;

(************************* END OF INPUT ********************************)


up = 1;  
do = 2; 

EFermiL = (pFLMajor^2 + pFLMinor^2)/(2*cML);           
EFermiW = (pFWellMajor^2 + pFWellMinor^2)/(2*cMW);            
EFermiR = (pFRMajor^2 + pFRMinor^2)/(2*cMR);             
EFermi = (EFermiL + EFermiW + EFermiR)/3;   

U1POS = EFermi + U1EF;          
U2POS = EFermi + U2EF;         

eV1POS[eV_] := (e2*Lt1)/(e2*Lt1 + e1*Lt2)*
  eV; (* Voltage at the Left barrier *)

eV2POS[eV_] := (e1*Lt2)/(e2*Lt1 + e1*Lt2)*
  eV     (* Voltage at the Right barrier *)

(**************************************************************)

kLsPOS[1] := pFLMajor;
kLsPOS[2] := pFLMinor; 

kWsPOS[1, eV_] := Sqrt[pFWellMinor*pFWellMinor + cMW*eV1POS[eV]];
kWsPOS[2, eV_] := Sqrt[pFWellMajor*pFWellMajor + cMW*eV1POS[eV]];

kRsPOS[1, eV_] := 
  Sqrt[pFRMinor*pFRMinor + cMR*(eV1POS[eV] + eV2POS[eV])];
kRsPOS[2, eV_] := 
  Sqrt[pFRMajor*pFRMajor + cMR*(eV1POS[eV] + eV2POS[eV])];

(* POSITIVE VOLTAGE *) 

\[Delta]LWPOS[m_, eV_] := kLsPOS[m]/kWsPOS[m, eV];       


\[Delta]LRPOS[m_, eV_] := kLsPOS[m]/kRsPOS[m, eV];      

xLWminPOS[m_, eV_] := 
  Re[Sqrt[1 - 1/(\[Delta]LWPOS[m, eV]*\[Delta]LWPOS[m, eV])]];

xLRminPOS[m_, eV_] := 
  Re[Sqrt[1 - 1/(\[Delta]LRPOS[m, eV]*\[Delta]LRPOS[m, eV])]];

xLLimitPOS[m_, eV_] := 
  Max[{xLWminPOS[m, eV], 
    xLRminPOS[m, eV]}];    (* integration limit  Positive voltage*)


xWPOS[m_, eV_, xL_] := 
  Sqrt[1.0 - \[Delta]LWPOS[m, eV]*\[Delta]LWPOS[m, eV]*(1.0 - xL*xL)];
xRPOS[m_, eV_, xL_] := 
  Sqrt[1.0 - \[Delta]LRPOS[m, eV]*\[Delta]LRPOS[m, eV]*(1.0 - xL*xL)];


kLPOS[m_, eV_, xL_] := kLsPOS[m]*xL;
kWPOS[m_, eV_, xL_, x_] := 
  kWsPOS[m, eV]*xWPOS[m, eV, xL] + x*xWPOS[m, eV, xL];
kRPOS[m_, eV_, xL_] := kRsPOS[m, eV]*xRPOS[m, eV, xL];


(**************************************************************)

(*   Airy function arguments Positive voltages *) 

t1POS[eV_] := (cI1* eV1POS[eV]/Lt1)^(1/3);

t2POS[eV_] := (cI2* eV2POS[eV]/Lt2)^(1/3);
xLt1[s_, eV_] := (Lt1*(1/cML*kLsPOS[s]*kLsPOS[s] - U1POS))/ eV1POS[eV];

xLt2[s_, eV_, 
   W_] := (Lt2*(1/cMW*kWsPOS[s, eV]*kWsPOS[s, eV] - U2POS))/
    eV2POS[eV] - Lt1 - W;  

q1POS[s_, eV_] := t1POS[eV]*xLt1[s, eV];
q2POS[s_, eV_] := t1POS[eV]*(Lt1 + xLt1[s, eV]);
q3POS[s_, eV_, W_] := t2POS[eV]*(Lt1 + W + xLt2[s, eV, W]);
q4POS[s_, eV_, W_] := t2POS[eV]*(Lt1 + W + Lt2 + xLt2[s, eV, W]);

p1POS[s_, eV_] := 
  t1POS[eV]*(AiryAiPrime[-q2POS[s, eV]] AiryBi[-q2POS[s, eV]] - 
     AiryAi[-q2POS[s, eV]] AiryBiPrime[-q2POS[s, eV]]);
p2POS[s_, eV_, W_] := 
  t2POS[eV]*(AiryAiPrime[-q4POS[s, eV, W]] AiryBi[-q4POS[s, eV, W]] - 
     AiryAi[-q4POS[s, eV, W]] AiryBiPrime[-q4POS[s, eV, W]]);

a1POS[s_, eV_] := AiryAi[-q2POS[s, eV]] AiryBi[-q1POS[s, eV]] -
   AiryAi[-q1POS[s, eV]] AiryBi[-q2POS[s, eV]];

b1POS[s_, eV_] := 
  t1POS[eV]*(AiryAiPrime[-q2POS[s, eV]] AiryBi[-q1POS[s, eV]] - 
     AiryAi[-q1POS[s, eV]] AiryBiPrime[-q2POS[s, eV]]);

g1POS[s_, eV_] := 
  t1POS[eV]*(AiryAiPrime[-q1POS[s, eV]] AiryBi[-q2POS[s, eV]] - 
     AiryAi[-q2POS[s, eV]] AiryBiPrime[-q1POS[s, eV]]);
x1POS[s_, eV_] := 
  t1POS[eV]*
   t1POS[eV]*(AiryAiPrime[-q2POS[s, eV]] AiryBiPrime[-q1POS[s, eV]] - 
     AiryAiPrime[-q1POS[s, eV]] AiryBiPrime[-q2POS[s, eV]]);

a2POS[s_, eV_, W_] := 
  AiryAi[-q4POS[s, eV, W]] AiryBi[-q3POS[s, eV, W]] -
   AiryAi[-q3POS[s, eV, W]] AiryBi[-q4POS[s, eV, W]];
b2POS[s_, eV_, W_] := 
  t2POS[eV]*(AiryAiPrime[-q3POS[s, eV, W]] AiryBi[-q4POS[s, eV, W]] - 
     AiryAi[-q4POS[s, eV, W]] AiryBiPrime[-q3POS[s, eV, W]]);
g2POS[s_, eV_, W_] := 
  t2POS[eV]*(AiryAiPrime[-q4POS[s, eV, W]] AiryBi[-q3POS[s, eV, W]] - 
     AiryAi[-q3POS[s, eV, W]] AiryBiPrime[-q4POS[s, eV, W]]) ;
x2POS[s_, eV_, W_] := 
  t2POS[eV]*
   t2POS[eV]*(AiryAiPrime[-q4POS[s, eV, W]] AiryBiPrime[-q3POS[s, eV, 
         W]] - AiryAiPrime[-q3POS[s, eV, W]] AiryBiPrime[-q4POS[s, eV,
          W]]);

(*****************************************************************************)
\

(* Amplitude of the transmission *)

AmpVerPOS[s_, eV_, xL_, W_, 
   x_] := (4 E^(
     I (-(Lt1 + Lt2 + W) kRPOS[s, eV, xL] + W kWPOS[s, eV, xL, x]))
      meff1*meff2*mR*mW*p1POS[s, eV]*p2POS[s, eV, W]*kLPOS[s, eV, xL]*
     kWPOS[s, eV, xL, 
      x])/((meff1*
         kLPOS[s, eV, 
          xL]*(I meff1*a1POS[s, eV]*kWPOS[s, eV, xL, x] + 
           mW*b1POS[s, eV]) + 
        mL*(meff1*g1POS[s, eV]*kWPOS[s, eV, xL, x] + 
           I*mW*x1POS[s, eV])) (meff2*
         kRPOS[s, eV, 
          xL]*(I*meff2*a2POS[s, eV, W]*kWPOS[s, eV, xL, x] + 
           mW*b2POS[s, eV, W]) + 
        mR*(meff2*g2POS[s, eV, W]* kWPOS[s, eV, xL, x] + 
           I*mW*x2POS[s, eV, W])) + 
     E^(2 I W kWPOS[s, eV, xL, 
        x]) (meff1*
         kLPOS[s, eV, 
          xL]*(meff1*a1POS[s, eV]* kWPOS[s, eV, xL, x] + 
           I*mW*b1POS[s, eV]) + 
        mL *(-I*meff1* g1POS[s, eV]*kWPOS[s, eV, xL, x] - 
           mW*x1POS[s, eV])) (meff2*
         kRPOS[s, eV, 
          xL]*(meff2*a2POS[s, eV, W]* kWPOS[s, eV, xL, x] + 
           I*mW *b2POS[s, eV, W]) + 
        mR *(-I meff2*g2POS[s, eV, W]*kWPOS[s, eV, xL, x] - 
           mW*x2POS[s, eV, W])));


(* transmission coefficient *)

DTranPOS[s_, eV_, xL_, W_, x_] := (mL*kRPOS[s, eV, xL])/(
   mR*kLPOS[s, eV, xL])*AmpVerPOS[s, eV, xL, W, x]*
   Conjugate [AmpVerPOS[s, eV, xL, W, x ]]; 


(* Voltage grid*)

NStep = 91;          
VminP = 0.001;           (* "zero" value of the applied bias-voltage *)  
VmaxP = 0.550;       


VstepP = (VmaxP - VminP)/NStep;



AverxLDTranPOS[s_, eV_, W_] := (Xx2 - Xx1)/nni* Sum[

    NIntegrate[
     xL*DTranPOS[s, eV, xL, W, x], {xL, xLLimitPOS[s, eV], 1.0}, 
     Method -> "LocalAdaptive"], {x, Xx1, Xx2, ddx}];


(********************************************************) 

W = 9;        (* Well/Nanoparticle Diameter *) 

dp1 = W/10; (*nm*)

c = (W*W)/(8.0*25812.8) ;           (*  Angstrom^2/Ohm  *)

GzeroPOS[s_] := c*kLsPOS[s]*kLsPOS[s];
time1 = AbsoluteTime[];
(* Spin channels transport for positive voltages*)

listAverxLDTranUpDoDoP = 
  ParallelTable[{eV, Re[GzeroPOS[up]*AverxLDTranPOS[1, eV, W]]}, {eV, 
    VminP, VmaxP, VstepP}];

listAverxLDTranDoUpUpP = 
  ParallelTable[{eV, Re[GzeroPOS[do]*AverxLDTranPOS[2, eV, W]]}, {eV, 
    VminP, VmaxP, VstepP}];

time2 = AbsoluteTime[];
timetaken = time2 - time1
$\endgroup$
  • 5
    $\begingroup$ It's the handling of underflow that has changed. See this comment and search the site for munfl for similar issues related to this change. Maybe write to WRI. I don't know of a simple fix that resets the behavior to pre V11.3. $\endgroup$ – Michael E2 May 1 '18 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.