3
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ 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)? $\endgroup$
    – evanb
    Dec 5, 2016 at 21:21
  • $\begingroup$ 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. $\endgroup$
    – Jim
    Dec 5, 2016 at 21:24
  • $\begingroup$ Are there any constraints on x1? $\endgroup$
    – bill s
    Dec 5, 2016 at 23:01
  • $\begingroup$ @jim what prevents this integral from blowing up when x2 == x1? $\endgroup$
    – evanb
    Dec 5, 2016 at 23:03
  • 2
    $\begingroup$ 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}. $\endgroup$
    – bill s
    Dec 6, 2016 at 2:24

1 Answer 1

6
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

$\endgroup$

Your Answer

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

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