Find $\lim\limits_{n \to \infty} \int_0^1 f_n(x) \, dx$ with $f_0(x) = x$ and $f_{n+1}(x) = \sin (\pi f_n(x))$

$\newcommand{\Z}{{\mathbb Z}}\newcommand{\Q}{{\mathbb Q}} \newcommand{\N}{{\mathbb N}} \newcommand{\R}{{\mathbb R}} \newcommand{\eps}{\varepsilon} \newcommand{\C}{{\mathbb C}} \newcommand{\ds}{\displaystyle} \newcommand{\BB}{{\mathcal B}}$ My answer is based on Winther's comments and exploits the idea presented also in metamorphy's answer. It took me a long time to find the answer -- also I consulted a lot of literature on maps of an interval (Best: The book of Collet/Eckmann), but it didn't help much for a solution to this specific question.

Claim 1: The limit of $\int_0^1 f_n(x)\,dx$ is $\int_0^1 x\,w(x)\,dx$, where $w$ is the limit of the sequence of functions $w_n:[0,1]\to\R,n\in\N$, determined by $w_0\equiv1$ and \begin{equation}\label{eqwn}\tag{1}w_{n+1}(x)=\frac1{\pi\sqrt{1-x^2}} \left(w_n(\tfrac1\pi\arcsin(x))+w_n(1-\tfrac1\pi\arcsin(x))\right)\mbox{ for }x\in[0,1]. \end{equation} The value of $\int_0^1 x\,w(x)\,dx$ can be numerically determined by approximating $w$. The first 11 decimals are 0.46796294994. (to be improved)

Edit: I added a second proof that also provides more information on the functions $w_n,w$. In particular, it turns out that $\sqrt{x(1-x)}w(x)$ can be extended to a holomorphic function on any disk with center $0.5$ and radius $<1.5$.
Additionally, it gives a theoretical basis for better approximations of $\int_0^1 x\,w(x)\,dx$. The first 101 decimals are \begin{equation}\nonumber\begin{array}{l} 0.4679629499441669225983425685149731926726\\ \ \ \ 4434670336201394146766196652779776580187\\ \ \ \ 502477579258564076297. \end{array}\end{equation} This confirms and improves the value given in metamorphy's comment.

The functions determined by (\ref{eqwn}) satisfy, for all continuous functions $f$, $$\int_0^1f(\sin(\pi z))w_n(z)\,dz= \int_0^{1/2} f(\sin(\pi z))(w_n(z)+w_n(1-z))\,dz= \int_0^1f(x)w_{n+1}(x)\,dx;$$ here the substitution $z=\frac1\pi\arcsin(x)$ was used. Since we also have $f_{n+1}(x)=f_n(\sin(\pi x))$ for all $n$, this leads to $\int_0^1f_n(x)\,dx=\int_0^1x\,w_n(x)\,dx$ for all $n$ as can be proved recursively. Therefore the first part of Claim 1 follows from
Lemma 1: The functions $x\mapsto \sqrt{x(1-x)}w_n(x)$ can be continued to continuous functions on $[0,1]$ and for $n\to\infty$, the sequence converges uniformly on $[0,1]$.
Definition: The limit of the sequence in the Lemma is written $\sqrt{x(1-x)}w(x)$ with some continuous function $w$ on $]0,1[$ which is integrable over $[0,1]$.
Remark: $w$ is the density of the invariant measure mentioned in Winther's comment. According to the literature, it is uniquely determined. It satisfies $$\begin{array}c \ds\int_0^1f(\sin(\pi x))w(x)\,dx=\int_0^1f(x)w(x)\,dx\mbox{ for all continuous functions and}\\ \ds w(x)=\frac1{\pi\sqrt{1-x^2}} \left(w(\tfrac1\pi\arcsin(x))+w(1-\tfrac1\pi\arcsin(x))\right)\mbox{ for }x\in]0,1[. \end{array}$$ Proof: The problem is that the denominator of (\ref{eqwn}) vanishes at $x=1$ and that $\arcsin$ has an infinite derivative at $x=1$. In a first step, we transform the sequence $w_n$ attached to $\arcsin$ into the analogous sequence attached to a function that has a bounded derivative. We use again the functions $$h(t)=\sin^2\left(\frac\pi2\,t\right), h^{-1}(t)=\frac2\pi\arcsin(\sqrt{ t})$$ that conjugated the map $x\to4x(1-x)$ to the tent map. They will conjugate the map $x\to\sin(\pi x)$ to a "deformed tent map".

Precisely, we put $v_n(t)=h'(t)w_n(h(t))$ so that we have $$v_{n+1}(t)=\frac{h'(t)}{\pi \sqrt{1-h(t)^2}} \left(w_n(\tfrac1\pi\arcsin(h(t)))+w_n(1-\tfrac1\pi\arcsin(h(t)))\right)$$ and hence \begin{equation}\label{eqvn}\tag{2} v_{n+1}(t)={\phi'(t)}\left(v_n(\phi(t))+v_n(1-\phi(t))\right) \end{equation} where $\phi(t)=h^{-1}(\frac1\pi\arcsin(h(t)))$. Here we used that $$\phi'(t)= (h^{-1})'(\tfrac1\pi\arcsin(h(t)))\cdot\frac1{\pi\sqrt{1-h(t)^2}}\cdot h'(t)= \frac{h'(t)}{h'(\phi(t))\pi\sqrt{1-h(t)^2}}$$ and that $h(1-t)=1-h(t)$, $h'(1-t)=h'(t)$. We leave it to the reader to show (using series expansions) that $\phi$ can be extended to a continuous function mapping bijectively $[0,1]$ to $[0,\frac12]$ and that the extension is real analytic, in particular $\phi'(0)=\frac1{\sqrt\pi}$ and $\phi'(1)=\frac{\sqrt2}\pi$. The inverse $\phi^{-1}$ of $\phi$ is the restriction to $[0,\frac12]$ of the map $\psi:[0,1]\to[0,1]$ conjugated to $\sin(\pi x)$ defined by $\psi(x)=h^{-1}(\sin(\pi h(x)))$. This map $\psi$ is a modified tent map and has the important advantage (compared with $\sin(\pi x)$) that the (one-sided) derivatives vanish nowhere. Here is a graph of $\psi$.

In a second step we show that the sequence $v_n$, $n\in\N$, of functions satisfying (\ref{eqvn}) and $v_0(t)=h'(t)$ converge uniformly to some limit we call $v$. By induction using that $\phi$ is real analytic, it follows that all $v_n$ are of class $C^1$. By (\ref{eqvn}), it also follows that $\int_0^1v_n(t)\,dt=1$ for all $n$.

Consider now the functions $d_n=v_{n+1}-v_n$. They satisfy similar to (\ref{eqvn}) $$d_{n+1}(t)={\phi'(t)}\left(d_n(\phi(t))+d_n(1-\phi(t))\right)$$ and $\int_0^1d_n(t)\,dt=0$. Unfortunately, this recursion seems not to permit sufficiently good estimates of $||d_n||=\max_{t\in[0,1]}|d_n(t)|$. Surprisingly the differentiated recursion \begin{equation}\label{eqdn}\tag{3} d_{n+1}'(t)={\phi'(t)}^2\left(d_n'(\phi(t))-d_n'(1-\phi(t))\right)+ {\phi''(t)}\left(d_n(\phi(t))+d_n(1-\phi(t))\right) \end{equation} is more convenient.

For the estimate, we use that $||d_n||\leq||d_n'||$. This follows from the facts that $\int_0^1d_n(t)\,dt=0$, hence $d_n$ vanishes in some point $t_n\in[0,1]$ and therefore $d_n(t)=\int_{t_n}^{t}d_n'(\tau)\,d\tau$. This implies with (\ref{eqdn}) that $$|d_{n+1}'(t)|\leq2||d_n'||(\phi'(t)^2+|\phi''(t)|)\mbox{ for all }t\in[0,1].$$ Now numerical estimates (see below) show that $$2(\phi'(t)^2+|\phi''(t)|)\leq \gamma<0.934.$$ This estimate is, of course, crucial for the sequel. Observe that estimating $\phi'(t)^2+|\phi''(t)|$ by $||\phi'||^2+||\phi''||$ would not be good enough because the latter quantity is $0.507...$. Thus we have found that \begin{equation}\nonumber||v_{n+1}-v_n||\leq||v_{n+1}'-v_n'||\leq \gamma||v_n'-v_{n-1}'||\mbox{ for all }n\in\N. \end{equation} This implies that \begin{equation}\nonumber||v_m-v_n||\leq||v_{m}'-v_n'||\leq\frac1{1-\gamma} \gamma^{n}||v_1'-v_0'||\mbox{ for all positive integers }n\leq m. \end{equation} This means that $v_n$, $n\in\N$, and $v_n'$, $n\in\N$, are Cauchy sequences with respect to uniform convergence. Therefore they converge to continuous functions on $[0,1]$. We call $v$ the limit of $v_n$; the one of $v_n'$ is then $v'$: $v$ is of class $C^1$.

Going back to the $w_n$ by $w_n(x)=v_n(h^{-1}(x))/h'(h^{-1}(x))=v_n(h^{-1}(x))/\sqrt{x(1-x)}$ now implies Lemma 1. We have $w(x)=v(h^{-1}(x))/\sqrt{x(1-x)}$ (or $v(x)=w(h(x))h'(x)$). Here are graphs of $v$ and $\frac1\pi w$ (determined as indicated below).

Numerical calculations: Except at the endpoints, the derivatives $\phi'$ and $\phi''$ were calculated by the formulas $$\phi'(t)=\frac{\phi(t+h/2)-\phi(t-h/2)}h,\ \phi''(t)=\frac{\phi(t+h)-2\phi(t)+\phi(t-h)}{h^2}$$ using a sufficiently small $h$. Here is a table of values. \begin{equation}\nonumber\begin{array}{r|c|c|c|c|c|c|c|c|c|c|c} x&0 & 0.1 & 0.2 & 0.3 & 0.4 & 0.5 & 0.6 & 0.7 & 0.8 & 0.9 & 1 \\\hline \phi'(x) & 0.564 & 0.560 & 0.547 & 0.530 & 0.511 & 0.493 & 0.478 & 0.466 & 0.457 & 0.452 & 0.450 \\\hline -\phi''(x)&0& 0.0901 & 0.155 & 0.186 & 0.186 & 0.167 & 0.138 & 0.104 & 0.0693 & 0.0345 & 0 \end{array}\end{equation}

The wanted limit $$L:=\lim_{n\to\infty}\int_0^1f_n(x)\,dx=\lim_{n\to\infty}\int_0^1 x\,w_n(x)\,dx=\int_0^1 x\,w(x)\,dx$$ can best be calculated (because $w(x)$ has singularities at 0 and 1) by substituting again $x=h(t)=\sin^2\left(\frac\pi2 t\right)$. We find $$L=\int_0^1 h(t)v(t)\,dt$$ with $v(t)=w(h(t))h'(t)$ the above limit of the $v_n$ which is a $C^1$ function.

The idea for the evaluation of $L$ is to approximate $v_n$ by piecewise linear function $\ell_n$, to use (\ref{eqvn}) to find a recursion for $\ell_n$ and in the limit to find an approximation for $v$ that is used to calculate $L$. We find the above graphs for $v$ and then for $w$ and $$L=0.46796294994...$$

Edit: We want to find out more about the solutions $v$ of the functional equation (compare (\ref{eqvn})) \begin{equation}\label{eqv}\tag{4} v(t)=(Tv)(t):={\phi'(t)}\left(v(\phi(t))+v(1-\phi(t))\right) \end{equation} where $\phi(t)=h^{-1}(\frac1\pi\arcsin(h(t)))$, $h(t)=\sin^2\left(\frac\pi2 t\right)$. Recall that $\phi$ can be extended to a continuous function mapping bijectively $[0,1]$ to $[0,\frac12]$ and that the extension is real analytic, in particular $\phi'(0)=\frac1{\sqrt\pi}$ and $\phi'(1)=\frac{\sqrt2}\pi$. We also have $\frac{\sqrt2}\pi\leq\phi'(t)\leq\frac1{\sqrt\pi}$ for all $t\in[0,1]$. Recall finally that $\phi=\left(\psi\mid_{[0,\frac12]}\right)^{-1}$ where $\psi(t)=h^{-1}(\sin(\pi(h(t))))$, $0\leq t\leq1$.

Claim 2: a) For all continuous functions $f:[0,1]\to\C$ we have $$\int_0^1 (Tf)(t)\,dt= \int_0^1f(t)\,dt\mbox{ and } \int_0^1 |(Tf)(t)|\,dt\leq \int_0^1|f(t)|\,dt.$$ b) If $f\neq0$ changes sign, i.e. there exist $t_1,t_2$ such that $f(t_1),f(t_2)\neq0$ and $f(t_1)/f(t_2)\not\in\R_+$, then there exists a positive integer $N$ such that $$\int_0^1 |(T^nf)(t)|\,dt< \int_0^1|f(t)|\,dt\mbox{ for }n\geq N.$$ c) If $f\neq0$ does not change sign then there exist $N\in\N$ and $c\in\C$ such that $c\,T^nf$ is positive throughout $[0,1]$ for $n\geq N$.

Before proving it, we give a corollary:

Corollary 1: Except for a scalar factor, there is at most one nonzero continuous solution $v$ of (\ref{eqv}). If it exists, it does not change sign and does not vanish on $[0,1]$.

Proof of the Corollary: Any solution of (\ref{eqv}) does not change sign according to Claim 2 b) since $T^nv=v$. Then Claim 2 c) implies that $v$ does not vanish.

Therefore, we can assume that two nonzero solutions $v,\tilde v$ both satisfy $\int_0^1v(t)\,dt=\int_0^1\tilde v(t)\,dt=1$. Then their difference $d=v-\tilde v$ has $\int_0^1d(t)\,dt=0$. It also is a solution of (\ref{eqv}) and the already proved statements imply a contradiction if $d\neq0$.

Proof of Claim 2: For the first part, we estimate directly $$\begin{array}{rcl} \ds\int_0^1 |(Tf)(t)|\,dt&\leq&\ds\int_0^{1}\phi'(t)|f(\phi(t))|dt+\int_0^{1}\phi'(t)|f(1-\phi(t))|dt\\ &\leq&\ds\int_0^{1/2}|f(z)|dz+\int_{1/2}^{1}|f(z)|dz= \int_0^1|f(z)|\,dz. \end{array}$$ Without the absolute value signs, we clearly have equality everywhere.

For b) and c), we rewrite the definition of $T$ $$(Tf)(t)=\sum_{\psi(s)=t}\frac1{|\psi'(s)|}f(s)\mbox{ for }t\in[0,1],$$ where the summation is taken over $s$ such that $\psi(s)=t$; here $s$ can be $\phi(t)$ or $1-\phi(t)$. By induction, this allows to express $T^n$ in a similar way: \begin{equation}\label{eqTn}\tag{5} \left(T^nf\right)(t)=\sum_{\psi^n(s)=t}\frac1{|(\psi^n)'(s)|}f(s)\mbox{ for }t\in[0,1]. \end{equation} Here $\psi^n$ denotes the composition $\underbrace{\psi\circ\psi\circ\cdots\circ\psi}_{n\text{ copies of }\psi}$. Now the graph of $\psi^n$ is similar to that of $f_n$ described in the question. $[0,1]$ can be divided into $2^n$ intervals $$I_1 = [a_0, a_1],\quad I_2 = [a_1, a_2],\quad \dotsc,\quad I_{2^n} = [a_{2^n-1}, a_{2^n}]$$ such that:
1. $\psi^n(a_k) = 1$ if $k$ is odd and $\psi^n(a_k) = 0$ if $k$ is even;
2. $\psi^n$ is monotonic in each $I_k$, increasing if $k$ is odd and decreasing if $k$ is even.

As the derivative of $\psi^n$ is between $\sqrt\pi^n$ and $\left(\frac\pi{\sqrt2}\right)^n$, the length of each subinterval is at most $\pi^{-n/2}$. Therefore, for any $t\in[0,1]$, the distance between two subsequent elements of $\{s\mid\psi^n(s)=t\}$ is at most $2\pi^{-n/2}$. Hence \begin{equation}\label{dist}\tag{6} \mbox{dist}(z,\{s\mid\psi^n(s)=t\})\leq 2 \pi^{-n/2}\mbox{ for any }t,z\in[0,1],\,n\geq1. \end{equation} For the proof of b) consider now a continuous function $f$ such there exist $z_1,z_2$ with $f(z_1),f(z_2)\neq0$ and $f(z_1)/f(z_2)\not\in\R^+$. Then by (\ref{dist}) and continuity, for all sufficiently large $n$ and any $t\in[0,1]$, there exist $s_1,s_2$ such that $\psi^n(s_1)=\psi^n(s_2)=t$, $f(s_1),f(s_2)\neq0$ and $f(s_1)/f(s_2)\not\in\R^+$. This implies that $$\left|\sum_{\psi^n(s)=t}\frac1{|(\psi^n)'(s)|}f(s)\right|< \sum_{\psi^n(s)=t}\frac1{|(\psi^n)'(s)|}|f(s)|$$ and hence by (\ref{eqTn}) that $|(T^nf)(t)|<(T^n|f|)(t)$ for all sufficiently large $n$. Integrating from 0 to 1 this yields using a) that $\int_0^1|(T^nf)(t)|dt<\int_0^1|f(t)|dt$: b) is proved.

For the proof of c) consider a nonzero continuous function $f$ not changing sign. Except for a constant factor, we can assume that $f(z)\geq0$ for all $z\in[0,1]$ and $f(z)>0$ for some $z$. By (\ref{dist}) and continuity, for all sufficiently large $n$ and all $t\in[0,1]$ there exists $s$ with $\psi^n(s)=t$ such that $f(s)>0$. Using (\ref{eqTn}), this implies that $(T^nf)(t)>0$ for all sufficiently large $n$ and $t\in[0,1]$ and c) is proved.

For $\delta>0$, consider now the compact set $C_\delta\subset\C$ consisting of all complex numbers having a distance at most $\delta$ from the interval $[0,1]$ and the Banach space $\BB_\delta$ containing all continuous functions on $C_\delta$ that are holomorphic in the interior of $C_\delta$, equipped with the $\sup$-norm $||\cdot||_\infty$. If we choose a sufficiently small $\delta$ then $\phi$ is in $\BB_\delta$ and $\gamma=||\phi'||_\infty<1$. Since $\phi(t),1-\phi(t)\in[0,1]$ for $t\in[0,1]$, this implies that $$\phi(z),1-\phi(z)\in C_{\gamma\delta}\mbox{ if }z\in C_\delta.$$ Therefore $T$ can be defined on $\BB_\delta$ - without changing the name - by defining for any $f\in\BB_\delta$ $$(Tf)(z)=\phi'(z)(f(\phi(z))+f(1-\phi(z))),\ \ z\in C_\delta.$$ Clearly $T:\BB_\delta\to\BB_\delta$ is a bounded linear operator. The fact that $\gamma<1$ implies the most important property of $T$ on $\BB_\delta$:
Claim 3: $T:\BB_\delta\to\BB_\delta$ is compact, i.e. it maps bounded sets to relatively compact sets.
Proof: Consider a subset $M\subset\BB_\delta$ bounded by $K$, i.e. $||f||_\infty\leq K$ for every $f\in M$. We have to show that $TM$ is relatively compact in $\BB_\delta$. By the Arzela-Ascoli Theorem, it is sufficient to show that $TM$ is uniformly bounded and that the elements of $TM$ have a common Lipschitz constant. This is not difficult. If $f\in M$ then $||Tf||_\infty\leq 2||\phi'||_\infty||f||_\infty\leq 2||\phi'||_\infty K.$ By Cauchy's inequalities, $|f'(t)|\leq \frac1{(1-\gamma\delta)}||f||_\infty$ for $t\in C_{\gamma\delta}$ and hence $$\begin{array}{rcl}|(Tf)'(z)|&\leq &|\phi'(z)|^2(|f'(\phi(z))|+|f'(1-\phi(z))|)+|\phi''(z)|(|f(\phi(z))|+|f(1-\phi(z))|)\\ &\leq& \frac{2||\phi'||_\infty^2}{(1-\gamma\delta)}K+2||\phi''||_\infty K=:L \end{array}$$ for $z\in C_\delta$ because $\phi(z),1-\phi(z)\in C_{\gamma\delta}$. It is well known that this implies that all $g\in TM$ have the common Lipschitz constant $L$. Claim 3 is proved.

Using the spectral theory of compact operators, we are now in a position to prove
Theorem: The Equation $Tv=v$ has a unique solution $v\in\BB_\delta$ satisfying $\int_0^1v(t)dt=1$. For every $f\in\BB_\delta$, we have $$\lim_{n\to\infty}T^nf=\int_0^1f(t)\,dt\cdot v.$$ In particular, the sequence $v_n$, $n\in\N$, defined by $v_0=h'$ and $v_{n+1}=Tv_n$ (also used in the first proof) which remains in $\BB_\delta$ since $v_0\in\BB_\delta$ converges to $v$ uniformly of $C_\delta$.
Proof of the Theorem: The spectral radius of $T$ is at least $r(T)=1$. If $r(T)<1$ then $||T^n||_\infty\to0$ as $n\to\infty$ and hence $T^nf\to0$ for all $f\in\BB_\delta$. This is impossible, because $\int_0^1(T^nv_0)(t)dt=1$ for $v_0=h'$ and all $n$. The spectral radius of $T$ is at most 1. Otherwise, there is an eigenvalue $\lambda$ of $T$ with $|\lambda|>1$. Consider a corresponding eigenfunction $f$ of $T$. Then we have $\int_0^1|(Tf)(z)|dz=|\lambda| \int_0^1|f(z)|dz > \int_0^1|f(z)|dz$ which is a contradiction to Claim 2 a). So we have proved that $r(T)=1$.

Consider now an eigenvalue $\lambda$ of $T$ with $|\lambda|=1$ and a corresponding eigenfunction $f$. As above, it satisfies $|T^nf(z)|=|f(z)|$ for all $n\in\N$. By Claim 2 b), $f$ does not change sign and we can assume that $f(z)\geq0$ for all $z\in[0,1]$. Then $\int_0^1(Tf)(z)dz=\lambda\int_0^1f(z)dz$ implies that $\lambda=1$. So $\lambda=1$ is an eigenvalue of $T$. As seen in the Corollary 1 below Claim 2, there is a unique (except for scalar factors) eigenfunction $v$ of $T$ corresponding to this eigenvalue. There cannot be a generalized eigenvector either: Otherwise there is a $g\in\BB_\delta$ such that $Tg=g+v$. Integrating from 0 to 1 using Claim 2 a) this implies $\int_0^1v(t)dt=0$: a contradiction. Therefore the generalized eigenspace of $T$ corresponding to the eigenvalue 1 has dimension 1. This proves the first part of the Theorem.

As we have seen above, 1 is the only eigenvalue of $T$ with modulus 1. All its other eigenvalues $\lambda$ must hence satisfy $|\lambda|<1$. Since by Claim 2 a) $\int_0^1f(z)dz=\int_0^1(Tf)(z)dz=\lambda\int_0^1f(z)dz$ for any corresponding eigenfunctions $f$, we must have $\int_0^1f(z)dz=0$ for corresponding eigenfunctions. This makes it useful to decompose $\BB_\delta=\langle v\rangle\oplus S$ where $S$ is the closed subspace of $\BB_\delta$ consisting of all $g$ such that $\int_0^1g(z)dz=0$: Write any $f\in\BB_\delta$ as $f=\int_0^1f(z)dz\cdot v+g$ with some $g\in S$. By Claim 2 a) we have $T(S)\subset S$, i.e. $S$ is an invariant subspace of $\BB_\delta$. The restriction $T\mid_S:S\to S$ is again a compact linear operator. As seen above, it cannot have eigenvalues $\lambda$ with $|\lambda|\geq1$. Hence $r(T\mid_S)<1$. Therefore, if we write any $f\in\BB_\delta$ as $f=\int_0^1f(z)dz\cdot v+g$, we have $T^nf=\int_0^1f(z)dz\cdot v+(T\mid_S)^ng$ and $||(T\mid_S)^ng||_\infty\leq||(T\mid_S)^n||_\infty\,||g||_\infty\to0$. This proves the second part of the Theorem.

We can obtain even more information about the unique solution $v$ in the Theorem. Observe that $\phi$ and $z\mapsto \phi(1-z)-\frac12$ are both odd; $\phi'$ and $z\mapsto\phi'(1-z)$ are even. This implies that $T:P\to P$, where $P$ is the closed subspace of $\BB_\delta$ consisting of all $f$ such that $f$ and $z\mapsto f(1-z)$ are even. By induction, it follows that all $\tilde v_n\in P$ if we start with $\tilde v_0\equiv1$ and define recursively $\tilde v_{n+1}=T\tilde v_n$. Hence also their limit $v$ is an element of $P$. Now a real analytic function $f$ such that $f$ is even and $z\mapsto f(1-z)$ is even must be 2-periodic!
Proposition 1: The essentially unique solution $v$ of $Tv=v$ is real analytic on the whole real axis, 2-periodic and even.

Numerically, we can determine values for $\delta$ such that $T:\BB_\delta\to\BB_\delta$ is a compact operator: We need that $\phi$ can be extended to a holomorphic function on $C_\delta$ and that the image $\phi(C_\delta)$ is contained in the interior of $C_\delta$. We find that $\delta=0.5$ is possible (see figure below). In the same way, we find that $\phi$ can be extended to a holomorphic function of the closed disk $\bar D(\frac12,0.6)$ containing the interval $[0,1]$ in its interior and that the image $\phi(\bar D(\frac12,0.6))$ is contained in the interior of $\bar D(\frac12,0.6)$ (see figure below).

The same proof as above yields that the essentially unique solution $v$ of $Tv=v$ can be extended to a holomorphic function on $C_{0.5}$ and on $\bar D(\frac12,0.6)$. Together with the symmetry and periodicity of $v$ we obtain
Proposition 2: The essentially unique solution $v$ of $Tv=v$ can be extended to a 2-periodic even holomorphic function on the strip $|\mbox{Im}(z)|<\frac12$ of width 1 around the real axis. Moreover, it can be analytically continued to the open disks $D(n+\frac12,0.6),\, n\in\Z.$

This result implies that $v$ can be represented by a single power series on the whole interval $[0,1]$ (as it is a subset of $D(\frac12,0.6)$). This suggest that $v$ can well be approximated by a polynomial on $[0,1]$. Such polynomials could be found numerically modifying the iteration $v_{n+1}=Tv_n$: Choose as $p_{n+1}$ the truncation to $M$ terms of the series determined by $\phi'(z)(p_n(\phi(z))+p_n(1-\phi(z)))$.

As a variant, Proposition 2 suggests to use Fourier series, but the software I use (pari/gp) seems not to have an easy mechanism to deal with compositions of Fourier series.

Another variant is to approximate $v$ by individual polynomials on a subdivision of $[0,1]$ into small intervals and to modify the iteration $v_{n+1}=Tv_n$ to determine them recursively. I used 100 subintervals, polynomials of degree 150 and a precision of 150 decimals to find the 101 decimals of $\int_0^1xw(x)dx=\int_0^1h(z)v(z)dz$ mentioned in the beginning.

In the last part, we carry back the above results to the original functions $w_n,w$ of Claim 1.

First of all, the functions $w_n$ of Claim 1 and $v_n$ of the Theorem are related by $v_n(x)=w_n(h(x))h'(x)$ with $h(t)=\sin(\frac\pi2t)^2$. Thus the Theorem implies that the functions defined by $\sqrt{x(1-x)}w_n(x)$ converge uniformly on $[0,1]$. Observe that $v_n$ of the Theorem are not in the subspace $\mathcal P$ used before Proposition 1 and $\sqrt{x(1-x)}w_n(x)$ do not define holomorphic functions near $x=0$.

Because of Proposition 2, we can write $v(t)=\tilde w(h(t))$ with some function $\tilde w$ real analytic on $[0,1]$. A comparison with the end of the proof of Lemma 1 and the subsequent Definition shows
Corollary 2 The equation $w(\sin(\pi x))\pi|\cos(\pi(x))|=w(x)+w(1-x)$ has a unique solution such that $\int_0^1w(x)dx=1$ and $\tilde w(x)=\sqrt{x(1-x)}w(x)$ is continuous on $[0,1]$. The corresponding $\tilde w$ is real analytic on $[0,1]$.
Now $\tilde w$ satisfies \begin{equation}\label{eqwt}\tag{7} \tilde w(x)=(U\tilde w)(x):=\frac{\sqrt x}{g(x)k(\frac1\pi\arcsin(x))} (\tilde w(\tfrac1\pi\arcsin(x))+\tilde w(1-\tfrac1\pi\arcsin(x))) \end{equation} for $x\in[0,1]$, where $g(x)=\pi\sqrt{1+x}$ and $k(t)=\sqrt{t(1-t)}$. Observe that $F(x)=\frac{\sqrt x}{g(x)k(\frac1\pi\arcsin(x))}$ can be continued to an analytic function at $x=0$ as the roots of $x$ in the numerator and in the denominator cancel each other. It can also be continued to an analytic funtion at $x=1$: We write $\tfrac1\pi\arcsin(x)=\frac12+\frac1\pi\int_1^x\frac{dt}{\sqrt{1-t^2}}= \frac12+r(x)\,\sqrt{1-x}$, where $r$ is analytic at $x=1$. Thus we can write $k(x)=\sqrt{\frac14-(1-x)r(x)^2}$ which reveals that $k$ and hence also $F$ can be continued to an analytic function at $x=1$. Altogether, $F$ is real analytic on $[0,1]$.

Finally, if $f(t)$ denotes a function which is analytic near $t=\frac12$, then $f(\frac12+s)+f(\frac12-s)$ is an even function of $s$ and hence can be written $g(s^2)$ with a function $g=g(u)$ analytic at $u=0$. This implies that $f(\frac1\pi\arcsin(x))+f(1-\frac1\pi\arcsin(x))=g((1-x)r(x)^2)$ (with the above function $r$) is analytic at $x=1$. Altogether, we have shown that with $f$ real analytic on $[0,1]$, also its image $Uf$ under the operator $U$ of (\ref{eqwt}) is real analytic on $[0,1]$.

Now the $\arcsin$ is holomorphic on every domain not containing $\pm1$. Numerically, it can be shown that $x\to\frac1\pi\arcsin(x)$ and $x\to1-\frac1\pi\arcsin(x)$ map the closed disk of any radius $r<1.5$ and center $0.5$ into its interior (see figure below for $r=1.45$). Observe that $x=1$ is a ramification point of $\arcsin$ and that by continuing $x\to\frac1\pi\arcsin(x)$ analytically once around $x=1$, we arrive at $x\to1-\frac1\pi\arcsin(x)$ and vice versa.

Therefore we can define the operator $U$ on ${\mathcal E}_r$, the Banachspace of all continuous functions $f$ on the closed disk $\bar D(\frac12,r)$ holomorphic in its interior by $$(Uf)(z)=F(z)(f(\tfrac1\pi\arcsin(z))+f(1-\tfrac1\pi\arcsin(z))),\ \ z\in \bar D(\tfrac12,r)$$ where $F(x)=\frac{\sqrt x}{g(x)k(\frac1\pi\arcsin(x))}$ and $g(x)=\pi\sqrt{1+x}$, $k(t)=\sqrt{t(1-t)}$. Again, we keep the name of the analogous operator on $C[0,1]$.

As $\psi$ and $x\mapsto\sin(\pi x)$ are conjugated via $h$, the density statement (\ref{dist}) can be carried over to the latter mapping. With this, the whole previous proof concerning $T$ can be carried over to $U$. We obtain
Proposition 3 The function $\tilde w(x)=\sqrt{x(1-x)}w(x)$ of Corollary 2 can be continued analytically to a holomorphic function on the open disk $D(\frac12,1.5)$. Its restriction to $\bar D(\frac12,r)$, $0.5<r<1.5$ is the ${\mathcal E}_r$-limit of the sequence defined by $\tilde w_0\equiv1$ and $\tilde w_{n+1}=U\tilde w_n$, $n\in\N$.
As before, this result implies that $\tilde w$ can be represented by a single convergent power series on $[0,1]$. It can therefore be well approximated by polynomials on $[0,1]$ which, as before, could be determined using a modification of $\tilde w_{n+1}=U\tilde w_n$. This seems close to the method used by metamorphy for the calculation of the limit in the question.

In a final remark, I would like to mention that equation (\ref{eqwt}) permits to continue $\tilde w$ analytically to bigger domains than $D(0.5,1.5)$. This cannot be pursued any further here...

This partial answer is merely the @Winther's idea from a different point of view.

We clearly have $f_{n+1}(x)=f_n(\sin\pi x)$. Now if $w:(0,1)\to\mathbb{R}$ is integrable, then $$\int_0^1 w(x)f_{n+1}(x)\,dx=\int_0^1\overline{w}(x)f_n(x)\,dx,$$ where $$\overline{w}(x)=\frac{1}{\pi\sqrt{1-x^2}}\left[w\left(\frac{\arcsin x}{\pi}\right)+w\left(1-\frac{\arcsin x}{\pi}\right)\right].$$ Thus, by induction, if $w_0\equiv 1$ and $w_{n+1}=\overline{w_n}$, then $$\int_0^1 f_n(x)\,dx=\int_0^1 xw_n(x)\,dx.$$ Experiments suggest that $w_\infty(x)=\lim\limits_{n\to\infty}w_n(x)$ exists (it then has to satisfy the functional equation given by @Winther) and has a shape of $v(x)/\sqrt{x(1-x)}$ with $v(x)$ decreasing from $\approx 0.18$ to $\approx 0.14$; more exactly, $v(1)=v(1/2)\cdot2\sqrt{2}/\pi$ and $v(0)=v(1)/(\sqrt{\pi}-1)$. I've computed the value given in my comment above by approximating $v(x)$ with a polynomial. I've also tried analysing Fourier series of $w_n(x)$, with no substantial progress so far.