Ask Your Question
2

Solving numerically a multiple integral

asked 2016-04-03 19:08:42 +0200

b0ira gravatar image

updated 2016-04-03 20:26:06 +0200

vdelecroix gravatar image
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

Could anyone help me please?

edit retag flag offensive close merge delete

Comments

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

kcrisman gravatar imagekcrisman ( 2016-04-05 03:06:22 +0200 )edit

1 Answer

Sort by ยป oldest newest most voted
1

answered 2016-07-28 14:09:48 +0200

mforets gravatar image

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.

edit flag offensive delete link more

Comments

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.

kcrisman gravatar imagekcrisman ( 2016-07-28 14:40:22 +0200 )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).

mforets gravatar imagemforets ( 2016-07-29 14:26:15 +0200 )edit

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: 2016-04-03 19:08:42 +0200

Seen: 334 times

Last updated: Apr 03 '16