This question would have been more properly posed on ask.sagemath.org...
The firs argument of numerical integral must be a function of one (real)
variable returning a real. It is indeed possible to nest such integrals,
but the second one *nests* the first one, hence long computation times.
Examples:
sage: (e^-(x^2/2)).integrate(x,-oo,oo)
sqrt(2)*sqrt(pi)
sage: (e^-(x^2/2)).integrate(x,-oo,oo).n()
2.50662827463100
sage: %time numerical_integral(lambda x1:e^-(x1^2/2),-oo,oo)
CPU times: user 29.8 ms, sys: 21 µs, total: 29.9 ms
Wall time: 29.8 ms
(2.5066282746309994, 2.5512942303338327e-08)
Both precision and computation time are acceptable. But consider:
sage: var("y")
y
sage: (e^-((x^2+y^2)/2)).integrate(y,-oo,oo).integrate(x,-oo,oo)
2*pi
sage: (e^-((x^2+y^2)/2)).integrate(y,-oo,oo).integrate(x,-oo,oo).n()
6.28318530717959
sage: %time numerical_integral(lambda x2:numerical_integral(lambda
x1:e^-((x1^2+x
....: 2^2)/2),-oo,oo)[0],-oo,oo)
CPU times: user 4.05 s, sys: 32 ms, total: 4.08 s
Wall time: 4.08 s
(6.283185307182902, 6.39626298408729e-08)
Acceptable precision, coputation time questionable. And this is only for a
double integral of a "well behaved" real function.
Furthermore, as pointed out, the values of your function are complex. You
must, at each step, separate real and imaginary part...
Therefore, I doubt that this solution is viable.
Monte-Carlo integration has been suggested. This solution has been
well-explored in some R packages dedicated to Bayesian statistics, but the
computation of the =value of the integral is still a delicate question
(and, being only of secondary interest to applied statistics, not as well
explored).
You may try to separate the real and imaginary part of your function, and
sample from these part. Look at Stan <https://mc-stan.org/> and the
companion R package bridgesampling
<https://cran.r-project.org/web/packages/bridgesampling/index.html> (which
could compute the integral from a pseudo-random sample obtained via Stan).
HTH,
Le jeudi 27 février 2020 04:21:21 UTC+1, saad khalid a écrit :
>
> I'm trying to compute/estimate a rather complicated looking integral. Here
> is the code I'm trying to run, with the necessary constants defined:
>
> var('t1,t2,u,w,k')
> T = 1
> m = 100
> E = 1
> v = 0
> y=1
> O = 1
> integral(integral(integral(
> integral(integral(
> e^(-t1^2/T^2)*e^(-t2^2/T^2)*e^(I*O*t1)*
> e^(-I*O*t2)*e^(-I*E*y^2*(1 - v)*t1^2/2)*
> e^(-I*E*y^2*(1 - v)*t2^2/2)*e^(-I*k*y*(1 - u)*t1)*
> e^(I*k*y*(1 - v)*t2)*
> e^((1 + I)*(sqrt(E)*y*w*t1 + w*k/sqrt(E)))*
> e^((1 - I)*(sqrt(E)*y*u*t2 + u*k/sqrt(E)))*
> e^(-w^2/2)*e^(-u^2/2)*w^(-1/2 + I*m^2/(2*E))*
> u^(-1/2 - I*m^2/(2*E)), (u, 0, Infinity)), (w, 0,
> Infinity)), (t2, -Infinity, Infinity)), (t1, -Infinity,
> Infinity)), (k, -Infinity, Infinity))
>
> I haven't been able to get a result from this code, it seems to run
> forever. I was hoping to be able to estimate the integral with some
> numerical methods, however I was having trouble getting a numerical
> integral set up properly. My first question is, can someone help me set up
> multivariable numerical integrals properly. I was trying something like
> numerical_integral(x*y,(x,0,1),(y,0,1))
> or
> numerical_integral(numerical_integral(x*y,(x,0,1)),(y,0,1))
>
> but neither seem to be the correct format, as they both give errors.
>
> My second question is, can anyone give some advice on how to approximate
> such an integral, where it's multivariable and the bounds are at infinity?
> I don't really know where to start.
>
> Thanks!
>
>
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-support/29ba4111-0ccc-4143-b68a-092e29ec6ceb%40googlegroups.com.
