This is inspired by the same beautiful integral expression for $\zeta(3)$ as this question, but goes in a slightly different direction. Writing the original integral in the form $$\int_0^1\frac{x(1-x)}{\sin\pi x}dx=7\frac{\zeta(3)}{\pi^3} ,$$ it turns out that for $n\in\mathbb N$ there is a unique monic polynomial $p_n$ of degree $n-1$ such that $$\int_0^1x^np_n(x)\frac{1-x}{\sin\pi x}dx=c_{2n+1}\frac{\zeta(2n+1)}{\pi^{2n+1} } $$ with rational $c_k=4(k-1)! \dfrac{2^k-1}{2^k}= (4-2^{2-k})(k-1)! $.
This follows for $\zeta(2n+1)$ from solving the linear system given by the blue lines numbered $n+1,…,2n$ in the other question. As Zurab Silagadze answered it by giving an explicit formula of the coefficients in the blue lines, the $p_n$'s can be calculated. I don't know however if it is possible to give a formula in closed form (meaning here that it should not contain a matrix inversion), but see below.
The first polynomials are $$\begin{align}
p_1(x)&=1 \\
p_2(x)&=3-x \\
p_3(x)&=25-20x+x^2 \\
p_4(x)&=455-707x+287x^2-x^3 \\
p_5(x)&=14301-34734x+29046x^2-8304x^3+x^4 \\
p_6(x)&=683067-2289309x+2949276x^2-1721434x^3+382547x^4-x^5 \\
\end{align}$$
The constant terms are supposedly the sequence A272482, thus, correcting the oeis typo $1/(2n)!$, $$[x^0]p_n(x)= {(2n)!}[x^{2n}y^n]\frac{\cos\frac{x(1-y)}{2}} {\cos\frac{x(1+y)}{2}} = \frac 1{4^n} {2n\choose n}\sum_{i=0}^n{n\choose i}E_i,$$ where $E_i$ are the Euler numbers. This seems to suggest something similar for the other coefficients, and thus possibly a closed form.
Are the $p_n$ known? How to find their closed form or generating function?
More generally now, define $$J(m,n,k)=J(n,m,k):= \int_0^1\frac{x^m(1-x)^n}{\sin^k\pi x}dx.$$For this to converge, we need $m,n\geqslant k$.
Experimentally, the situation for $k=2$ is quite similar to the $k=1$ case in that $J(m,n,k)$ is a rational combination of values $\dfrac{\zeta(i)}{\pi^{i+1}}$ with $i$ running over all odd numbers between $\min(m,n)$ and $m+n-1$, e.g. $$J(7,4,2)=\dfrac{105}2\left(-\dfrac{\zeta(5)}{\pi^6}+51\dfrac{\zeta(7)}{\pi^8}-405\dfrac{\zeta(9)}{\pi^{10}}\right).$$
This leads to new possibilities of representing odd zeta values as integrals, this time with $\sin^2\pi x$ in the denominator. Writing as a shortcut $h_m:=J(m,m,2)$, we can for example express $\dfrac{\zeta(2n-1)}{\pi^{2n}}$ as a rational combination of $h_2,\dots,h_n$, i.e. as an integral $$\dfrac{\zeta(2n-1)}{\pi^{2n}}=\int_0^1q_n(x-x^2)\frac{x^2(1-x)^2}{\sin^2\pi x}dx,$$ where $q_n$ is a unique polynomial of degree $n-2$. The first of them are:
$$\begin{align} \zeta(3)&=\frac{\pi^4}{6}h_2 \\
\zeta(5)&=\frac{\pi^6}{90}(h_2 +2h_3),\qquad \text{ i. e. } q_2(z)=\frac{1}{90}(1+2z) \quad \text{ etc. }\\
\zeta(7)&=\frac{\pi^8}{1890}(2h_2 +4h_3+3h_4)\\
\zeta(9)&=\frac{\pi^{10}}{28350}(3h_2 +6h_3+5h_4+2h_5)\\
\zeta(11)&=\frac{\pi^{12}}{935550}(10h_2 +20h_3+17h_4+8h_5+2h_6)\\
\zeta(13)&=\frac{\pi^{14}}{638512875}(\color{blue}{691}h_2 +1382h_3+1180h_4+574h_5+175h_6+30h_7)\\ \end{align}$$
Experimentally, in $\dfrac{\zeta(2n-1)}{\pi^{2n}}$ the last coefficient (i.e. the one of $h_n$ and the leading term of $q_n$) is $\dfrac{2^{2n-2}}{(2n)!}$ and the one preceding it is $\dfrac{n(n-2)}6\dfrac{2^{2n-2}}{(2n)!}$, while for the first coefficient (equally, the constant term of $q_n$), the occurrence of $\color{blue}{691}$ in the expression for $\zeta(13)$ suggests that it involves the Bernoulli number $B_{2n-2}$.
Any ideas about these polynomials?
Finally, for $k\geqslant 3$ there does not seem to exist any closed form, at least not in terms of zeta values.
What about $J({3,3,3})= \int\limits_0^1\dfrac{x^3(1-x)^3}{\sin^3\pi x}dx$?
Best Answer
Using the reflection formula followed by the recurrence formula and the Beta integral representation (DLMF) \begin{align} \frac{x(1-x)}{\sin \pi x}&=\frac{1}{\pi}x(1-x)\Gamma(x)\Gamma(1-x)\\ &=\frac{1}{\pi}\Gamma(x+1)\Gamma(2-x)\\ &=\frac{2}{\pi}B(1+x,2-x)\\ &=\frac{2}{\pi}\int_0^\infty \frac{t^x}{(1+t)^3}\,dt\\ &=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{e^{u(2x-1)}}{\cosh^3u}\,du \end{align} last expression is obtained with $t=e^{2u}$. Then, for calculating \begin{equation} I_f=\int_0^1\frac{x(1-x)}{\sin \pi x}f(x)\,dx \end{equation} one can express \begin{align} I_f&=\frac{1}{2\pi}\int_0^1\int_{-\infty}^{\infty}\frac{e^{u(2x-1)}}{\cosh^3u}\,duf(x)\,dx\\ &=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{e^{-u}}{\cosh^3u}\,du\int_0^1e^{2ux}f(x)\,dx\\ &=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{e^{-u}}{\cosh^3u}F(u)\,du \label{eq:intf} \end{align} assuming that the change of the integration order is valid and denoting \begin{equation} F(u)=\int_0^1e^{2ux}f(x)\,dx \end{equation} If $F(u)$ is analytic in the half-plane $\Im(u)>0$ and $\left|F(u)\right|=o\left(\frac{e^{4u}}{u} \right)$ for $\left|u\right|\to \infty$,$I_f$ is evaluated by integrating along the real axis closed by the semi-large circle $\Im(u)>0$, using the residue method. Poles are situated at $u_n=i(2n+1)\pi/2$ with $n=0,1,2...$. Residues are $1/2F''(i(2n+1)\pi/2)-F'(i(2n+1)\pi/2)$, where $F''(z)$ and $F'(z)$ are respectively the first and second derivative of $F(z)$. As the half-circle contribution vanishes, it comes \begin{equation} I_f=i\sum_{n=0}^\infty\left[\frac{1}{2}F''(i(2n+1)\frac{\pi}{2})-F'(i(2n+1)\frac{\pi}{2})\right] \end{equation} When $f(x)=1$, to express the original integral, $F(u)=\frac{e^{2u}-1}{2u}$, a simple calculation shows, as expected, that \begin{equation} I=\frac{8}{\pi^3}\sum_{n=0}^\infty \frac{1}{(2n+1)^3}=\frac{7\zeta(3)}{\pi^3} \end{equation}
Another expression for the result is obtained by derivation under the integral \begin{equation} I_f=-2i\sum_{n=0}^\infty \int_0^1x(1-x)f(x)e^{i(2n+1)\pi x}\,dx \end{equation} For a real function $f$, as the summation should be real, it is sufficient to keep the imaginary contribution to the integral: \begin{equation} I_f=2\sum_{n=0}^\infty \int_0^1x(1-x)f(x)\sin\left( (2n+1)\pi x \right)\,dx \end{equation}
Now, suppose that a function $f_p$ is known such that the integrals \begin{equation} J_p=\int_0^1x(1-x)f_p(x)\sin\left( (2n+1)\pi x \right)\,dx=\frac{A_p}{(2n+1)^{2p+1}} \end{equation} which gives the relation \begin{equation} I_{f_p}=2A_p\sum_{n=0}^\infty\frac{1}{(2n+1)^{2p+1}}=2A_p\left( 1-2^{-2p-1} \right)\zeta(2p+1) \end{equation} Denoting the function $Q^0(x)=x(1-x)f_p(x)$ and $Q^1,Q^2(x)$ its first and second antiderivative. Two successive integrations by parts can be performed: \begin{align} J_p&=-(2n+1)\pi \int_0^1Q^1(x)\cos\left( (2n+1)\pi x \right)\,dx\\ &=(2n+1)\pi\left[Q^2(1)+Q^2(0)\right]-(2n+1)^2\pi^2\int_0^1Q^2(x)\sin\left( (2n+1)\pi x \right)\,dx \label{eq:jp} \end{align} The free parameters in $Q^2(x)$ can be chosen in order that $Q^2(1)=Q^2(0)=0$. With \begin{equation} Q^2(z)=\int_0^zdt\int_0^tQ^0(u)\,du+az+b \end{equation} one may chose $b=0$ and $a=-\int_0^1\,dt\int_0^tQ^0(u)\,du$. Thus \begin{equation} Q^2(z)=\int_0^zdt\int_0^tQ^0(u)\,du-z\int_0^1\,dt\int_0^tQ^0(u)\,du \end{equation} If $f_p(x)$ is a polynomial, then $Q^2(x)$ also. By construction, $x=0$ and $x=1$ are among its roots. It can be written as \begin{equation} Q^2(x)=x(1-x)f_{p+1}(x) \end{equation} or \begin{equation} f_{p+1}(x)=\frac{\int_0^xdt\int_0^tu(1-u)f_p(u)\,du-x\int_0^1\,dt\int_0^tu(1-u)f_p(u)\,du}{x(1-x)} \end{equation} $J_p$ can be written as \begin{equation} J_p=-(2n+1)^2\pi^2\int_0^1x(1-x)f_{p+1}(x)\sin\left( (2n+1)\pi x \right)\,dx \end{equation} One obtain \begin{equation} \int_0^1x(1-x)f_{p+1}(x)\sin\left( (2n+1)\pi x \right)\,dx=-\frac{1}{\pi^2}\frac{A_p}{(2n+1)^{2p+3}} \end{equation} and thus \begin{equation} I_{f_{p+1}}=2A_{p+1}\left( 1-2^{-2p-3} \right)\zeta(2p+3) \end{equation} with \begin{equation} A_{p+1}=-\frac{A_p}{\pi^2} \end{equation} Starting from $f_1(x)=1$ one obtains \begin{align} f_2(x)&=\frac{1}{12}(x^2-x-1)\\ f_3(x)&=\frac{1}{360}(x^4-2x^3-2x^2+3x+3)\\ f_4(x)&=\frac{1}{20160}(x^6-3x^5-3x^4+11x^3+11x^2-17x-17)\\ ... \end{align} which gives \begin{align} \int_0^1\frac{x(1-x)}{\sin\pi x}f_2(x)\,dx&=-\frac{31}{4}\frac{\zeta(5)}{\pi^5}\\ \int_0^1\frac{x(1-x)}{\sin\pi x}f_3(x)\,dx&=\frac{127}{16}\frac{\zeta(7)}{\pi^7}\\ \int_0^1\frac{x(1-x)}{\sin\pi x}f_4(x)\,dx&=-\frac{511}{64}\frac{\zeta(9)}{\pi^9}\\ ... \end{align} Starting from $f_1(x)=x(3-x)$, other series can be obtained. For example \begin{align} &f_2(x)=\frac{1}{60}(2x^4-10x^3+5x^2+5x+5)\\ &f_3(x)=-\frac{1}{5040}(3x^6-21x^2(x^3-x^2-x-1)-49(x+1))\\ &\int_0^1\frac{x(1-x)}{\sin\pi x}f_2(x)\,dx=-\frac{381}{4}\frac{\zeta(7)}{\pi^7}\\ &\int_0^1\frac{x(1-x)}{\sin\pi x}f_3(x)\,dx=-\frac{1533}{16}\frac{\zeta(9)}{\pi^9} \end{align} Other starting points can be obtained by choosing other members of the list proposed in the question above. For example, starting from $f_1(x)=x^2P_3(x)$ above leads to an apparent different expression for $\zeta(9)$: \begin{align} &f_2(x)=-\frac{1}{56}(x^6-27x^5+57x^4-13x^3-13x^2-13x-13)\\ &\int_0^1\frac{x(1-x)}{\sin\pi x}f_2(x)\,dx=-\frac{22995}{8}\frac{\zeta(9)}{\pi^9} \end{align} Obtained polynomials are not of the form $x^pP_p(x)$ as discussed in the question, however the above method may perhaps be adapted in this case.
EDIT 04/06/2017 : (sorry for the length of this answer...)
One may characterize more precisely the family of these polynomials. It helps to symmetrize the expressions: \begin{equation} I_f=\int_0^1\frac{x(1-x)}{\sin \pi x}f(x)\,dx=\frac{1}{8}\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}g(y)\,dy \end{equation} with $g(y)=f(x)$ and $x=(1+y)/2$. In this form it is clear that odd contribution of the polynomial $g(y)$ vanishes. The same symmetrization for the proposed decomposition above reads: \begin{equation} I_f=\frac{(-1)^n}{4}\sum_{n=0}^\infty \int_{-1}^1(1-y^2)g(y)\cos\left( (2n+1) \frac{\pi y}{2} \right)\,dy \end{equation} One may adapt the method developed above. If $g_p(y)$ is an even polynomial such as \begin{equation} \sum_{n=0}^\infty \int_{-1}^1(1-y^2)g_p(y)\cos\left( (2n+1) \frac{\pi y}{2} \right)\,dy=\frac{A_p}{(2n+1)^{2p+3}} \end{equation} then, by integrating twice by part, the polynomial \begin{equation} g_{p+1}(y)=\frac{\int_{-1}^ydt\int_{-1}^t(1-u^2)g_p(u)\,du-\frac{y+1}{2}\int_{-1}^1\,dt\int_{-1}^t(1-u^2)g_p(u)\,du}{1-y^2} \end{equation} is such that \begin{equation} \sum_{n=0}^\infty \int_{-1}^1(1-y^2)g_{p+1}(y)\cos\left( (2n+1)\frac{\pi y}{2} \right)\,dy=\frac{A_{p+1}}{(2n+1)^{2p+5}} \end{equation} with $A_{p+1}=-4A_p/\pi^2$. thus \begin{equation} \frac{1}{8}\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}g_{p+1}(y)\,dy=2A_{p+1}\left( 1-2^{-2p-3} \right)\zeta(2p+5) \end{equation} One may show that $g(y)$ is an even polynomial of $y$. For $g_0(y)=1$ one has, as expected \begin{equation} \frac{1}{8}\int_{-1}^1\frac{1-y^2}{\cos\pi y/2}\,dy= \frac{7\zeta(3)}{\pi^3} \end{equation} Then, the recurrence above produces a series of even polynomials $g_p(y)$ of degree $2p$ giving successive integral expressions for $\zeta(2p+3)$. Due to the parity remark, one can conclude that any polynomial $Q(y)$, with its even power coefficient identical to that of $g_p(y)$, is such that \begin{equation} \int_{-1}^1\frac{1-y^2}{\cos\pi y/2}Q(y)\,dy=\left( -1 \right)^p8\left( 2^{2p+3}-1 \right)\frac{\zeta(2p+3)}{\pi^{2p+3}} \end{equation} The condition reads \begin{equation} Q(y)+Q(-y)=2g_p(y) \end{equation} The first polynomials (written with $Y=y^2$) are: \begin{align} g_0(y)&=1\\ g_1(y)&=\frac{1}{12}\left( Y-5 \right)\\ g_2(y)&=\frac{1}{360}\left( Y^2-14Y+61 \right)\\ g_3(y)&=\frac{1}{20160}\left( Y^3-27Y^2+323Y-1385 \right)\\ g_4(y)&=\frac{1}{1814400}\left( Y^4-44Y^3+1006Y^2-11804Y+50521 \right)\\ g_5(y)&=\frac{1}{239500800}\left( Y^{5}-65Y^4+2410Y^3-53954Y^2+631621Y-2702765\right)\\ g_6(y)&=\frac{1}{43589145600}\left(Y^6-90Y^5+4915Y^4-178268Y^3+3980887Y^2-46590634Y+199360981 \right) \end{align} In terms of the non-symmetrized function, any polynomial of the form \begin{equation} f(x)=g_p\left( 2x-1 \right)+P(2x-1) \end{equation} where $P(z)$ is an arbitrary odd polynomial, gives a result proportional to $\zeta(2p+1)$ when integrated as in $I_f$ defined above.