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
munflfor 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. – Michael E2 May 01 '18 at 12:07