Here comes a possible solution using (directly) the numerical_integral
. (Without getting errors.)
First, let us introduce the function
h(x,y)=12π⋅(1+|1+x+y|3)exp−x2+y22 .
We will immediately also compute numerically the integral
∫20−20h(0,y)dy ,
because the (type for the) sage result is possibly unexpected.
sage: h(x,y) = exp( -(x^2+y^2)/2 ) / 2 / pi / ( 1+abs(x+y+1)^3 )
sage: numerical_integral( lambda y : h(0,y), (-20,20) )
(0.2032706496721935, 2.562242049607144e-09)
sage: numerical_integral( lambda y : h(0,y), (-20,20) )[0]
0.2032706496721935
So the expression numerical_integral( lambda y : h(0,y), (-20,20) )
is not a float number, but a tuple containing an approximative value for the integral, and the error made during the computation.
So a possible way to get the integral is...
sage: h(x,y) = exp( -(x^2+y^2)/2 ) / 2 / pi / ( 1+abs(x+y+1)^3 )
sage: numerical_integral( lambda x : numerical_integral( lambda y : h(x,y), (-20,20) )[0], (-20,20) )
(0.45483611545273767, 1.9556993802183342e-09)
This is in concordance with gp/pari, for instance:
? \p 50
realprecision = 57 significant digits (50 digits displayed)
? h(x,y) = exp( -(x^2+y^2)/2 ) / 2 / Pi / ( 1 + abs(x+y+1.0)^3 )
%1 = (x,y)->exp(-(x^2+y^2)/2)/2/Pi/(1+abs(x+y+1.0)^3)
? intnum( x=-20,20, intnum( y=-20,20, h(x,y) ) )
%2 = 0.45483723187162310868693510119720179912113863975808
Note: The given integral is roughly the expected value of the random variable
11+|1+X+Y|3 ,
where X,Y are two independent standard normal random variables. (Distributed following the N(0,12) law...)
(We integrate on the full real line instead of [−20,+20], the exponential decay is already strong.)
We can of course change variables using instead Z,W=1√2(X±Y). (Use the plus for Z, the minus for W.) The corresponding passage matrix is orthogonal. So the integral should be the same as
E[11+|1+X+Y|3]=E[11+|1+Z√2|3]=∫R11+|1+z√2|31√2πexp−z22dz .
One more numerical check:
sage: numerical_integral( 1/(1+abs(1+z*sqrt(2))^3) * exp(-z^2/2), [-20,20] )[0] / sqrt(2*pi).n()
0.454836900304870
sage: numerical_integral( 1/(1+abs(1+z*sqrt(2))^3) * exp(-z^2/2) / sqrt(2*pi), [-20,20] )
(0.4548369003048696, 5.4312008620095026e-08)