Maximising an Integral

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.

• 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"?) Mar 20 '18 at 18:14

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}]


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}}

• 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? Mar 20 '18 at 18:32
• 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. Mar 20 '18 at 18:45
• Yes, that's the sort of restriction I was thinking of. Mar 20 '18 at 20:47
• I have checked that something like v > 1000000000001/1000000000000 is sufficient. Mar 20 '18 at 20:54
• 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. Mar 20 '18 at 22:10