Loading [MathJax]/jax/output/HTML-CSS/jax.js
Ask Your Question
2

2D integrate fails.

asked 11 years ago

anonymous user

Anonymous

updated 10 years ago

FrédéricC gravatar image

Hi Guys: I have a function as bellow:

sage: h_1_1_1
0.500000000000000*e^(-1/2*p1^2 - 1/2*p2^2)/((abs(p1 + p2 + 1)^3 + 1)*pi)

I need to performa an integrate as bellow:

h_1_1_1.integrate((p1, -20, 20)).integrate((p2, -20, 20))

Got followings:

0.159154943092*integrate(integrate(e^(-1/2*p1^2 - 1/2*p2^2)/(abs(p1 + p2 + 1)^3 + 1), p1, -20, 20), p2, -20, 20)

Then i need to know the value, using following

N(_)

Then got errors: .

    199         # We only return the result.
--> 200         return numerical_integral(f, a, b)[0]
    201 
    202     def _tderivative_(self, f, x, a, b, diff_param=None):

/Volumes/sage-5.9-OSX-64bit-10.7-x86_64-Darwin/Sage-5.9-OSX-64bit-10.7.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/gsl/integration.so in sage.gsl.integration.numerical_integral (sage/gsl/integration.c:1802)()

ValueError: Integrand has wrong number of parameters
Preview: (hide)

1 Answer

Sort by » oldest newest most voted
2

answered 7 years ago

dan_fulea gravatar image

updated 7 years ago

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)expx2+y22 . We will immediately also compute numerically the integral 2020h(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=12(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+Z2|3]=R11+|1+z2|312πexpz22dz . 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)
Preview: (hide)
link

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 11 years ago

Seen: 461 times

Last updated: Jul 13 '17