# Solving numerically a multiple integral

var("r1, r2, v1, v2")
tau(r1,r2,v1,v2)=((r1-r2)-0.8)/(v1-v2)
E(r1,r2,v1,v2)= exp(-tau(r1,r2,v1,v2)/3)*1.5/tau(r1,r2,v1,v2)^2
f(r1,r2,v1,v2)=exp(-(v1^2+v2^2)/2)*exp(-E(r1,r2,v1,v2)/2)

integral(integral(integral(integral(f(r1,r2,v1,v2),r1,0,10),r2,0,10),v1,-inf,inf),v2,-inf,inf)


RuntimeError: ECL says: In function GCD, the value of the second argument is 1.0 which is not of the expected type INTEGER

edit retag close merge delete

Most likely related to http://trac.sagemath.org/ticket/14821 and keepfloat:true in Maxima

( 2016-04-05 03:06:22 +0100 )edit

Sort by » oldest newest most voted

Have you tried scipy's nquad function? I think that in your problem special care has to be taken in the line $v_1=v_2$, but it can be handled with the opts parameter.

more

Can you give specific syntax? (Maxima and GSL both have numerical integration that Sage syntax automatically supports, you could give all three!) Not sure how these support infinite bounds though.

( 2016-07-28 14:40:22 +0100 )edit

The infinite bounds can be passed as [-np.inf, np.inf]. However I didn't manage to obtain an answer for OP's integral -- it takes ages to converge. Perhaps a Monte Carlo integration is more reliable in this case (see mcint).

( 2016-07-29 14:26:15 +0100 )edit