Thank you for reading my questions.
It is known that the stationary phase method has this form: $$ \int_a^bg(t)e^{jf(t)dt}\approx\sum\limits_{t_0\in\Sigma}g(t_0)e^{jf(t_0)+j*\text{sign}(f''(t_0))\frac{\pi}{4}}\frac{\sqrt{2\pi}}{\sqrt{|f''(t_0)|}} \tag{1} $$ You can check the wiki for this: https://en.wikipedia.org/wiki/Stationary_phase_approximation
This works well when the zero point is at bounded numbers. But what if the zero point is at infinity?
Here is the integration I want to deal with
$$ \int_{-\infty}^{\infty} \frac{\beta_0 e^{jk \left( \sqrt{x_1^2+y_1^2 - 2xy_1 + x^2} - \sqrt{x_2^2+y_2^2 - 2xy_2 + x^2} \right)}}{\sqrt{x_1^2+y_1^2 - 2xy_1 + x^2}\sqrt{x_2^2+y_2^2 - 2xy_2 + x^2}} dx \tag{2} $$
where $x_1 = x_2, y_1\neq y_2$
Then $$ \begin{align} g(x)=&\frac{\beta_0}{\sqrt{x_1^2+y_1^2 - 2xy_1 + x^2}\sqrt{x_2^2+y_2^2 - 2xy_2 + x^2}}\\ \overset{M\rightarrow\infty}{\approx}&\frac{\beta_0}{x^2} \end{align} $$ $$ f(x)=k \left( \sqrt{x_1^2+y_1^2 - 2xy_1 + x^2} - \sqrt{x_2^2+y_2^2 - 2xy_2 + x^2} \right) $$ $$ f'(x)=\frac{{k (x - y_1)}}{{\sqrt{x_1^2+y_1^2 + x^2 - 2 y_1 x }}} - \frac{{k (x - y_2)}}{{\sqrt{x_2^2+y_2^2 + x^2 - 2 y_2 x}}} $$ $$ \begin{align} f''(x)=&kx_1^2\left(\frac{1}{(x_1^2+y_1^2 + x^2 - 2y_1x)^{\frac{3}{2}}} - \frac{1}{(x_2^2+y_2^2 + x^2 - 2y_2x)^{\frac{3}{2}}}\right) \end{align} $$ $f''(x)$ is equal to zero iff $y_1^2+x^2-2y_1 x=y_2^2+x^2-2y_2 x$, which is actually $x=\frac{y_1+y_2}{2}$.
And it can be proved that $\lim_{x\rightarrow \infty} f'(x)=0$ and $\lim_{x\rightarrow -\infty} f'(x)=0$. Hence, there are only two zero points for $f'(x)$ and both of them are at infinity. Hence, I use the following stationary phase method: $$ \lim_{x\rightarrow \infty}g(x)e^{jf(x)+j*\text{sign}(f''(x))\frac{\pi}{4}}\frac{\sqrt{2\pi}}{\sqrt{|f''(x)|}}+\lim_{x\rightarrow -\infty}g(x)e^{jf(x)+j*\text{sign}(f''(x))\frac{\pi}{4}}\frac{\sqrt{2\pi}}{\sqrt{|f''(x)|}} \tag{3} $$ We need the limits below $$ \begin{align} \lim_{x\rightarrow \infty}\frac{g(x)}{\sqrt{|f''{(x)}|}}=&\lim_{x\rightarrow \infty}\sqrt{\frac{g^2(x)}{|f''{(x)}|}}\\ =&\lim_{x\rightarrow \infty}\sqrt{\frac{\beta_0^2}{kx_1^2\left|\frac{x^4}{(x_1^2+y_1^2 + x^2 - 2xy_1)^{\frac{3}{2}}} - \frac{x^4}{(x_2^2+y_2^2+ x^2 - 2xy_2)^{\frac{3}{2}}}\right|}}\\ =&\sqrt{\frac{\beta_0^2}{3kx_1^2|y_1-y_2|}} \end{align} $$
\begin{align} \lim_{x\rightarrow \infty}\text{sign}(f''(x))=&\lim_{x\rightarrow \infty}\text{sign}(x_2^2+y_2^2-x_1^2-y_1^2+2x(y_1-y_2))\\ =&\text{sign}(y_1-y_2) \end{align}
\begin{align} \lim_{x\rightarrow\infty}f(x)=&\lim_{x\rightarrow \infty}k \left( \sqrt{x_1^2+y_1^2 - 2xy_1 + x^2} - \sqrt{x_2^2+y_2^2 - 2xy_2 + x^2} \right)\\ =&k(y_2-y_1) \end{align}
By substitution in (3), we have (1) is equal to $$ 2\sqrt{\frac{2\pi\beta_0^2}{3kx_1^2|y_1-y_2|}}\cos{(k(y_2-y_1)+\text{sign}(y_1-y_2)\frac{\pi}{4})} \tag{4} $$ I write a code to compare (2) and (4)
clc;
clear;
close all;
c = 3e8; %speed of light
frequency = 10e9; % The wavelength of the original paper is unknown
lambda = c/frequency;
x1 = 1.5;
x2 = x1;
y1 = -1;
y2 = 1;
k = 2*pi/lambda;
beta0 = 1;
c = 3e8; %speed of light
frequency = 100e9; % The wavelength of the original paper is unknown
lambda = c/frequency;
f = @(x) k * (sqrt((x1^2+y1^2) - 2*x*y1 + x.^2) ...
- sqrt(x2^2+y2^2 - 2*x*y2 + x.^2));
g = @(x) beta0 ./ sqrt(x1^2+y1^2 - 2.*x.*y1 + x.^2)...
./ sqrt(x2^2+y2^2 - 2.*x.*y2 + x.^2);
fdd = @(x) k * x1^2/ (x1^2+y1^2 + x^2 - 2 * y1*x)^(3/2) - ...
k * x2^2/ (x2^2+y2^2 + x^2 - 2 * y2*x)^(3/2);
integrand = @(x) g(x).*exp(1j*f(x));
x = 1:1e6:1e8;
integrate_result = zeros(1,length(x));
SPM_result = integrate_result;
for i = 1:length(x)
integrate_result(i) = abs(integral(integrand, -x(i), x(i)));
SPM_result(i) = g(x(i))*exp(1j*f(x(i))+1j*sign(fdd(x(i)))*pi/4)*sqrt(2*pi)/sqrt(abs(fdd(x(i))))+...
g(-x(i))*exp(1j*f(-x(i))+1j*sign(fdd(-x(i)))*pi/4)*sqrt(2*pi)/sqrt(abs(fdd(-x(i))));
disp(i);
end
plot(x,abs(integrate_result),'-^','markerindices',1:ceil(length(x)/10):length(x));
hold on;
plot(x,abs(SPM_result),'-+','markerindices',1:ceil(length(x)/10):length(x));
legend('Integration','Analytical limit')
title('Comparison between limit and the simulation');
disp(['The ratio square:',num2str(abs(SPM_result(end)/integrate_result(end))^2)]);
I find they are not equal but their ratio is always $\frac{4}{3}$. When I change the parameters, the ratio stays the same.
I guess substituting $x$ by $\frac{1}{x}$ won't make a difference. I tried it. The ratio then is $4$.
clc;
clear;
close all;
c = 3e8; %speed of light
frequency = 10e9; % The wavelength of the original paper is unknown
lambda = c/frequency;
x1 = 1.5;
x2 = x1;
y1 = -1.2;
y2 = 5;
k = 2*pi/lambda;
beta0 = 1;
c = 3e8; %speed of light
frequency = 100e9; % The wavelength of the original paper is unknown
lambda = c/frequency;
f_orig = @(x) k * (sqrt((x1^2+y1^2) - 2*x*y1 + x.^2) ...
- sqrt(x2^2+y2^2 - 2*x*y2 + x.^2));
g_orig = @(x) beta0 ./ sqrt(x1^2+y1^2 - 2.*x.*y1 + x.^2)...
./ sqrt(x2^2+y2^2 - 2.*x.*y2 + x.^2);
fdd_orig = @(x) k * x1^2/ (x1^2+y1^2 + x^2 - 2 * y1*x)^(3/2) - ...
k * x2^2/ (x2^2+y2^2 + x^2 - 2 * y2*x)^(3/2);
f_recp = @(x) k * (sqrt((x1^2+y1^2) - 2./x*y1 + 1./x.^2) ...
- sqrt(x2^2+y2^2 - 2./x*y2 + 1./x.^2));
g_recp = @(x) -beta0 ./ sqrt(x1^2+y1^2 - 2./x.*y1 +1./ x.^2)...
./ sqrt(x2^2+y2^2 - 2./x.*y2 + 1./x.^2)./x.^2;%We have -1./x.^2 due to d(1/x)
fdd_recp = @(x) k * ((2+3*(x1^2+y1^2)*x^2+x*y1*(-6-2*(x1^2+y1^2)*x^2+3*y1*x))/ ...
(x^6 * ((x1^2+y1^2) + 1/x^2 - 2*y1/x)^(3/2)) + ...
(-2-3*(x2^2+y2^2)*x^2+x*y2*(6+2*(x2^2+y2^2)*x^2-3*x*y2))/ ...
(x^6 * ((x2^2+y2^2) + 1/x^2 - 2*y2/x)^(3/2)));
integrand = @(x) g_orig(x).*exp(1j*f_orig(x));
x = 1:1e3:1e5;
integrate_result = zeros(1,length(x));
SPM_result = integrate_result;
for i = 1:length(x)
integrate_result(i) = abs(integral(integrand, -x(i), x(i)));
x_r = 1./x(i);
SPM_result(i) = g_recp(x_r)*exp(1j*f_recp(x_r)+1j*sign(fdd_recp(x_r))*pi/4)*sqrt(2*pi)/sqrt(abs(fdd_recp(x_r)))+...
g_recp(-x_r)*exp(1j*f_recp(-x_r)+1j*sign(fdd_recp(-x_r))*pi/4)*sqrt(2*pi)/sqrt(abs(fdd_recp(-x_r)));
disp(i);
end
plot(x,abs(integrate_result),'-^','markerindices',1:ceil(length(x)/10):length(x));
hold on;
plot(x,abs(SPM_result),'-+','markerindices',1:ceil(length(x)/10):length(x));
legend('Integration','Analytical limit')
title('Comparison between limit and the simulation');
disp(['The ratio square:',num2str(abs(SPM_result(end)/integrate_result(end))^2)]);
Is it just a coincidence? I checked my deduction many times but I cannot find the error. Could you tell me if is this the right way to use the stationary phase method at an infinite zero point?
I think you can focus on the case where zero points are at infinity because I have given the formulas. But feel free if you think it is interesting to substitute $x$ by $\frac{1}{x}$ since I may make some mistakes due to my ignorance.
Thank you! Any help is appreciated.