Thanks to Maxim's comment for this answer. I've slightly rewritten it to try and do this without steepest descent, but the contour integration we'll need is essentially the same thing anyway.
The contribution missing from the singularity at $t=\pi/2$ is
$$K = \int_{\pi/2 - \epsilon}^{\pi/2}f_k(t)e^{ix\cos{t}} dt \approx \int_0^\epsilon (\frac{2t}{\pi})^k e^{ixt} dt$$
Integration by parts once and applying Riemann Lebesgue's lemma shows the $[\epsilon,\infty)$ contribution to be $O(1/x)$ (assuming $k<0$), which we will neglect.
So
$$K \sim \int_0^\infty(\frac{2t}{\pi})^k e^{ixt} dt$$
Following the same idea as a Fresnel integral you can deform the contour of integration to up the imaginary axis (show the quarter circle contribution vanishes) to get
$$K \sim i \int_0^\infty(\frac{2it}{\pi})^k e^{-xt} dt = (\frac{2}{\pi})^k i^{k+1} \frac{\Gamma(k+1)}{x^{k+1}}$$
(note $0<k+1<1$ so we were safe to neglect $O(1/x)$ terms). So all in all
$$I(x) \sim \cos{(x-\pi/4)} \sqrt{\frac{\pi}{2x}} + (\frac{2}{\pi})^k \frac{\Gamma(k+1)}{x^{k+1}} \cos{\frac{\pi(k+1)}{2}}$$
It is more convenient to consider $\displaystyle I(n,\alpha)=\int_0^1t^{-\alpha}e^{i\pi nt}dt$; then our integral is the areal part of $I(n,\alpha)$. Making the substitution $x=\frac{t}{n}$
$$I(n,\alpha)=\frac{1}{n^{1-\alpha}}\int_0^nx^{-\alpha}e^{i\pi x}dx=\frac{1}{n^{1-\alpha}}\int_0^\infty x^{-\alpha}e^{i\pi x}dx-\frac{1}{n^{1-\alpha}}\int_n^\infty x^{-\alpha}e^{i\pi x}dx=I_1+I_2$$
The first integral is well-known and can be evaluated, for example, by means of integration in the complex plane.
Let's consider the integral along the following contour:
$$\frac{1}{n^{1-\alpha}}\oint z^{-\alpha}e^{i\pi z}dz=\int_r^R z^{-\alpha}e^{i\pi z}dz+I_R+\int_R^r \big(e^\frac{\pi i}{2}z\big)^{-\alpha}e^{-\pi z}idz+I_r\tag{0}$$
There are no singularities inside the contour; therefore, $\displaystyle \oint=0$.
It is not difficult to show that the integrals along the big and small quarter-circles $\,\displaystyle I_r, I_R\to 0\,\,\text{at}\, r\to 0\,\, \text{and}\,\, R\to\infty\,$ correspondingly.
Therefore,
$$I_1=\frac{1}{n^{1-\alpha}}\int_0^\infty x^{-\alpha}e^{i\pi x}dx=\frac{ie^{-\frac{i\pi\alpha}{2}}}{n^{1-\alpha}}\int_0^\infty x^{-\alpha}e^{-\pi x}dx=\frac{e^\frac{i\pi(1-\alpha)}{2}}{n^{1-\alpha}}\int_0^\infty x^{-\alpha}e^{-\pi x}dx$$
$$I_1=\frac{e^\frac{i\pi(1-\alpha)}{2}}{(\pi n)^{1-\alpha}}\Gamma(1-\alpha)\tag{1}$$
$\displaystyle I_2$ we evaluate via integration by part
$$I_2=-\frac{1}{n^{1-\alpha}}\int_n^\infty x^{-\alpha}e^{i\pi x}dx=-\frac{1}{n^{1-\alpha}}\frac{e^{i\pi z}x^{-\alpha}}{i\pi}\bigg|_{x=n}^\infty-\frac{\alpha}{n^{1-\alpha}i\pi}\int_n^\infty x^{-\alpha-1}e^{i\pi x}dx$$
$$=i\frac{(-1)^{n+1}}{\pi n}+\frac{(-1)^{n+1}\alpha}{(\pi n)^2}+\frac{\alpha (1+\alpha)}{\pi^2 n^{1-\alpha}}\int_n^\infty x^{-\alpha-2}e^{i\pi x}dx\tag{2}$$
Taking together (1) and (2)
$$I=\frac{e^\frac{i\pi(1-\alpha)}{2}}{(\pi n)^{1-\alpha}}\Gamma(1-\alpha)+i\frac{(-1)^{n+1}}{\pi n}+\frac{(-1)^{n+1}\alpha}{(\pi n)^2}+\frac{\alpha (1+\alpha)}{\pi^2 n^{1-\alpha}}\int_n^\infty x^{-\alpha-2}e^{i\pi x}dx$$
Taking the real part
$$\int_0^1 t^{-\alpha} \cos (n \pi t)\, dt=\frac{\sin\frac{\pi\alpha}{2}}{(\pi n)^{1-\alpha}}\Gamma(1-\alpha)+\frac{(-1)^{n+1}\alpha}{(\pi n)^2}+O\Big(\frac{1}{n^4}\Big)\,,\,\, n\gg1$$
Integrating several times by part and using $\displaystyle \Gamma(\alpha)\alpha(1+\alpha)...(m+\alpha)=\Gamma(\alpha+m+1)$ we can get a complete asymptotic series:
$$\boxed{\,\,\int_0^1 t^{-\alpha} \cos (n \pi t)\, dt\sim\frac{\sin\frac{\pi\alpha}{2}}{(\pi n)^{1-\alpha}}\Gamma(1-\alpha)+\frac{(-1)^n}{\Gamma(\alpha)}\sum_{k=1}^\infty (-1)^k\frac{\Gamma(\alpha+2k-1)}{(\pi n)^{2k}}\,\,}$$
Best Answer
Let's start with defining the integral $I(x) = I_1(x) + iI_2(x)$ with
$$ I_1(x) = \int_0^1 \cos(g(t))\exp(2\pi ixt)dt, \quad I_2(x) = \int_0^1 \sin(g(t))\exp(2\pi ixt)dt, \quad x\to \infty $$
which then gives that the integral of interest is $I(x)$ with $g(t) = -2\pi\int_0^t \vartheta_3(-\pi t',\lambda)dt'$. Note that for each integral, there are no stationary points, this means that only the points $a=0$ and $b=1$ will contribute to the answer.
Given an integral on the form $I(x) = \int_{a}^b q(t)\exp(ixp(t))dt$ with no stationary points, it's value is asymptotically
$$ I(x) \sim \frac{ie^{ixp(a)}q(a)}{xp'(a)} - \frac{ie^{ixp(b)}q(b)}{xp'(b)}, \quad x\rightarrow \infty $$
So let's look at the theta function. We have
$$ \int_0^t \vartheta_3(-\pi t',\lambda)dt' = \int_0^t \left(1+2\sum_{n=1}^{\infty}\cos(2\pi nt') \right)dt = t+\sum_{n=1}^{\infty}\int_0^t \cos(2\pi nt')dt' = t +\sum_{n=1}^{\infty}\lambda^{n^2} \frac{\sin(2\pi nt)}{\pi n} $$
so
$$ g(t) = -2\pi t-2 \sum_{n=1}^{\infty}\lambda^{n^2} \frac{\sin(2\pi nt)}{n} $$
Clearly, $g(0)=0$ and $g(1) = -2\pi$. It is quite interesting to note that asymptotically, it is only the term $-2\pi t$ which contributes to the asymptotic expansion. I.e. if we plot the numerical solution where we truncate at different number of terms as a function of $x$, all will approach the same value for sufficiently large $x$.
For $I_1(x)$, let $p(t) = 2\pi t, q(t) = \cos(g(t))$. Note that $g(0) = 0 \Rightarrow q(0)=1$ and $g(1)=-2\pi\Rightarrow q(1) = 1$. Thus, one has
$$ I_1(x) \sim \frac{i}{2\pi x} \left(1 - e^{2\pi ix} \right) $$
with (one may check that this is true for both $x\to\pm\infty$)
$$ \Re(I_1(x)) \sim \frac{\sin(2\pi x)}{2\pi x}, \quad \Im(I_1(x)) \sim \frac{1}{2\pi x} (1-\cos(2\pi x)), \quad x\to \pm \infty $$
Now, we also look at $I_2(x)$, which is the same above except having $q(t)=\sin(g(t))$ instead. But at the endpoints, $q(0)=0$ and $q(1)=0$ in this case, so $I_2(x) \sim 0$ to leading order. Looking at this numerically, one sees that there are oscillations which decay faster than $O(x^{-1})$.
The asymptotics for $I(x)$ are seen below, compared to the numerical solution (here $\lambda = 0.5$).
Now, I'd like to add that I'm a bit unsure if I did things correctly. Clearly, the numerical integral I've plotted fits nicely, but I'm not certain if it's correctly implemented.o it would be nice to see if you can verify the results on your part as well.