Ask Your Question

2D Gaussian integrals

asked 2023-05-29 12:18:37 +0200

Unwise gravatar image

updated 2023-05-29 12:36:32 +0200

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)

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 flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2023-05-30 01:53:27 +0200

achrzesz gravatar image

Your sage code can be replaced by:

var('a b')
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
edit flag offensive delete link more

Your Answer

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

Add Answer

Question Tools

1 follower


Asked: 2023-05-29 12:18:37 +0200

Seen: 74 times

Last updated: May 30