Help with Hartree-Fock integral

I am trying to integrate an effective potential function for a particle in a box which takes into account the repulsion between two electrons. The integral is

Integrate[E^(x2 - 1)^2 1/(x1 - x2) E^(x2 - 1)^2, {x2, 0, 2}]


and I keep getting the integral returned again with no answer. Any help would be appreciated.

• What answer do you expect to get? Is x1 in (0,2)? Any particular reason you are you writing the exponential twice, rather than just E^(2 (x2-1)^2)? Dec 5, 2016 at 21:21
• No particular reason, it was just how I wrote it. x1 should be between 0 and 2 because it is the position of a particle on a distance between 0 and 2.
– Jim
Dec 5, 2016 at 21:24
• Are there any constraints on x1? Dec 5, 2016 at 23:01
• @jim what prevents this integral from blowing up when x2 == x1? Dec 5, 2016 at 23:03
• If you add the assumption: Integrate[E^(x2 - 1)^2 1/(x1 - x2) E^(x2 - 1)^2, {x2, 0, 2}, Assumptions -> 0 < x1 < 2] then you get the message that the integral does not converge on {0,2}. Dec 6, 2016 at 2:24

Would a numeric approach be sufficient?

hf[x1_?NumericQ] :=
E^(2 x1^2) (Log[(x1 + 1)/(1 - x1)]) +
NIntegrate[(E^(2 x2^2) - E^(2 x1^2))/(x1 - x2), {x2, -1, 1}];

Plot[hf[x1], {x1, -1, 1}]


Here I made a few transformations. I translated the problem from the interval {0, 2} to {-1, 1}, which made it easier for me to read (especially the symmetry). I also cleaved the singular part off and assumed that the Cauchy Principal Value was sought.

E^(2 x2^2)/(x1 - x2) == E^(2 x1^2)/(x1 - x2) + (E^(2 x2^2) - E^(2 x1^2))/(x1 - x2)


The integral of the first term is

Integrate[E^(2 x1^2)/(x1 - x2), {x2, -1, 1},
Assumptions -> -1 < x1 < 1, PrincipalValue -> True]
(*  E^(2 x1^2) Log[(1 + x1)/(1 - x1)]  *)


Alternatively, one can integrate a series expansion of the second term to get a rapidly convergent series approximation to the integral.

ClearAll[integral, integralterm, hf2];

(* integral of a term of the power series *)
integralterm[x1_ /; x1 == 0, k_] := x1;
integralterm[x1_, k_] =
Block[{x1, x2,
k}, -(2^k/k!) x1^(2 k - 1)
Sum[(x1 (x2/x1)^(j + 1))/(
j + 1), {j, 0, 2 k - 1}] /. {{x2 -> -1}, {x2 -> 1}} //
Differences // First // FullSimplify[#, -1 < x1 < 1] &
];

(* sum of the integrals of the terms up to order n *)
integral[x_, n_] := Total@Table[integralterm[x, k], {k, n}];

(* the total integral, approximated with n terms of the series *)
hf2[x1_ /; x1 == 0, n_] := x1;
hf2[x1_] := hf2[x1, 24];      (* order 24 gives ~MachinePrecision, less rounding error *)
hf2[x1_, n_Integer] :=
(hf2[x_, n] = Block[{x}, -E^(2 x^2) Log[(1 - x)/(1 + x)] + integral[x, n]];
hf2[x1, n]);


It's faster and arbitrary accuracy is possible by including more terms.

Plot[hf2[x1], {x1, -1, 1}] // AbsoluteTiming
Plot[hf[x1], {x1, -1, 1}] // AbsoluteTiming


Machine precision begins to be obtained around n == 20. The convergence at the ends of the interval is of course slower. At n == 24, convergence stops improving (using machine precision inputs) because of rounding error. I use a high-precision version of the NIntegrate method above for checking the accuracy of the series expansion.

hf[x_?NumericQ, wp_] :=(* high precision *)
With[{x1 = SetPrecision[x, wp + 0.5]},
E^(2 x1^2) (Log[(x1 + 1)/(1 - x1)]) +
NIntegrate[(E^(2 x2^2) - E^(2 x1^2))/(x1 - x2), {x2, -1, 1}, WorkingPrecision -> wp]]

ListPlot[Transpose@Table[
RealExponent@{(#1 - #2)/#2, #1 - #2} &[hf2[x1, 24], hf[x1, 32]],
{x1, Range[-1., 1., 1./1024][[2 ;; -2 ;; 2]]}],
PlotRange -> {-17., -13.}, Frame -> True,
GridLines -> {None, {RealExponent@\$MachineEpsilon}},
DataRange -> {-1, 1}, FrameLabel -> {x1, HoldForm[Log10[error]]}]


A plot of the relative error (blue) and absolute error (yellow).