Statistical Mechanics - Canonical Partition Function - An harmonic Oscillator 
*

*The problem statement, all variables and given/known data


With the Hamiltonian here:
Compute the cananonical ensemble partition function given by $\frac{1}{h} \int dq dp \exp^{-\beta(H(p,q)}$
for 1-d , where $h$ is planks constant


*Relevant equations

*The attempt at a solution
I am okay for the $p^2/2m$ term and the $aq^2$ term via a simple change of variables and using the gaussian integral result $\int e^{-x^2} dx = \sqrt{\pi}$
I am stuck on the $ \int dq e^{\beta   b q^{3}}$  and $ \int dq e^{\beta c q^{3}}$ terms.
If these were of the form $ \int dq e^{-\beta b q^{3}}$  I could evaluate via $\int dx e^{-x^n} = \frac{1}{n} \Gamma (1/n) $ where $ \Gamma(1/n) $ is the gamma function;
however because it is a plus sign I have no idea how to integrate forms of $ \int dq e^{x^n}$
Or should I be considering the integral over $q$ all together and there is another way to simply:
$\int dq e^{-\beta(aq^2-bq^3-cq^4)}$
PART B)
To compute the grand canonical ensemble of $N$ oscillators where
$ E_n=((n+\frac{1}{2})−x(n+\frac{1}{2})^2)\bar{h}ω $
to leading order in $x$
I have that
$\zeta(β,μ)=∑\lim^{N=∞}_{N=0}z^NZ(β,N)$
where $\zeta(β,μ)$ is the grand canonical ensemble, $Z(β,N)$ is the canonical ensemble and $z=e^{βμ} $ is the fugacity.
I have attempted to answer this question via summing first over $N$ in hope of getting a simplified expression for $Z_1$ to leading order in $x$ and then I will raise this to the power of $N$ , dividing by $N!$ (Gibbs factor) and then use (1) above to get the grand canonical partition function.
My attempt is as follows:
$Z_1 = \sum_n e^{-\beta((n+\frac{1}{2})\bar{h}\omega-x(n+\frac{1}{2})^2\bar{h}\omega)}=e^{-\frac{\beta \bar{h} \omega}{2}}\sum_n e^{-\beta \bar{h} \omega n} e^{\beta \bar{h} \omega x (n+\frac{1}{2})^2}$
I will expand out the exponential $x$ term to get
$ e^{-\frac{\beta \bar{h} \omega}{2}} \sum_n e^{-\beta \bar{h} \omega n} (1+\beta\bar{h}\omega x(n+\frac{1}{2})^2) + O(x^2)$
$e^{-\frac{\beta \bar{h} \omega}{2}}\sum_n e^{-\beta \bar{h} \omega n} + e^{-\frac{\beta \bar{h} \omega}{2}}e^{-\beta \bar{h} \omega n} (\beta \bar{h} \omega x) \sum_{n=0}^{n=\infty} (n^2+ \frac{1}{4} + n)$
which diverges...
Many thanks in advance

This is my first answer, so I hope I'm doing it right.
As pointed out in an earlier comment, I think you need to start of by getting the limits straight, which will answer a couple of your questions. The integral over $p$ is independent and easily done as you've stated yourself. The integral over $q$ goes from $-\infty$ to $+\infty$, as it is the position in one dimension. 
Note in passing that it is $$\int_0^{\infty} e^{-x^n} = \frac{1}{n}\Gamma\left(\frac{1}{n}\right)$$ but your lower limit is $-\infty$, and so this cannot be used.
(Incidentally, $\int_{-\infty}^{+\infty}e^{\pm x^3}dx$ does not converge to the best of my knowledge).
But all of this is beside the point: unless I've misunderstood you (please correct me if I'm wrong!), you're claiming that 
$$\int_{-\infty}^{+\infty}dq \,\,e^{-\beta a q^2 + \beta b q^3 + \beta c q^4} = \int_{-\infty}^{+\infty}dq \,\,e^{-\beta a q^2} \int_{-\infty}^{+\infty}dq \,\,e^{\beta b q^3} \int_{-\infty}^{+\infty}dq \,\,e^{ \beta c q^4}$$
which is clearly not true. So performing the integrals separately is not the way to go and you must consider the integral over all the functions of $q$ together. If the extra terms had been linear in $q$ you could have used the "completing the square" trick, but I don't there is anything similar for higher powers. For an exact (and possibly useless) answer, you might be able to use this formula that is completely beyond me, but I don't think it's helpful.
With a little research I found the original question online. If it's the same question, it says 

You should work in the approximation where the anharmonicity is small

In other words, when 

$$\frac{b}{a} \ll 1    \quad \quad \quad \frac{c}{a^2}\ll 1$$

Which leads me to believe that they require you to perform some sort of series in the two "problematic" terms. 
If I were to solve this question, I would do the following. I'd begin by making a simple variable substitution: $u = \sqrt{a} q$, which makes my integral
$$\frac{1}{\sqrt{a}} \int_{-\infty}^\infty du\,\,\exp{\left(-\beta u^2 + \frac{\beta b}{a\sqrt{a}}u^3 + \frac{\beta c}{a^2} u^4\right)} = \frac{1}{\sqrt{a}} \int_{-\infty}^\infty du\,\, e^{-\beta u^2}e^{\frac{\beta b}{a\sqrt{a}}u^3} e^{\frac{\beta c}{a^2} u^4}$$
I would then approximate the anharmonic terms using $e^{ax}\approx 1 + ax + \frac{a^2 x^2}{2}...$
$$I = \frac{1}{\sqrt{a}} \int_{-\infty}^\infty du\,\, e^{-\beta u^2} \left(1 + \frac{\beta b}{a\sqrt{a}}u^3 + \frac{\beta^2 b^2}{2 a^3}u^6\right) \left(1+ \frac{\beta c}{a^2} u^4\right)$$
You'll notice I took three terms for the first approximation and only two for the second. Why I did this will become clear in a moment. This is now just a sum of integrals of the form
$$\int_{-\infty}^\infty du\,\, e^{-\beta u^2}u^n$$
Clearly, when $n$ is an odd integer, these integrals are zero, since the function is odd and its positive and negative contributions cancel each other. This is why when expanding the series earlier I stopped at the first even power, to get the leading dependencies in $b$ and $c$.
We are thus left with three integrals (the fourth term involving $\frac{b^2 c}{a^5}$ can be neglected since it is of a much higher order.)
$$I = \frac{1}{\sqrt{a}}\left( \int_{-\infty}^\infty du\,\, e^{-\beta u^2}  +  \frac{\beta c}{a^2}  \int_{-\infty}^\infty du\,\, e^{-\beta u^2}u^4 +  \frac{\beta^2 b^2}{2 a^3}\int_{-\infty}^\infty du\,\, e^{-\beta u^2}u^6\right)$$
The remaining (even) integrals can be easily computed realising that whenever $2z-1$ is even, 
$$\int_{-\infty}^{\infty} u^{2z-1} e^{-\beta u^2}du = \beta^{-z} \Gamma[z]$$
Thus,
$$I \approx \sqrt{\frac{\pi}{\beta a}}\left( 1 + \frac{3 c}{4 a^2 \beta} + \frac{15 b^2}{16 a^3 \beta}\right)$$

EDIT:
Upon some introspection I've decided that it does indeed make sense to have a term in $\frac{b^2}{a^3}$. I've made the necessary changes.
Furthermore, the complete canonical partition function (including the integral over $p$) is thus
$$Z = \frac{\pi}{\beta}\sqrt{\frac{2 m}{a}}\left( 1 + \frac{3 c}{4 a^2 \beta} + \frac{15 b^2}{16 a^3 \beta}\right)$$
