Numerical Integration Results

The numerical integration of (4) using a low degree Gauss rule is very accurate when $\vert\rho\vert \simeq 1$, except when $h$ is close to but not equal to $sk$. In order to avoid the problems when $h\simeq sk$, Drezner and Wesolowsky rewrote (4) as

$$L(h,k,\rho) = L(h,k,s) - \frac{s}{2\pi} \int_0^{\sqrt{1-\rho... ...}}e^{-\frac{(h-sk)^2}{2x^2}}e^{-\frac{shk}{1+\sqrt{1-x^2}}}dx.$$

The integrand in this expression is smooth except for the term $e^{-\frac{(h-sk)^2}{2x^2}}$ that is the source of relatively large errors when $h$ is close to $sk$. Using the Taylor expansion

$$\frac{1}{\sqrt{1-x^2}} e^{-\frac{shk}{1+\sqrt{1-x^2}}}= e^{-\frac{shk}{2}}(1+(4-shk)x^2/8)+O(x^4),$$ (5)

$L(h,k,\rho)$ becomes

$$L(h,k,\rho) = L(h,k,s) - \frac{s}{2\pi} \int_0^{\sqrt{1-\rho^2}} e^{-\frac{(h-sk)^2}{2x^2}}e^{-\frac{shk}{2}}(1+(4-shk)x^2/8)dx$$ (6)

$$- \frac{s}{2\pi} \int_0^{\sqrt{1-\rho^2}} e^{-\frac{(h-sk)^2}{2x^2}... ...k}{1+\sqrt{1-x^2}}}}{\sqrt{1-x^2}}-e^{-\frac{shk}{2}}(1+(4-shk)x^2/8)\bigg) dx.$$

The Drezner and Wesolowsky algorithm computes the second integral numerically, and the integral of the term $e^{-\frac{(h-sk)^2}{2x^2}}e^{-\frac{shk}{2}}(1+(4-shk)x^2/8)$ analytically, using the formula

$$\int_0^{a}e^{-b^2/(2x^2)}(1+cx^2)dx=a(1-c(b^2-a^2)/3) e^{-b^2/2a^2}-b(1-cb^2/3)\sqrt{2\pi}\Phi(-b/a),$$

with $a=\sqrt{1-\rho^2}$, $b=\vert h-sk\vert$ and $c=(4-shk)/8$.