From a steepest descent approach:
Recast the integral into the following form:
$$ \int\limits_{0}^{\infty} cos(x^2)dx = Re\int\limits_{0}^{\infty}e^{ix^2}dx$$
From here on, I'll drop the Real operator and it will be implicit that you will take the real part at the end. Performing a change of variables and making x a complex variable, the above integral can be recast in the following format:
$$\sqrt{s}\int\limits_{0}^{1}e^{isz^2}dz$$
which is in the expected form for a steepest descent method:
$$g(z)\int e^{sf(z)}dz$$
with $g(z) = \sqrt{s}$ and $f(z) = iz^2$. It can be shown that the path of steepest descent cuts through the origin at an angle of $\frac{\pi}{4}$ degrees. Using Arfken and Weber's notation, $z_0=0$ and $\alpha = \frac{\pi}{4}$. Now the tricky part is drawing the contour. This is where I believe the OP made an error. In order to draw a contour that crosses the origin at $\pi/4$, part of that contour would have to come from the bottom left quadrant. However we also need the contour to run along the real axis from 1 to the origin. You can try it for yourself, you'll find drawing this contour would be very difficult. One way to get around this is to take only 1/2 of the path of steepest descent. In other words when drawing your contour start at the origin then proceed in the $\pi/4$ direction rather than start in the bottom left quadrant and move to the top right. This leads to the OP's missing factor of 1/2. Once you've done this, you can close the contour by coming back to the origin from 1. This leads you to the correct answer:
$$\int\limits_{C}(\cdot)= \frac{1}{2}\int\limits_{S.D.}(\cdot) - \int\limits_{0}^{1}(\cdot) = 0$$
$$ \sqrt{s}\int\limits_{0}^{1}e^{isz^2}dz = \frac{1}{2}\frac{\sqrt{2\pi}g(z_0)e^{sf(z_0)}e^{i\alpha}}{|sf''(z_0)|^{1/2}}$$
plugging in the values from earlier and taking the real part, you should get the correct answer of $$\sqrt{\frac{\pi}{8}}$$
OP is basically correct. Recall that the Euler function is
$$ \varphi(\tau)~:=~\prod_{m\in\mathbb{N}}(1-q^m)=\frac{1}{\sum_{n\in\mathbb{N}_0}p(n)q^n},
\qquad q~\equiv~e^{2\pi i \tau},$$
and that the Dedekind eta function is
$$ \eta(\tau)~:=~q^{\frac{1}{24}}\varphi(\tau). $$
We calculate the partition function
$$\begin{align} p(n)~=~&\oint_{|q|<1}\frac{dq}{2\pi i}\! \frac{q^{-(n+1)}}{\varphi(\tau)}\cr
~=~&\int_{i\delta+\left[-\frac{1}{2},\frac{1}{2}\right]} \! d\tau~ \frac{q^{\frac{1}{24}-n}}{\eta(\tau)},
\qquad \delta~>~0, \cr
~=~&\int_{i\delta+\left[-\frac{1}{2},\frac{1}{2}\right]} \! d\tau~ \sqrt{-i\tau}\frac{q^{\frac{1}{24}-n}}{\eta(\tilde{\tau})}, \qquad
\tilde{\tau}~\equiv~-\tau^{-1}, \cr
~=~&P^{\prime}(n),&\end{align} $$
where
$$\begin{align} P(n)
~:=~&\int_{i\delta+\left[-\frac{1}{2},\frac{1}{2}\right]} \! \frac{d\tau}{2\pi}~ \frac{q^{\frac{1}{24}-n}\tilde{q}^{-\frac{1}{24}} }{\sqrt{-i\tau}\varphi(\tilde{\tau})}, \qquad
\tilde{q}~\equiv~e^{2\pi i \tilde{\tau}},\cr
~=~&\frac{1}{\sqrt{\lambda}}\int_{i\lambda\delta+\left[-\frac{\lambda}{2},\frac{\lambda}{2}\right]} \! \frac{d\bar{\tau}}{2\pi}~ \frac{e^{-\lambda S(\bar{\tau})}}{\sqrt{-i\bar{\tau}\rule[0pt]{0pt}{6pt}}\varphi(-\lambda\bar{\tau}^{-1})}, \qquad
\bar{\tau}~\equiv~\lambda\tau,
\end{align}$$
where
$$ \lambda~\equiv~\sqrt{n-\frac{1}{24}}~>~0 $$
is large, and where
$$\begin{align} S(\bar{\tau})~:=~&2\pi i\left(\bar{\tau}-\frac{1}{24\bar{\tau}}\right),\cr
S^{\prime}(\bar{\tau})~=~&2\pi i\left(\frac{1}{24\bar{\tau}^2}+1\right),\cr
S^{\prime \prime}(\bar{\tau})~=~&-\frac{\pi i}{6\bar{\tau}^3}.\end{align}$$
Note that
$$ \varphi(-\lambda\bar{\tau}^{-1}) ~\to ~1
\qquad \text{for}\qquad\lambda~\to~\infty.$$
There are 2 stationary points:
$$\begin{align}
\bar{\tau}_{\pm}~=~& \frac{\pm i}{2\sqrt{6}} \cr
\Updownarrow&\cr
\tau_{\pm}~=~& \frac{\pm i}{2\sqrt{6}\lambda} \cr
\Updownarrow&\cr
\tilde\tau_{\pm}~=~& \pm i 2\sqrt{6}\lambda. \cr
\end{align}$$
We calculate
$$\begin{align} S(\bar{\tau}_{\pm})~=~&\mp\pi \sqrt{\frac{2}{3}}, \cr
S^{\prime \prime}(\bar{\tau}_{\pm})~=~&\pm 8\sqrt{6}\pi.\end{align}$$
Clearly only the stationary point $\bar{\tau}_+$ in the upper halfplane contributes. The method of steepest descent yields
$$\begin{align} P(n)
~\sim~&\frac{1}{\lambda\sqrt{2\pi S^{\prime \prime}(\bar{\tau}_+)}} \frac{e^{-\lambda S(\bar{\tau}_+)}}{\sqrt{-i\bar{\tau}_+}} \cr
~=~&\frac{e^{\lambda \pi \sqrt{\frac{2}{3}}}}{\lambda2\pi\sqrt{2}}\qquad \text{for}\qquad\lambda~\to~\infty,
\end{align}$$
and therefore
$$ p(n)~=~P^{\prime}(n)
~\sim~\frac{e^{\lambda \pi \sqrt{\frac{2}{3}}}}{\lambda^2 4\sqrt{3}} \qquad \text{for}\qquad\lambda~\to~\infty.$$
For more details, see Ref. 1.
References:
- E.M. Stein & R. Shakarchi, Complex Analysis, 2003; Appendix A.4, Theorem 4.1, p. 334.
Best Answer
As @Maxim mentioned, the Laplas method can be applied to the integral, and the steepest descend is along the real axis, near $t=1$. The asymptotics can be obtained as the asymptotics of special function, because $I(x)=2K_{-n-1}(2x)$.
We can also get the asymptotics (few several terms) directly. Let $$I_n(x)=\int_{0}^{\infty}t^{n}e^{-x(t+\frac{1}{t})}dt$$ We can consider the generating function instead, namely $$J(\alpha,x)=\sum_{n=0}^\infty\frac{\alpha ^n}{n!}I_n(x)=\int_{0}^{\infty}e^{\alpha t}e^{-x(t+\frac{1}{t})}dt=\int_{0}^{\infty}e^{\alpha t}e^{-xf(t)}dt$$ where $f(t)=\frac{1}{t}+t$. Then $\,\,I_n(x)=\frac{\partial^n }{\partial\alpha^n}J(\alpha,x)\Big|_{\alpha=0}$
$f'(t)=0$ at $t=1$; therefore near this point $$f(t)=f(1)+f'(1)(t-1)+f''(1)\frac{(t-1)^2}{2!}+f'''(1)\frac{(t-1)^3}{3!}+...$$ $$f(t)=2+(t-1)^2-(t-1)^3+(t-1)^4-=...$$ $$J(\alpha,x)=e^{-2x}\int_{0}^{\infty}e^{\alpha t}e^{-x\big((t-1)^2-(t-1)^3+(t-1)^4-+...\big)}dt=e^{-2x}\int_{-1}^{\infty}e^{\alpha (t+1)}e^{-x\big(t^2-t^3+t^4-+...\big)}dt$$ $$=\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\sqrt x}^{\infty}e^{\frac{\alpha t}{\sqrt x}}e^{-t^2}e^{\frac{t^3}{\sqrt x}-\frac{t^4}{x}+ -...}dt=\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\infty}^{\infty}\,-\,\,\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\infty}^{-\sqrt x}$$ The second integral contains an additional exponentially small factor ($\sim e^{-x}$), so we drop it and focus on the first term. Expanding the integrand near $t=0$ into the series $$J(\alpha,x)=\frac{e^{-2x+\alpha}}{\sqrt x}\int_{-\infty}^{\infty}\Big(1+\frac{\alpha t}{\sqrt x}+\frac{\alpha^2 t^2}{2x}+...\Big)\Big(1+\frac{t^3}{\sqrt x}+\frac{t^6}{2x}-\frac{t^4}{x}+O\big(x^{-\frac{3}{2}}\big)\Big)e^{-t^2}dt$$
Using $\int_{-\infty}^{\infty}e^{-t^2}t^{2k}dt=(-1)^k\frac{\partial^k}{\partial s^k}\int_{-\infty}^{\infty}e^{-st^2}dt\Big|_{s=1}=(-1)^k\sqrt\pi\frac{\partial^k}{\partial s^k}\frac{1}{\sqrt s}\Big|_{s=1}$, integrating term by term and grouping the terms with the same power of $x$, we get $$J(\alpha,x)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\,e^{\alpha}\Big(1+\frac{\alpha^2+3\alpha+\frac{3}{4}}{4x}+O\big(\frac{1}{x^2}\big)\Big)$$ $$I_n(x)=\frac{\partial^n }{\partial\alpha^n}J(\alpha,x)\Big|_{\alpha=0}$$ $$I(0)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{3}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$ $$I(1)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{15}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$ $$I(2)= \sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{35}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$ $$....$$
$$I_n(x)=\sqrt\pi\,\,\frac{e^{-2x}}{\sqrt x}\Big(1+\frac{(2n+3)(2n+1)}{16x}+O\big(\frac{1}{x^2}\big)\Big)$$
The main asymptotic term of $I_n(x)$ does not depend on $n$.