2
$\begingroup$

I do not understand why I cannot maximise the following integral by changing v

rrrintegral[egap_?NumericQ, v_?NumericQ] := 
  Module[{}, 
   NIntegrate[energy^2/(Exp[(energy - v)] - 1), {energy, egap, 10}]];

Maximize[{v*rrrintegral[1, v], v > 0}, v]

I get the following error:

The integrand has evaluated to non-numerical values for all sampling points in the region with boundaries

I have tried using "?NumericQ" but with no luck.

$\endgroup$
1
  • 1
    $\begingroup$ I get no error, but Maximize does not evaluate. Maximize is a symbolic solver, which requires a symbolic expression it can analyze. I think you need a numeric solver, such as FindMaximum or NMaximize. (But you also need to do something about the singularity in the integral as energy == v. Maybe Method -> "PrincipalValue"?) $\endgroup$
    – Michael E2
    Mar 20 '18 at 18:14
6
$\begingroup$

Perhaps use FindMaximum or NMaximize?

There's a singularity in the integrand when egap < v < 10 that needs to be addressed. My guess is Method -> "PrincipalValue" is what is intended. For that method, one needs to specify the singular points in the iterator: {energy, egap, v, 10}.

rrrintegral[egap_?NumericQ, v_?NumericQ] := If[egap < v < 10,
   NIntegrate[energy^2/(Exp[(energy - v)] - 1), {energy, egap, v, 10},
     Method -> "PrincipalValue"],
   NIntegrate[energy^2/(Exp[(energy - v)] - 1), {energy, egap, 10}]];

NMaximize[{v*rrrintegral[1, v], v > 0}, v]
(*  {33.4991, {v -> 3.12388}}  *)

Check the big picture (or do this first and use FindMaximum, which head-to-head is usually faster than NMaximize):

Plot[{v*rrrintegral[1, v]}, {v, 0, 5}]

Mathematica graphics

$\endgroup$
6
$\begingroup$

Quite a straightforward approch would be evaluating the integral symbolically, and then maximizing the appropriate function. Using Assumptions we can restrict the integral to an appropriate domain, e.g. (it takes a bit to evaluate):

Integrate[ energy^2/(Exp[(energy - v)] - 1), {energy, egap, 10}, 
           Assumptions -> 10 > egap > v > 0]
1/3 (-1000 + egap^3 - 3 I egap^2 Pi + 300 Log[1 - E^(10 - v)] - 
  3 egap^2 Log[-1 + E^(egap - v)] + 60 PolyLog[2, E^(10 - v)] - 
  6 egap PolyLog[2, E^(egap - v)] - 6 PolyLog[3, E^(10 - v)] + 
  6 PolyLog[3, E^(egap - v)])

Next we define

Int[egap_, v_] := 1/3 (-1000 + egap^3 + 300 Log[1 - E^(10 - v)] 
    -3 egap^2 Log[1 - E^(egap - v)] + 60 PolyLog[2, E^(10 - v)] 
    -6 egap PolyLog[2, E^(egap - v)] - 6 PolyLog[3, E^(10 - v)] 
    +6 PolyLog[3, E^(egap - v)])

Existing complex expression in the integral doesn't hurt since we can get rid off by changing the sign under Log. I have to use Re because of very small imaginary perturbations being a numerical artefact.

NMaximize[{v Int[1, v] // Re, 10 > v > 0}, v]
{33.4991, {v -> 3.12388}}
$\endgroup$
5
  • $\begingroup$ Actually, there may be a problem with both our answers: Doesn't the integral go to Infinity as v goes to 1 (if egap = 1)? So, no max., unless there are further restrictions on v? $\endgroup$
    – Michael E2
    Mar 20 '18 at 18:32
  • $\begingroup$ It seems you are right, however If I replace my approach by e.g. NMaximize[{v Int[1, v] // Re, 10 > v > 1.1}, v] it'll be correct. $\endgroup$
    – Artes
    Mar 20 '18 at 18:45
  • $\begingroup$ Yes, that's the sort of restriction I was thinking of. $\endgroup$
    – Michael E2
    Mar 20 '18 at 20:47
  • $\begingroup$ I have checked that something like v > 1000000000001/1000000000000 is sufficient. $\endgroup$
    – Artes
    Mar 20 '18 at 20:54
  • $\begingroup$ The unsatisfactory part of this closed form approach is the presence of terms of the form PolyLog[k, Exp[u]]; this is just asking for numerical trouble. Unfortunately, Bose-Einstein or Fermi-Dirac functions aren't yet built-in, so this will have to do. $\endgroup$
    – J. M.'s torpor
    Mar 20 '18 at 22:10

Your Answer

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

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