How to calculate a hard definite integral?

Question

I have the following function of $4$ - variables

f1111[x1_, x11_, x2_, x22_] := (-(x1 - x11)^4 - 6(x1 - x11)^2 (x2 - x22)^2
                                 + 3(x2 - x22)^4)/(4π((x1 - x11)^2 + (x2 - x22)^2)^3)

the variables are x1, x11, x2 and x22.

Wolfram Mathematica can easily calculate the indefinite integral of this function.

FullSimplify[ Integrate[ f1111[x1, x11, x2, x22], x1, x11, x2, x22]]

The result is a very long function. However, my problem is a definite integral, which needs to be calculated in a square domain,

Integrate[ f1111[x1, x11, x2, x22], 
             {x1, -0.5, -0.25}, {x11, -0.25, 0}, {x2, -0.5, -0.25}, {x22, -0.5, -0.25}]

The above code doesn't work, waiting for about 30 minutes nothing happened.

I tried numerical integration with different methods (LocalAdaptive and GlobalAdaptive), and I also tried to increase WorkingPrecision and PrecisionGoal, but the problem of singularity and low convergence is always present.

Is there a way to get the correct solution in a reasonable time?

Because I will need to perform this integration several times.

Answer 1

We couldn't be well-satisfied taking the black box system results when having in mind that no computer system is free of bugs and that its various aspects related to symbolic integration may seem not quite perfect (see e.g. this remarks). And so our purpose is to calculate this integral with more laborious approach, at least with a few steps which could be checked with a pencil at hand. Nevertheless the system will tackle the problem clearly in a few times more efficient way.

The integrand appears to be quite symmetric and it is convenient to change variables and define a new function. We can see that

f1111[x + w, -x + w, y + z, -y + z] // Simplify

 -((x^4 + 6 x^2 y^2 - 3 y^4)/(16 π (x^2 + y^2)^3))

and so we define new variables $$(x_1,x_{11}, x_2, x_{22}) \mapsto (x,w,y,z)=(\frac{x_1+x_{11}}{2}, \frac{x_1-x_{11}}{2}, \frac{x_2+x_{22}}{2},\frac{x_2-x_{22}}{2})$$ and

f[x_, y_] := -((x^4 + 6 x^2 y^2 - 3 y^4)/(16 π (x^2 + y^2)^3))

we calculate the jacobian of this linear transformation

Det @ D[{ x + w, w - x, y + z, z - y}, {{x, w, y, z}}]

 4

The function f doesn't depend on w and z and so integration simplifies, however one has to calculate the integration region appropriately. We will integrate the function in two steps, first integrating f with respect x and w (the result will be stored in if) and then if with respect to y and z. Now we find appropriate ranges of variables.

x1x11 = {{-1/2, -1/4}, {-1/4, -1/4}, {-1/4, 0}, {-1/2, 0}, {-1/2, -1/4}};
x2x22 = {{-1/2, -1/2}, {-1/4, -1/2}, {-(1/4), -(1/4)}, {-1/2, -(1/4)}, {-1/2, -1/2}};
tsf = {(#1 - #2)/2, (#1 + #2)/2} &;

now we have

xw = tsf @@@ x1x11
yz = tsf @@@ x2x22

{{-(1/8), -(3/8)}, {0, -(1/4)}, {-(1/8), -(1/8)}, {-(1/4), -(1/4)}, {-(1/8), -(3/8)}}
{{0, -(1/2)}, {1/8, -(3/8)}, {0, -(1/4)}, {-(1/8), -(3/8)}, {0, -(1/2)}}

GraphicsRow[{Graphics[{Magenta, Thickness[0.025], Line[x1x11]}], 
             Graphics[{Blue, Thickness[0.025], Line[xw]}]}]

This picture for x1x11 and xw is helpful for appropriate specifications of integration regions. Quite analogical is the region for x2x22 and yz. Integration has to be divided in two subsets for the both regions xw and yz, and so (a simple exercise is left for the reader why integration over w has been changed into multiplication by (2 x + 1/2) for -(1/4)< x <-(1/8)) and (-2 x) for -(1/8) < x < 0, as well as by (2 y + 1/4) for -1/8 < y < 0 and -2 y + 1/4 for 0 < y < 1/8)

and finally we evaluate:

AbsoluteTiming[
  if[y_] = Assuming[{-(1/8) < y < 1/8}, 
             4 (Integrate[f[x, y] (2 x + 1/2), {x, -(1/4), -(1/8)}] 
                 + Integrate[f[x, y] (-2 x), {x, -(1/8), 0}])];
  int = Simplify[ Integrate[if[y] (2 y + 1/4), {y, -1/8, 0}] 
                  + Integrate[if[y] (-2 y + 1/4), {y, 0, 1/8}]]]

{11.3938, (-4π  + 8 ArcCot[2] + 8 ArcTan[1/2] + 28 Log[2] - 11 Log[5])/(128 π )}

This timing is for my old and slow laptop (i3, 1.9GHz, 4 GB RAM )

N[int, 20]

 -0.0085637914008856298463

Another way is exploiting appropriately ImplicitRegion, however we have performed this way because the transformation tsf is frequently used in physics e.g. retarded and advanced variables.

Answer 2

You can also do step by step indefinite integration.

With Rubi https://rulebasedintegration.org/ i easily got the right result. Unfortunately indefinite integration with Integrate failed. Didn't examine further why.

assum = Assumptions -> {-1/2 < x1 < -1/4, -1/4 < x11 < 0, -1/2 < 
 x2 < -1/4, -1/2 < x22 < -1/4};
int1[x1_] = 
Int[(-(x1 - x11)^4 - 6 (x1 - x11)^2 (x2 - x22)^2 + 3 (x2 - x22)^4)/(
4 [Pi] ((x1 - x11)^2 + (x2 - x22)^2)^3), x1];
lim1 = Limit[int1[x1], x1 -> -1/4, Direction -> 1, assum] - 
Limit[int1[x1], x1 -> -1/2, Direction -> -1, assum];
int2[x2_] = Int[lim1, x2];
lim2 = Limit[int2[x2], x2 -> -1/4, Direction -> 1, assum] - 
Limit[int2[x2], x2 -> -1/2, Direction -> -1, assum] // 
Simplify[#, assum] &;
int3[x11_] = Int[lim2, x11];
lim3 = Limit[int3[x11], x11 -> 0, Direction -> 1, assum] - 
Limit[int3[x11], x11 -> -1/4, Direction -> -1, assum] // 
Simplify[#, assum] &;
int4[x22_] = Int[lim3, x22];
lim4 = Limit[int4[x22], x22 -> -1/4, Direction -> 1] - 
Limit[int4[x22], x22 -> -1/2, Direction -> -1] //   Simplify[#] &;
lim4 // ComplexExpand // N
(*   -0.008563791400885638`   *)