# Numerical computation of a double integral

I wish to numerically evaluate a function of the form

$$P(v_0) = \frac{e^{\int_0^{v_0}F(v)dv}}{\int_{-\infty}^{+\infty}e^{\int_0^{v_1}F(v)dv}dv_1}$$

In the denominator is the normalization constant which gives me the trouble.

The function $$F$$ is given by

$$F(v) = -v+\frac{1}{1+\frac{v}{e^v-1}}-\frac{1}{1+\frac{ve^v}{e^v-1}}$$

I can compute the nominator numerically:

P[v0_] :=
Exp[NIntegrate[-v + 1/(1 + v/(-1 + E^v)) - 1/(
1 + (E^v v)/(-1 + E^v)), {v, 0, v0}]];
Plot[P[v0], {v0, -20, 20}]


but can't manage to do it with the normalization constant included.

Any suggestions?

There's no need for finesse or truncation here:

F[v_] = -v + 1/(1 + v/(E^v - 1)) - 1/(1 + (v E^v)/(E^v - 1));
G[v1_?NumericQ] := Exp[NIntegrate[F[v], {v, 0, v1}]]
Gnorm = NIntegrate[G[v1], {v1, -\[Infinity], \[Infinity]}]
(*    2.87582    *)
P[v0_?NumericQ] := G[v0]/Gnorm

• Can you explain the importance of your use of Set here, as opposed to SetDelayed? – CA Trevillian Mar 29 '20 at 18:49
• @CATrevillian Set (=) is for definitions that I want done only once, like defining F[v] or evaluating Gnorm, whereas SetDelayed (:=) is for definitions that I want done afresh at every call. I'd say that for defining F[v] we can do either = or := ad libitum, but for Gnorm we absolutely must use an immediate = so that its numerical integral is evaluated only once, not afresh at every call. – Roman Mar 29 '20 at 19:08

The function F decays so rapidly that it should be okay to truncate the domain of integration. For example, we have

NIntegrate[Q[v1], {v1, -30, 30}, PrecisionGoal -> 14]/ NIntegrate[Q[v1], {v1, -100, 100}, PrecisionGoal -> 14] - 1


2.66454*10^-15

So the following should have quite a few digits of precision:

F[v_] := -v + 1/(1 + v/(-1 + E^v)) - 1/(1 + (E^v v)/(-1 + E^v))
Q[v1_?NumericQ] := Exp[NIntegrate[F[v], {v, 0, v1}]];
const = NIntegrate[Q[v1], {v1, -30, 30}, PrecisionGoal -> 14]
P[v1_] := Q[v1]/const;


isn't the normalization constant the integral starting at 0 (otherwise, you need to find the function at -infinity where it blows up). So use the final value of F as the normalization constant and divide by it at the end.