2
$\begingroup$

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.

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ NIntegrate throws up a couple of warnings but gives approx -0.00855780339 as the result quickly. Did you try this ? $\endgroup$ Commented Feb 27, 2020 at 10:04
  • $\begingroup$ Yes I did, but with Trapezoidal method I get different result, - 0.00709118. $\endgroup$ Commented Feb 27, 2020 at 10:07
  • $\begingroup$ So it's not a problem if Mathematica gives warnings, the result can still be correct? $\endgroup$ Commented Feb 27, 2020 at 10:09
  • 8
    $\begingroup$ Timing[Integrate[f1111[x1,x11,x2,x22],{x1,-1/2,-1/4},{x11,-1/4,0},{x2,-1/2,-1/4},{x22,-1/2,-1/4}]] returns {67.340763, (-11*Log[5]+4*(Pi-4*ArcTan[2]+Log[128]))/ (128*Pi)} with no warnings or errors and which should be exactly correct. That number is approximately -0.0085637914008856298463472158822972 $\endgroup$ Commented Feb 27, 2020 at 10:37
  • 4
    $\begingroup$ Mathematica 101: Never mix inexact numbers with functions meant to be used to obtain analytical exact result. $\endgroup$ Commented Feb 28, 2020 at 22:24

2 Answers 2

Reset to default
7
$\begingroup$

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

enter image description here

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ The integrand is a rational function containing no symbolic parameters, and the bounds are likewise numeric. If Integrate is going to have issues it will be with assessing convergence of a multiple integral, not with the sorts of things that arise in more general cases. (This comment only pertains to the opening sentence of the response; the rest appears to be quite useful). $\endgroup$ Commented Feb 28, 2020 at 16:44
  • 1
    $\begingroup$ @DanielLichtblau That opening sentence is a very general remark, not intended to be an explanation of the issue here. My intention was to demonstrate a pedestrian approach providing the correct answer, while the introduction says more or less: "Working with Integrate you have to be careful". I hope you haven't taken that statement as a critics of the internal implementation of Integrate, being one of my favourite functions in Mathematica as can be seen from my activity on this forum. $\endgroup$ Commented Feb 28, 2020 at 18:08
0
$\begingroup$

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` *)
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.