Taylor series and approximations to complicated integrals So, I feel like I should technically know how to do this, but I'm really not sure. I have a certain integral to calculate or rather, approximate since it's quite hard (if you'd like to know, it's the one in this question and it's been quite painful. Now, I've managed, with a considerable amount of manipulation, to reduce it to
$$I(k) \approx \int_0^\infty \frac{p^6}{(1 + a^2 p^2)(1 + b^4 p^4)} \frac{1}{(1 + a^2 (p^2+k^2-p k))(1 + b^4 (p^2+k^2-p k)^2)}dp$$
which was still too hard. I then decided to make a further approximation, writing it 
$$I(k) \approx \int_0^\infty \frac{p^6}{(1 + a^2 p^2)(1 + b^4 p^4)} \frac{1}{(1 + a^2 (p^2+k^2))(1 + b^4 (p^2+k^2)^2)}dp$$
which is slightly (but only ever so slightly) more doable. The problem is, the end result is an extremely complicated result and I'd like to approximate it even further. (The reason I actually had to calculate the integral was I didn't know how it depended on the parameters $a$ and $b$).
Now, I've got the solution and its dependence on the parameters, but it's extraordinarily ugly (it's about two pages long, so I won't include it here!). Furthermore, this is just the start of a bunch of complicated calculations I'll have to perform, and so breaking it up into some manageable chunks of some form would be very helpful.  I'm mainly interested in the "small" and "large" regimes, as well as where the function "changes" the overall power law. For example, in the following Log-Log plot, you see there is a "knee" around the region of "5" to "10", at different points for different parameters of the problem. 

Questions:
1) Does anyone have any ideas on how to break this up into a more manageable function say, using a Taylor Series? An earlier study showed that it sort of went as $$\sim \frac{1}{1 + c_1 k^2 + c_2 k^5}$$ but in that case the integral was enormously simpler (it still needed to be approximated by manipulating $c_1$ and $c_2$ though) and there were no parameters like $a$ and $b$ involved.
I can get rid of one of the parameters by rescaling the variable but I'm still left with another. I can then do a Taylor series approximation around $0$ to get the leading order terms, and another around $\infty$ to get the asymptotic dependence. However, in the intermediate region I'm basically working in the dark. 
2) Can the Taylor series around some pivotal point could help me in approximating the function in this region. I can easily perform such a series, but it gives me results in powers of $(k-k_0)$, which work locally, but not globally (i.e. I couldn't just "attach" it to the other initial and asymptotic values.
3) Furthermore, if I "knew" there should be (say) a $k^{-1}$ behaviour around a certain point (say $k_0$), would I be able to verify this using a Taylor series? Would I be able to get it's coefficient somehow?
I'm sorry if these seem like evident questions, I really feel like I ought to know the answers (or atleast if it's possible!), given how much I've used these series in Physics! I'd really appreciate any help, I've wasted an insane amount of time on this, but I can't move forward without even a roughly analytic form of the graph.
Let me know if there's anything I haven't explained clearly!

EDIT:
Couple of points to add on the values of the parameters $a,b$ and the variable $k$. $a$ and $b$ are small, $b\sim 1/3$ and $|a|<1$ for physically acceptable situations. However, $k$ can in theory stretch from $0$ to $\infty$ (though we usually stop around $O(1000)$.
 A: If $a$ and $b$ are neither small nor large, the best you can do in terms of approximation is to consider $k\gg 1$ and expand around $k=\infty$.
In this regime the important region of the integral is when $p$ is of the order of $k$, so write $p=xk$. Then, the ''$1+$'' part in the denominators can be neglected and you get (from the first form of the problem)
$$\frac{1}{a^4b^8k^5}\int_0^\infty \frac{1}{(1-x+x^2)^3}dx$$
Corrections can be found via geometric series for the denominators
A: I noticed that in the question that you previously posed, you need a closed form solution because you intend to summon it (presumably with arguments $a$ and $b$) for use elsewhere.
Assuming that my understanding is correct, I would suggest creating a response surface of the values of the integral $U(a,b)$ in 2 dimensions based on numerically evaluating the integral over the entire spectrum of $a$ and $b$.
As far as an analytical approach goes, I don't have any suggestions other than the excellent ones already put forth (along with their drawbacks). I would caution that any attempt using a Taylor's series would entail the computation of the derivatives at some reference value ($p=0$ seems logical). Since your integral extends out to $\infty$, you'll need a LOT of them! While you may be able to use Mathematica to compute these derivatives, you'll very likely encounter an even more hideous mess than you have so far already! And the truncation errors may be unacceptably large in the end rendering the effort futile.
All the best.
