Too long for a comment, so I thought might as well write it in detail.
Instead of the Laplace transform, which seems complicated to calculate, if we use the Mellin transform as Mariusz has done, the integral becomes more manageable.
The Mellin transform of a function $\ f(x)$ with respect to $\ x$ is defined as:
$$\mathcal M(f(x))(s) = \int_0^{\infty} x^{s-1} \cdot f(x) \cdot dx$$
The Mellin transform also has some nice properties that allows us to extract an expression for $\ I(a,b)$, though this makes me wonder - shouldn't there be some relationship between the Laplace and Mellin transforms?
I have only recently began learning about this transform, so I don't have any concrete ideas yet, but their relationship is something worth investigating later.
Anyways, here is a solution via Mellin transform. This still doesn't resolve the original doubt of proceeding via Laplace transform, but perhaps this one gives an insight to a solution using the latter.
If we take the Mellin transform of $\ I(a,b)$ with respect to $\ a$, we get:
$$\mathcal M(I)(s) = \int_0^{\infty} a^{n-1}\cdot I(a,b) \cdot da$$
$$= \int_{a=0}^{\infty} a^{n-1}\cdot \int_{x=0}^{\infty} \frac{\cos(ax)}{x^b+1} dx \cdot da$$
Swapping the order of integration, we get:
$$\mathcal M(I)(s) = \int_{0}^{\infty} \frac{1}{x^b+1} \cdot \mathcal M({\cos(ax)}) dx$$
Side note: You can derive that the Mellin transform of cos(nx) with respect to x is
$$\mathcal M(cos(nx))(s) = \frac{1}{n^s}\cdot \Gamma(s) \cdot \cos(\frac{s\pi}{2})$$
by exploiting the fact that
$$\mathcal M(e^{-px})(s) = \frac{1}{p^s}\cdot \Gamma(s)$$
Which itself follows straight from a simple substitution and the definition of the Gamma function.
Then we simply take the real part of $\mathcal M(e^{-inx})$ [i.e. $\ p = in$], using Euler's formula to arrive at the result.
Continuing forth, we get:
$$\mathcal M(I)(s) = \int_0^{\infty} \frac{1}{x^b+1} \cdot \frac{1}{x^s}\cdot \Gamma(s) \cdot \cos(\frac{s\pi}{2}) dx$$
$$ = \Gamma(s) \cdot \cos(\frac{s\pi}{2}) \cdot \int_0^{\infty} \frac{1}{x^b+1} \cdot \frac{1}{x^s} dx$$
Now, we arrive at the relatively hard part. We shall substitute
$$\ t = \frac{1}{x^b+1}$$
This implies that
$$\ x = (\frac{1-t}{t})^{\frac{1}{b}}$$
and
$$\ dx = \frac{-1}{b} \cdot (1-t)^{\frac{1}{b}-1} \cdot t^{\frac{-1}{b}-1}dt$$
After changing the upper and lower bounds to suitable values and swapping the limits of integration to get rid of the negative sign, we have:
$$\mathcal M(I)(s) = \Gamma(s) \cdot \cos(\frac{s\pi}{2}) \cdot \int_0^1 (\frac{1-t}{t})^{\frac{-s}{b}} \cdot t \cdot \frac{1}{b} \cdot (1-t)^{\frac{1}{b}-1} \cdot t^{\frac{-1}{b}-1}dt$$
$$\mathcal M(I)(s) = \Gamma(s) \cdot \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \int_0^1 t^{(\frac{s-1}{b})}\cdot (1-t)^{(\frac{1-s}{b}-1)} dt$$
$$\mathcal M(I)(s) = \Gamma(s)\cdot \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \beta(\frac{1-s}{b}, \frac{s-1}{b}+1)$$
Where $\beta(x,y)$ is the Beta function.
Using the identities of the Beta function, we further obtain:
$$\mathcal M(I)(s) = \Gamma(s)\cdot \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \Gamma(\frac{1-s}{b}) \cdot \Gamma(\frac{s-1}{b}+1)$$
$$ = \Gamma(s)\cdot \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \frac{\pi}{\sin[\pi\cdot (\frac{1-s}{b})]}$$
By Ramanujan's Master theorem, if a function $\ f(x)$ can be written as:
$$\ f(x) = \sum_{k=0}^{\infty} \frac{\phi(k)}{k!} \cdot (-x)^k$$
Then
$$\mathcal M(f(x))(s) = \Gamma(s) \cdot \phi(-s)$$
Note that in our case, we have
$$\ M(I)(s) = \Gamma(s) \cdot \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \frac{\pi}{\sin[\pi\cdot (\frac{1-s}{b})]}$$
By comparison,
$$\phi(-s) = \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \frac{\pi}{\sin[\pi\cdot (\frac{1-s}{b})]}$$
This implies
$$\phi(s) = \frac{1}{b} \cdot \cos(\frac{s\pi}{2}) \cdot \frac{\pi}{\sin[\pi\cdot (\frac{s+1}{b})]}$$
This means that our original integral $\ I(a,b)$ can be written as:
$$\ I(a,b) = \sum_{k=0}^{\infty} \frac{\phi(k)}{k!} \cdot (-a)^k$$
Which finally gives us:
$$\ I(a,b) = \frac{\pi}{b} \cdot \sum_{k=0}^{\infty} \frac{(-a)^k}{k!} \cdot \frac{\cos(k\frac{\pi}{2})}{\sin(\frac{\pi}{b}\cdot (1+k))}$$
Note that for $\ b = 2$, a special case, the $\cos(k\frac{\pi}{2})$ and $\sin(\frac{\pi}{b}\cdot (1+k))$ terms cancel each other, since they are equal by their properties, and we are left with:
$$\ I(a,2) = \frac{\pi}{2} \cdot \sum_{k=0}^{\infty} \frac{(-a)^k}{k!}$$
Which is the Maclaurin series of $\ e^x$ evaluated at $\ x= -a$, thus giving us the well known result:
$$\ I(a,2) = \int_0^{\infty} \frac{cos(ax)}{x^2+1} dx = \frac{\pi}{2}\cdot e^{-a}$$
Similarly, for $\ b =4$, we have:
$$\ I(a,4) = \frac{\pi}{4} \cdot \sum_{k=0}^{\infty} \frac{(-a)^k}{k!} \cdot \frac{\cos(k\frac{\pi}{2})}{\sin(\frac{\pi}{4}\cdot (1+k))}$$
A few things to observe:
- For $\ k = 4n + 1$, the term $\frac{\cos(k\frac{\pi}{2})}{\sin(\frac{\pi}{4}\cdot (1+k))}$ becomes 0.
- For $\ k = 4n + 2$, the term $\frac{\cos(k\frac{\pi}{2})}{\sin(\frac{\pi}{4}\cdot (1+k))}$ becomes -$ \sqrt2$.
- For $\ k = 4n + 3$, the term $\frac{\cos(k\frac{\pi}{2})}{\sin(\frac{\pi}{4}\cdot (1+k))}$ is indeterminate.
We take the limit as $\ k$ approaches $3$ and find it to be -$2$, while the limit as it approaches $7$ is $2$ and it oscillates between these two values for greater $n$.
- For $\ k = 4n + 4$ the term $\frac{\cos(k\frac{\pi}{2})}{\sin(\frac{\pi}{4}\cdot (1+k))}$ becomes $ \sqrt2$.
Still confused how to solve for a general $\ b$ though, apart from computing that terrible sum.
Best Answer
Making change of variable $s=t^2$, we have for $x\ne 1$, $$\begin{align*} I&=\int_0^{\infty} \frac{t}{(t^2+1)(t^2+x^2)} \mathrm{d}t\\&=\frac12\int_0^\infty \frac1{(s+1)(s+x^2)}\mathrm ds\\&=\frac1{2(x^2-1)}\int_0^\infty\left( \frac1{s+1}-\frac1{s+x^2}\right)\mathrm ds\\&=\frac1{2(x^2-1)}\left[\ln\left(\frac{s+1}{s+x^2}\right)\right]^\infty_0\\&=\frac{\ln x}{x^2-1}. \end{align*}$$ For $x=1$, we have $I =\lim_{x\to 1}\frac{\ln x}{x^2-1}=\frac12$. We can also do it explicitly: $$ \frac12\int_0^\infty \frac1{(s+1)^2}\mathrm ds=\left[-\frac1{2(s+1)}\right]^\infty_0=\frac12. $$