8
$\begingroup$

I have this code below, which is calculating the Binding energy of an electron in a Quantum Well Wire with a hydrogenic impurity in it. Well you don't have to care much about what kind of calculation it does, because it is returning the right number, my only problem is, that it's taking about 20 minutes for it to return a single value for the Eb function (you can try Eb[0.7, 1, 0.01]). I'm wondering, if there's a way to make this code run faster. As you can see I have written everything almost the same way as one would write on paper. I've searched and tried many different approaches to make it faster, but nothing has helped so far. 

e = 4.803*10^-10;
m = 0.067*9.109*10^-28;
h = 1.054*10^-27;
c = 2.997*10^10;
e0 = 13.18;

O1[Om_] = 10^13*Om;

oH[h0_] = (e*10000*h0)/(m*c);

Oc[h0_, Om_] = Sqrt[oH[h0]^2 + 4*O1[Om]^2];

aH[h0_, Om_] = Sqrt[h/(m*Oc[h0, Om])];

r0[rho_, phi_, z_, rhoi_] = 
  Sqrt[rho^2 + rhoi^2 - 2*rho*rhoi*Cos[phi] + z^2];

Psi[rho_, h0_, Om_] = E^(-(1/2)*(rho^2*aH[1, Om]^2)/aH[h0, Om]^2);

MGamma[rho_, phi_, z_, rhoi_, lambda_, Om_] = 
  E^(-lambda*r0[rho, phi, z, rhoi]);

CPhi[rho_, phi_, z_, rhoi_, lambda_, h0_, Om_] = 
  Psi[rho, h0, Om]*MGamma[rho, phi, z, rhoi, lambda, Om];

intCPhiCPhi[rhoi_, lambda_, h0_, Om_] := 
  NIntegrate[
   Abs[CPhi[rho, phi, z, rhoi, lambda, h0, Om]]^2*rho, {rho, 0, 
    Infinity}, {phi, 0, 2*\[Pi]}, {z, -Infinity, +Infinity}];

leftover[rho_, phi_, z_, rhoi_, lambda_, h0_, Om_] = 
  CPhi[rho, phi, z, rhoi, lambda, h0, Om]*
    h^2/(2*m*(aH[1, Om])^2)*(Psi[rho, h0, Om]/rho*
       D[MGamma[rho, phi, z, rhoi, lambda, Om], rho] + 
      2*D[MGamma[rho, phi, z, rhoi, lambda, Om], rho]*
       D[Psi[rho, h0, Om], rho] + 
      Psi[rho, h0, Om]*
       D[MGamma[rho, phi, z, rhoi, lambda, Om], {rho, 2}] + 
  Psi[rho, h0, Om]/rho^2*
   D[MGamma[rho, phi, z, rhoi, lambda, Om], +{phi, 2}] + 
  Psi[rho, h0, Om]*
   D[MGamma[rho, phi, z, rhoi, lambda, Om], {z, 2}]) + (e^2*
  Abs[CPhi[rho, phi, z, rhoi, lambda, h0, Om]]^2)/(e0*
  r0[rho, phi, z, rhoi]*aH[1, Om]);

IntLeftover[rhoi_?NumericQ, lambda_?NumericQ, h0_?NumericQ, 
  Om_?NumericQ] := (2*e0^2*h^2)/(e^4*
     m)*((1/intCPhiCPhi[rhoi, lambda, h0, Om] )*
    NIntegrate[
     leftover[rho, phi, z, rhoi, lambda, h0, Om]*rho, {rho, 0, 
      Infinity}, {phi, 0, 2*\[Pi]}, {z, -Infinity, +Infinity}])

Eb[rhoi_, h0_, Om_] := 
  FindMaximum[IntLeftover[rhoi, lambda, h0, Om], lambda];
Improve this question
$\endgroup$
7
  • $\begingroup$ the only obvious changes I'd make are 1. Change the 'calls' to D in leftover to eg DMGamma[rho_, phi_, z_, rhoi_, lambda_, 1] = D[MGamma[rho, phi, z, rhoi, lambda], rho] etc 2. make sure that MGamma is/is not dependent on Om (this shouldn't affect performance) 3. perhaps rewrite ah Oc with numeric quantities directly integrated into the definitions instead of relying on calling O1,oH (this could possibly affect performance) $\endgroup$
    – user42582
    Jan 7, 2018 at 11:36
  • $\begingroup$ @user42582 I have actually tried all of that suggestions before. That changes do not show any noticeable change in speed. $\endgroup$
    – tyler_house
    Jan 7, 2018 at 11:46
  • $\begingroup$ Does "SymbolicProcessing"->0 help? $\endgroup$
    – xzczd
    Jan 7, 2018 at 13:32
  • $\begingroup$ You would have to improve the time of the integrals. 3D integrals are often difficult. Your IntLeftoever takes 1 - 15 sec. to evaluate depending on lambda. -- What's the value for lamda for your trial example? I got {lambda -> 2.28525} in 17 sec. $\endgroup$
    – Michael E2
    Jan 7, 2018 at 13:40
  • $\begingroup$ @xzczd Adding the option Method -> {Automatic, "SymbolicProcessing" -> 0} doesn't help. $\endgroup$
    – tyler_house
    Jan 7, 2018 at 13:43

1 Answer 1

Reset to default
8
$\begingroup$

One way to speed up multidimensional integrals is to reduce the PrecisionGoal. That will be acceptable in some cases but not in others.

Modified OP's code:

Clear[e, m, h, c, e0]  (* I moved the parameter initialization to later (unimportant *)

O1[Om_] = 10^13*Om;
oH[h0_] = (e*10000*h0)/(m*c);
Oc[h0_, Om_] = Sqrt[oH[h0]^2 + 4*O1[Om]^2];
aH[h0_, Om_] = Sqrt[h/(m*Oc[h0, Om])];
r0[rho_, phi_, z_, rhoi_] = Sqrt[rho^2 + rhoi^2 - 2*rho*rhoi*Cos[phi] + z^2];
Psi[rho_, h0_, Om_] = E^(-(1/2)*(rho^2*aH[1, Om]^2)/aH[h0, Om]^2);
MGamma[rho_, phi_, z_, rhoi_, lambda_, Om_] = E^(-lambda*r0[rho, phi, z, rhoi]);
CPhi[rho_, phi_, z_, rhoi_, lambda_, h0_, Om_] = 
  Psi[rho, h0, Om]*MGamma[rho, phi, z, rhoi, lambda, Om];

intCPhiCPhi[rhoi_, lambda_, h0_, Om_] := 
  NIntegrate[Abs[CPhi[rho, phi, z, rhoi, lambda, h0, Om]]^2*rho, (* NB: Abs[..]^2=(..)^2 *)
   {phi, 0, 2*π}, {rho, 0, Infinity}, {z, -Infinity, +Infinity},
   PrecisionGoal -> 3];

leftover[rho_, phi_, z_, rhoi_, lambda_, h0_, Om_] = 
  CPhi[rho, phi, z, rhoi, lambda, h0, Om]*
    h^2/(2*m*(aH[1, Om])^2)*(Psi[rho, h0, Om]/rho*
       D[MGamma[rho, phi, z, rhoi, lambda, Om], rho] + 
      2*D[MGamma[rho, phi, z, rhoi, lambda, Om], rho]*D[Psi[rho, h0, Om], rho] + 
      Psi[rho, h0, Om]*D[MGamma[rho, phi, z, rhoi, lambda, Om], {rho, 2}] + 
      Psi[rho, h0, Om]/rho^2*D[MGamma[rho, phi, z, rhoi, lambda, Om], +{phi, 2}] + 
      Psi[rho, h0, Om]*D[MGamma[rho, phi, z, rhoi, lambda, Om], {z, 2}]) +
      (e^2*Abs[CPhi[rho, phi, z, rhoi, lambda, h0, Om]]^2) / 
       (e0*r0[rho, phi, z, rhoi]*aH[1, Om]);

IntLeftover[rhoi_?NumericQ, lambda_?NumericQ, h0_?NumericQ, 
  Om_?NumericQ] := (2*e0^2*h^2)/(e^4*
     m)*((1/intCPhiCPhi[rhoi, lambda, h0, Om])*
    NIntegrate[leftover[rho, phi, z, rhoi, lambda, h0, Om]*rho,
     {phi, 0, 2*π}, {rho, 0, Infinity}, {z, -Infinity, +Infinity},
      PrecisionGoal -> 3]);

Eb[rhoi_, h0_, Om_] := FindMaximum[IntLeftover[rhoi, lambda, h0, Om], lambda];

OP's example:

PrintTemporary@Dynamic@Clock@{0, Infinity};  (* running timer (helps preserve my sanity *)
Block[{  (* parameter initialization *)
  e = 4.803*10^-10,
  m = 0.067*9.109*10^-28,
  h = 1.054*10^-27,
  c = 2.997*10^10,
  e0 = 13.18},
 Eb[0.7, 1, 0.01] // AbsoluteTiming
 ]

FindMaximum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient increase in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

(*  {11.0647, {1.16175, {lambda -> 2.28525}}}  *)

Playing with it, I sometimes did not get a FindMaximum::lstol warning, but I always got 2.28525 or 2.28526 for an answer. I suspect that the noise from the numerical integration is the source.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Wow, thank you very much! I just checked the values of lambda, actually not setting PrecisionGoal->3 returns lambda -> 2.28403, but there's no significant difference in the Energy value, so this is exactly what I was looking for. Quick question, does the option Method -> {"GlobalAdaptive"} have any significant effect on the speed? $\endgroup$
    – tyler_house
    Jan 7, 2018 at 14:13
  • $\begingroup$ @tyler_house No, Method -> {"GlobalAdaptive"} virtually has no effect; it's the default, anyway. I meant to remove it. I was just trying things out. $\endgroup$
    – Michael E2
    Jan 7, 2018 at 14:17
  • $\begingroup$ Just one more question, is there a way to know what settings of PrecisionGoal will lead to better performance? I remember that I had tried different settings for PrecisionGoal with no luck. $\endgroup$
    – tyler_house
    Jan 7, 2018 at 14:20
  • 3
    $\begingroup$ @tyler_house The lower the faster. For single integrals, the default is around 8; for multivariate ones, it's 6. (Often the problem with slow integral are singularities in the integrand. Your integrand looks too complicated for me to make sure that is not the case here, but it doesn't seem to be the case.) $\endgroup$
    – Michael E2
    Jan 7, 2018 at 14:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.