$\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...
