# Plot of 3D function defined through NIntegrate over infinite intervals

I want to speed-up the plotting of the function like this

rad[s_] := 1/(1 + s^4);
func[x_?NumericQ, y_?NumericQ] := NIntegrate[rad[Norm[{u, v}]], {u, x, \[Infinity]}, {v, y, \[Infinity]}, PrecisionGoal -> 5, AccuracyGoal -> 3, Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}];
Plot3D[func[x, y], {x, -5, 5}, {y, -5, 5}] // AbsoluteTiming


that currently costs about 8.6s on my quite fast laptop.

I read about the trick of reducing this to plotting of a numerical solution to a differential equation, but it worked for one variable and finite interval of integration. Here it would be a PDE

DSolve[D[z[x, y], x, y] == rad[Norm[{x, y}]], z, {x, y}]


but I need NDSolve and initial conditions at infinity which are unclear for me.

Dmitri

Here an other way, that is also faster.

Construct an interpolation function of the -5 < u < 5 and -5 < v < 5 area and integrate the remaining 3 stripes only once.

(rad[s_] = 1/(1 + s^4);

f55 = NIntegrate[rad[Norm[{u, v}]], {u, 5, \[Infinity]}, {v, 5, \[Infinity]}];


Edit: Adding the two missing stripe areas

fu5inf = NIntegrate[rad[Norm[{u, v}]], {u, 5, \[Infinity]}, {v, -5, 5}];

fv5inf = NIntegrate[rad[Norm[{u, v}]], {u, -5, 5}, {v, 5, \[Infinity]}];

nint[a_?NumericQ, b_?NumericQ] := NIntegrate[rad[Norm[{u, v}]], {u, a, 5}, {v, b, 5}];

fi = FunctionInterpolation[nint[a, b] + f55 + fu5inf + fv5inf, {a, -5, 5}, {b, -5, 5}];

Plot3D[fi[x, y], {x, -5, 5}, {y, -5, 5}, ImageSize -> 400]) // Timing • I save a little time by adding integrand = Simplify[rad[Norm[{u, v}]], u ∈ Reals && v ∈ Reals]; and using it as the integrand. – Michael E2 May 20 '18 at 15:48