# 2D Gaussian integrals

Hello. I have the following Gaussian function $$G(\mathbf r) = \frac{1}{\sqrt \pi R_0}e^{-\frac{r^2}{2R_0^2}} , \qquad R_0>0$$ of a 2D vector $\mathbf r=(x, y)$ as well as a 'transfer' function $$T_k(\mathbf r, \mathbf r', d)=\frac{e^{ikd+\frac{ik}{2d}|\mathbf r-\mathbf r'|^2}}{id(2\pi/k)}, \qquad d,k>0.$$ My objective is to performe the symbolic integration $$G'(\mathbf r_1)=\int d^2\mathbf r_0G(\mathbf r_0) T_k(\mathbf r)$$ in Sage. First of all: this integral is actually extremely simple to compute by hand, by applying the standard Gaussian formula. However, I'm using this example to learn how to do symbolic integration in Sage.

My first question: is there some particular syntax to do this 2D integral directly, or is it always necessary to do the integrals in $dx$ and $dy$ separately? Following this route, after expanding the modulus squared we get $$G'(\mathbf r_1) = \frac{1}{\sqrt pi R_0}\frac{e^{ikd}e^{\frac{ik}{2d}(x_1^2+y_1^2)}}{id(2\pi/k)}\int dx_0 \ e^{-\left(\frac{1}{2R_0^2}-\frac{il}{2d}\right)x_0^2-\frac{ik}{d}x_0x_1}\int dy_1 \ e^{-\left(\frac{1}{2R_0^2}-\frac{il}{2d}\right)y_0^2-\frac{ik}{d}y_0y_1}.$$ Since the two integrals are analogous, I focus on the first. The code I've written to try and solve it is the following:

#define all variables
var('x_0, x_1, y_0, y_1, d, k, R_0')
assume(d>0, k>0, R_0>0)

#define constant prefactors
a = 1/(2*R_0^2) - i*k/(2*d)
b = i*k*x_1 / d

#define integrand
I(x_0) = exp(-a*x_0^2 + b*x_0)

#perform Gaussian integral
result = I.integral(x_0, -oo, oo)
show(result)


However, this only returns several warnings of the type

 Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.


on top of the error Error trying to find limit of -sqrt(pi)*(-2*i)/(2*i*sageVARR_0^4*sageVARd*sageVARk/(-2*sageVARR_0^2*sageVARd^2+2*sqrt(sageVARR_0^8*sageVARd^2*sageVARk^2+sageVARR_0^4*sageVARd^4))+1)/sqrt(-sageVARR_0^2*sageVARd^2+sageVARR_0^2*sageVARd*sqrt(sageVARR_0^4*sageVARk^2+sageVARd^2))*sageVARR_0^2*sageVARd*exp((-i)*sageVARR_0^2*sageVARk^2*sageVARx_1^2/(2*sageVARR_0^2*sageVARd*sageVARk+2*i*sageVARd^2))/2*erf(-1/(-2*i)*(2*i*sageVARR_0^4*sageVARd*sageVARk/(-2*sageVARR_0^2*sageVARd^2+2*sqrt(sageVARR_0^8*sageVARd^2*sageVARk^2+sageVARR_0^4*sageVARd^4))+1)*sqrt(-sageVARR_0^2*sageVARd^2+sageVARR_0^2*sageVARd*sqrt(sageVARR_0^4*sageVARk^2+sageVARd^2))/sageVARR_0^2/sageVARd*(sageVARx_0+i*sageVARk*sageVARx_1/sageVARd/(i*sageVARR_0^2*sageVARk-sageVARd)*2*sageVARR_0^2*sageVARd/2)).

With a finite choice of extremal points, the integration is performed correctly. I assume the problem here is that Sage has trouble finding the limits of a complex functions, but how do I circumvent the issue?

edit retag close merge delete

Sort by » oldest newest most voted

Your sage code can be replaced by:

var('a b')
assume(a>0)
I0=integrate(exp((-a*x^2+b*x)),x,-oo,oo)
var('x_0, x_1, y_0, y_1, d, k, R_0')
I1=I0.subs({a:1/(2*R_0^2) - i*k/(2*d),b:i*k*x_1 / d})
%display latex
I1

more