How can I speed up NIntegrate computation?

I want to Plot a definite integral

myIntegral[x_] := NIntegrate[
Sqrt[(a - b)^2 + (c - d)^2]/(a - b)
, {a, 0, x}
, {b, x, 1}
, {c, 0, 1}
, {d, 0, 1}
]

Plot[myIntegral[x], {x, 0, 1}]


I am using Mathematica 11.3. When I evaluate the code, Mathematica immediately gives some warnings about numerical integration converging too slowly, then it seems stuck in computation. I have tried some options of NIntegrate, for example as suggested here, but have had no success.

Is there a way to obtain the result of Plot in an acceptable time (e.g. some minutes)?

• Method -> "MonteCarlo" allows you to plot the function in less than a minute, but the plot is not perfectly smooth: i.imgur.com/sEuL6oN.png increasing the number of integration points makes the plot slower but nicer. Jun 19 '18 at 14:56
• Thank you for the hint. It seems that Method -> "QuasiMonteCarlo" can give a smooth plot in this case. In addition, it is slightly faster (few seconds less). Jun 19 '18 at 16:15
• @renato Self-answers are encouraged on StackExchange, so you would be welcome to post your comment as an answer! Jun 19 '18 at 17:46

The integral is within reach of Integrate

Assuming[0 < x < 1,
Simplify[
Integrate[Sqrt[(a-b)^2 + (c-d)^2]/(a-b), {a,0,x}, {b,x,1}, {c,0,1}, {d,0,1}]]]


which gives you

(-11*Sqrt[2] + 11*Sqrt[2+(-2+x)*x] + x*(-8-7*Sqrt[2+(-2+x)*x] + 13*Sqrt[1+x^2] +
2*x*(6-3*Sqrt[2+(-2+x)*x] + x*(-4+2*x+Sqrt[2+(-2+x)*x] - Sqrt[1+x^2]))) +
18*x^2*ArcCoth[Sqrt[2+(-2+x)*x]] + 9*ArcSinh[1] - 3*ArcSinh[1-x] -
3*ArcSinh[x] + 6*Log[1-x] + 9*x*Log[-1+Sqrt[2+(-2+x)*x]] + 3*((-2+x)*Log[1 +
Sqrt[2 + (-2 + x)*x]] - 4*x*Log[-(((-1 + x)*(1 + Sqrt[1 + x^2]))/x)] +
x^3*(Log[2+x^2+2*Sqrt[1+x^2]] + 2*Log[(1-x)/(x+x*Sqrt[2-2*x+x^2])])))/36


and that will be far faster than doing thousands of NIntegrate inside Plot

You can first calculate the integral over c and d to get:

res = Integrate[{Sqrt[t + (c - d)^2]}, {c, 0, 1}, {d, 0, 1}];
out=FullSimplify[res, {a >= 0, b >= 0, t > 0}]


{1/3 (2 t^(3/2) + Sqrt[1 + t] - 2 t Sqrt[1 + t] + 3 t ArcCoth[Sqrt[1 + t]])}

Now put (a-b)^2 for t and dividing by (a-b).

out2=FullSimplify[(out /. t -> (a - b)^2)/(a - b), {a >= 0,b >= 0, t >= 0}];


So you can then write a simplified version of you integral as:

myIntegral[x_] := NIntegrate[out2, {a, 0, x}, {b, x, 1}]

Plot[myIntegral[x], {x, 0, 1}]


However, to get the first integral i put in some numbers for t=10,100,1000... and then guessed the solution. That was a bit faster, but mathematica is also able to find this solution if you wait for a while.