In the post, I tackled the integral by power series and integration by parts and obtained that
$$
\int_0^{\frac{\pi}{2}} \frac{x^2}{\sin x} d x=2\pi G-\frac{7}{2}\zeta(3)
$$
where $G$ is the Catalan’s constant.
Similarly, we can express the integrand as a power series of exponential functions.
$$
\begin{aligned}
I_n &=\int_0^{\frac{\pi}{2}} \frac{x^n}{\sin x} d x\\&=2 i \int_0^{\frac{\pi}{2}} \frac{x^n}{e^{x i}-e^{-x i}} d x \\
& =2 i \int_0^{\frac{\pi}{2}} x^{n-1} e^{-x i} \sum_{k=0}^{\infty} e^{-2 k x i} d x \\
& =2 i \sum_{k=0}^{\infty} \int_0^{\frac{\pi}{2}} x^{n-1} e^{-(2 k+1) x i} d x \\
& =2 i \sum_{k=0}^{\infty}\left[\int_0^{\frac{\pi}{2}} x^n \cos (2 k+1) x d x-i \int_0^{\frac{\pi}{2}} x^n \sin (2 k+1) x\right.
\end{aligned}
$$
Comparing their real parts yields
$$
I_n=2 \sum_{k=0}^{\infty} \underbrace{ \int_0^{\frac{\pi}{2}} x^n \sin (2 k+1) x d x}_{J_n(k)}
$$
Using integration by parts twice, we get a reduction formula for $J_n(k)$.
$$
\begin{aligned}
J_n(k) & =-\frac{1}{2(2+1}\left[x^n \cos (2 k+1) x\right]_0^{\frac{\pi}{2}}+\frac{n}{2 k+1} \int_0^{\frac{\pi}{2}} x^{n-1} \cos (2 k+1) x d x \\
& =\frac{n}{(2 k+1)^2}\int_0^{\frac{\pi}{2}} x^{n-1} d(\sin (2 k+1) x) \\
& =\frac{n}{(2 k+1)^2}\left[x^{n-1} \sin (2 k+1) x\right]_0^{\frac{\pi}{2}}-\frac{n(n-1)}{(2 k+1)^2 }\int_0^{\frac{\pi}{2}} x^{n-2} \sin (2 k+1) xdx \\
& =\frac{(-1)^k n}{(2 k+1)^2}\left(\frac{\pi}{2}\right)^{n-1}-\frac{n(n-1)}{(2 k+1)^2} J_{n-2}(k)
\end{aligned}
$$
$$J_n(k) =\frac{(-1)^k n}{(2 k+1)^2}\left(\frac{\pi}{2}\right)^{n-1}-\frac{n(n-1)}{(2 k+1)^2} J_{n-2}(k) \tag*{(1)} $$
Plugging back summation into the formula yields
$$
\boxed{ I_n=\frac{n \pi^{n-1}}{2^{n-2}} G-2 n(n-1) \sum_{k=0}^{\infty} \frac{J_{n-2}(k)}{(2 k+1)^2}} \tag*{(2)}
$$
Using the formula $(2)$, we can get
$$
\begin{aligned} \int_0^{\frac{\pi}{2}} \frac{x^2}{\sin x} d x& =2 \pi G-4 \sum_{k=0}^{\infty} \frac{J_0(k)}{(2 k+1)^2} \\
& =2 \pi G-4 \sum_{k=0}^{\infty} \frac{1}{(2 k+1)^3} \\
& =2 \pi G-\frac{7}{2} \zeta(3)
\end{aligned}
$$
$$
\begin{aligned}
\int_0^{\frac{\pi}{2}} \frac{x^3}{\sin x} d x & =\frac{3 \pi^2}{2} G-12 \sum_{k=0}^{\infty} \frac{J_1(k)}{(2 k+1)^2} \\
& =\frac{3 \pi^2}{2} G-12 \sum_{k=0}^{\infty} \frac{(-1)^k}{(2 k+1)^4} \\
& =\frac{3 \pi}{2} G-12\beta(4)
\end{aligned}\\
$$
where $\beta(*)$ is the Dirichlet beta function.
$$
\int_0^{\frac{\pi}{2}} \frac{x^4}{\sin x} d x=\pi^3 G-24 \sum_{k=0}^{\infty} \frac{J_2(k)}{(2 k+1)^2}
$$
Using the reduction formula $(1) $ for $J_n(k)$, we have
$$
\begin{aligned}
J_2(k) & =\frac{2(-1)^k}{(2 k+1)^2}\left(\frac{\pi}{2}\right)-\frac{2}{(2 k+1)^2} J_0 \\
& =\frac{\pi(-1)^k}{(2 k+1)^2}-\frac{2}{(2 k+1)^3}
\end{aligned}
$$
Plugging back yields
$$
\begin{aligned}
\int_0^{\frac{\pi}{2}} \frac{x^4}{\sin x} d x=&\pi^3 G-24 \pi \sum_{k=0}^{\infty} \frac{(-1)^k}{(2 k+1)^4}+48 \sum_{k=0}^{\infty} \frac{1}{(2 k+1)^5} \\
= & \pi^3 G-24 \pi\beta(4) +\frac {93}2\zeta(5) \\
\end{aligned}
$$
checked by WA.
Theoretically, we can walk with $I_n$ as far as we like using both formula $(1)$ and $(2)$ though the work is tedious and long.
Can we go further with $n\ge m\ge 2$?
Comments and methods are highly appreciated.
Best Answer
Consider integration on a rectangular contour $C$ $$ f(z)=\frac{z^n}{\sin z}\\ C:0\to \frac\pi2 \to \frac\pi2 +i\infty\to i\infty \to0 $$ sorry for not presenting the picture, but I think it is easy to see
Note that $\sin z$ only has real roots, so $f(z)$ has no poles in the contour, then $$ \oint f(z)dz=\lim_{R\to\infty}\int_R^0 f(iy)i~dy+\int_0^{\pi/2}f(x)dx+\int_0^R f\left(\frac\pi2+iy\right)i ~dy+\int_{\pi/2}^0 f(x+ iR)dx=0 $$ Easily $|f(x+i R)|=O(R^n e^{-R})\to0$, so the last term vanishes. Hence $$ \int_0^{\pi/2}f(x)dx=\int_0^\infty f\left(iy\right)i ~dy-\int_0^\infty f\left(\frac\pi2+iy\right)i ~dy =\int_0^\infty \frac{(iy)^n}{\sinh y}-\frac{i\left(\pi/2+iy\right)^n}{\cosh y} ~dy $$ Finally, ultilize the classic result $$ \int_0^\infty \frac{y^n}{\sinh y}dy=2=(2-2^{-n})n!\zeta(n+1)\\ \int_0^\infty \frac{y^n}{\cosh y}dy=2n!\beta(n+1) $$ we arrive at $$ I_n=\int_0^{\pi/2}\frac{x^n}{\sin x}dx=i^n(2-2^{-n})n!\zeta(n+1)-2\sum_{k=0}^n\frac{n!}{(n-k)!}\left(\frac\pi2\right)^{n-k}i^{k+1}\beta(k+1) $$ Since the integral is real, taking the real and imaginary part of RHS provides the value of the integral and an identity (seemingly useless, though). Some explicit results: $$ \begin{align} &I_1= 2 G \\ &I_2= 2 \pi G-\frac{7 \zeta (3)}{2} \\ &I_3= \frac{3 \pi ^2 G}{2}-12 \beta(4) \\ &I_4= -24 \pi \beta(4)+\pi ^3 G+\frac{93 \zeta (5)}{2} \\ &I_5= -30 \pi ^2 \beta(4)+240 \beta(6)+\frac{5 \pi ^4 G}{8} \\ \end{align} $$ With slight modification, the method generalizes to all sorts of integrals of the form $\text{Rational of } e^{ix}\times \text{polynomial}$ on any intervals, only involving special functions up to polylogarithms.
For example $$ \int_0^{\pi/2} \frac{x^3}{(\sin x+1) (\cos x+1)^2} dx=\\ \frac{5 \pi ^3}{24}+\frac{\pi ^2}{2}+\pi -\left(12+14 \pi +3 \pi ^2\right)G+\left(-2+3 \pi +\frac{13 \pi ^2}{4}-\frac{\pi ^3}{4}\right) \log (2)+\left(\frac{21}{2}+\frac{63 \pi }{8}\right) \zeta (3) $$
$$ \int _0^{\pi/4}\frac{x^4}{\sin^4x}dx=\frac{\pi ^2 G}{4}+2G-2 \beta (4)-\frac{3 \pi \zeta (3)}{32}-\frac{\pi ^2}{8}-\frac{\pi ^3}{48}-\frac{\pi ^4}{192}+\frac{1}{48} \pi ^3 \log (2)+\frac{1}{2} \pi \log (2) $$
$$ \int_0^\infty\frac{\arctan^7x}{x^6}dx=\left(-\frac{147}{8}+\frac{315 \pi ^2}{16}-\frac{63 \pi ^4}{128}\right) \zeta (3)\\ +\left(\frac{945 \pi ^2}{128}-\frac{3255}{32}\right) \zeta (5)-\frac{8001 \zeta (7)}{256}+\frac{21 \pi ^6}{1280}-\frac{7 \pi ^4}{64}+\left(\frac{21 \pi ^2}{4}-\frac{35 \pi ^4}{16}+\frac{7 \pi ^6}{320}\right) \log (2) $$
Update
OP requested a more general form. The close form does exist.
Let $$ S(m,n)=\int_0^{\pi/2}\frac{x^n}{\sin^mx}\text d x $$ Same contour integration yields $$ S(m,n) =\int_0^\infty \frac{(iy)^n}{\sinh^m y}-\frac{i\left(\pi/2+iy\right)^n}{\cosh^m y} ~dy $$ Now it suffice to evaluate the Mellin transform $$ \int_0^\infty \frac{x^{s-1}}{\sinh^mx}dx\qquad \int_0^\infty \frac{x^{s-1}}{\cosh^mx}dx $$ The basic idea goes like this:
For all rational function with the only pole $t=1$, namely $\dfrac{P(t)}{(1-t)^q}$ for some polynomial $P$ such that $P(1)\ne0$ and some positive integer $q$, one can always expand it into a linear combination of $\text{Li}_{-q}, q\in\mathbb N$. Note that $$ \text{Li}_{-q}=\left( t\frac{d}{dt}\right)^q\frac1{1-t}=\sum_{k\ge0}k^q~t^k $$ is always a rational function.
Then combine partial fraction decomposition and utilize (assume absolute convergence) $$ \int_0^\infty x^{s-1}\text{Li}_r(ze^{-x}) dx=\sum_{k\ge0}\frac{z^k}{k^r}\int_0^\infty x^{s-1}e^{-lx} dx=\Gamma(s)\sum_{k\ge0}\frac{z^k}{k^{r+s}}=\Gamma(s)\text{Li}_{r+s}(z) $$ showing that the Mellin transform of any rational function of $e^{-x}$ is expressible in a linear combination of polylogarithms multiplied by $\Gamma(s)$ and, in this case, it can be converted into $\zeta$ and $\beta$ s.
Here I would simply state that $$ \frac1{\sinh^mx}=\sum_{l\ge0}\left(\sum_{k=0}^m(1+(-1)^{l+m})c_k^{(m)}l^k\right)e^{-lx}\\ \frac{(2t)^m}{(1-t^2)^m}=\sum_{k=0}^mc_k^{(m)}\Big(\text{Li}_{-k}(t)+(-1)^m\text{Li}_{-k}(-t)\Big) $$ where $t=e^{-x}$ and $\{c_k^{(m)}\}$ are rational coefficients whose explicit expression contains triple summation involving Stirling numbers.
I am tired of it and do not want to write down all that stuff here. Perhaps some one out there are willing to do calculation for whom I will leave it.
Hence, $$ \int_0^\infty \frac{x^{s-1}}{\sinh^mx}dx=\Gamma(s)\sum_{k=0}^mc_k^{(m)}\Big(1+(-1)^m(1-2^{s-k})\Big)\zeta(s-k) $$ As for the other one, let $t\to it$ to get $$ \frac{(2t)^m}{(1+t^2)^m}=i^{-m}\sum_{k=0}^mc_k^{(m)}\Big(\text{Li}_{-k}(it)+(-1)^m\text{Li}_{-k}(-it)\Big) $$
$$ \int_0^\infty \frac{x^{s-1}}{\cosh^mx}dx=i^{-m}\Gamma(s)\sum_{k=0}^mc_k^{(m)}\Big(\text{Li}_{s-k}(i)+(-1)^m\text{Li}_{s-k}(-i)\Big) $$
Utilizing $$ \text{Li}_s(i)=(2^{1-2s}-2^{-s})\zeta(s)+i \beta(s) $$ one can verify the result reduces to $\zeta$ s when $m$ is even, and $\beta$ s when $m$ is odd.
Finally, put everything together. $$ S(m,n)=n!\sum_{k=0}^{n+1}c_{n-k+1}^{(m)}\Big(1+(-1)^m(1-2^{k})\Big)\zeta(k)\\ -\sum_{r=0}^n\frac{n!}{(n-r)!}\left(\frac\pi2\right)^{n-r}i^{r+1-m}\sum_{k=0}^mc_k^{(m)}\Big(\text{Li}_{r+1-k}(i)+(-1)^m\text{Li}_{r+1-k}(-i)\Big) $$ Writing the $c^{(m)}_k$ s explicitly results in a monstrous quintuple summation. Indeed the explicit result anyway.
Frankly I would use CAS to find those coefficients rather than do the tedious summation, so I wrote a Mathematica code to evaluate any integral of this form.
Here's a list for $0<m\le n\le5$ $$ \int _0^{\pi/2}x \csc (x)dx=2 G\\ \int_0^{\frac{\pi }{2}}x^2 \csc (x)dx=2 \pi G-\frac{7 \zeta (3)}{2}\\ \int _0^{\pi/2}x^2 \csc ^2(x)dx=\pi \log (2)\\ \int _0^{\pi/2}x^3 \csc (x)dx=\frac{3 \pi ^2 G}{2}-12 \beta(4)\\ \int _0^{\pi/2}x^3 \csc ^2(x)dx=\frac{3}{4} \pi ^2 \log (2)-\frac{21 \zeta (3)}{8}\\ \int _0^{\pi/2}x^3 \csc ^3(x)dx=-6 \beta(4)+\frac{3 \pi ^2 G}{4}+6 G-\frac{3 \pi ^2}{8}\\ \int _0^{\pi/2}x^4 \csc (x)dx=-24 \pi \beta(4)+\pi ^3 G+\frac{93 \zeta (5)}{2}\\ \int _0^{\pi/2}x^4 \csc ^2(x)dx=\frac{1}{2} \pi ^3 \log (2)-\frac{9 \pi \zeta (3)}{4}\\ \int _0^{\pi/2}x^4 \csc ^3(x)dx=-12 \pi \beta(4)+\frac{\pi ^3 G}{2}+12 \pi G-21 \zeta (3)+\frac{93 \zeta (5)}{4}-\frac{\pi ^3}{4}\\ \int _0^{\pi/2}x^4 \csc ^4(x)dx=-\frac{3 \pi \zeta (3)}{2}-\frac{\pi ^3}{12}+\frac{1}{3} \pi ^3 \log (2)+2 \pi \log (2)\\ \int _0^{\pi/2}x^5 \csc (x)dx=-30 \pi ^2 \beta(4)+240 \beta(6)+\frac{5 \pi ^4 G}{8}\\ \int _0^{\pi/2}x^5 \csc ^2(x)dx=-\frac{45 \pi ^2 \zeta (3)}{16}+\frac{465 \zeta (5)}{32}+\frac{5}{16} \pi ^4 \log (2)\\ \int _0^{\pi/2}x^5 \csc ^3(x)dx=-15 \pi ^2 \beta(4)-120 \beta(4)+120 \beta(6)+\frac{5 \pi ^4 G}{16}+15 \pi ^2 G-\frac{5 \pi ^4}{32}\\ \int _0^{\pi/2}x^5 \csc ^4(x)dx=-\frac{15 \pi ^2 \zeta (3)}{8}+\frac{155 \zeta (5)}{16}-\frac{35 \zeta (3)}{4}-\frac{5 \pi ^4}{96}+\frac{5}{24} \pi ^4 \log (2)+\frac{5}{2} \pi ^2 \log (2)\\ \int _0^{\pi/2}x^5 \csc ^5(x)dx=-\frac{45}{4} \pi ^2 \beta(4)-100 \beta(4)+90 \beta(6)+\frac{15 \pi ^4 G}{64}+\frac{25 \pi ^2 G}{2}+10 G-\frac{5 \pi ^2}{8}-\frac{55 \pi ^4}{384} $$