I want to solve this problem by using $\varepsilon-\delta $ definition of a limit, but failed.How to compute this limit?
$f\in C\left( \left[ 0,1 \right] \times \left[ 0,1 \right] \right) $
$$\mathop {\lim }\limits_{n \to \infty } {\left( {\frac{{\left( {2n + 1} \right)!}}{{{{\left( {n!} \right)}^2}}}} \right)^2}\int_0^1 {\int_0^1 {{{\left[ {xy\left( {1 - x} \right)\left( {1 - y} \right)} \right]}^n}f\left( {x,y} \right)} } {\text{d}}x{\text{d}}y$$
@HansLundmark has already given away the answer. To prove it we divide the region $R=\{(x,y)|0\le x\le1,0\le y\le1\}$ into $R_1=\left\{(x,y)|\frac12-r<x<\frac12+r,\frac12-r<y<\frac12+r\right\}$ and $R_2=R\setminus R_1$. Then $$\begin{align}I=\int\int_R\left[xy(1-x)(1-y)\right]^nf(x,y)d^2A&=\int\int_R\left[xy(1-x)(1-y)\right]^nf\left(\frac12,\frac12\right)d^2A\\ &\quad+\int\int_{R_1}\left[xy(1-x)(1-y)\right]^n\left[f(x,y)-f\left(\frac12,\frac12\right)\right]d^2A\\ &\quad+\int\int_{R_2}\left[xy(1-x)(1-y)\right]^n\left[f(x,y)-f\left(\frac12,\frac12\right)\right]d^2A\\ &=I_0+I_1+I_2\end{align}$$ Now, $$I_0=f\left(\frac12,\frac12\right)\left[B(n+1,n+1)\right]^2=\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2f\left(\frac12,\frac12\right)$$ Since $f(x,y)$ is continuous at $(x,y)=\left(1/2,1/2\right)$ we know that for any $\epsilon_1>0$ there is a $\delta_1(\epsilon_1)>0$ such that $\left|f(x,y)-f\left(\frac12,\frac12\right)\right|<\epsilon_1$ whenever $$\sqrt{\left(x-\frac12\right)^2+\left(y-\frac12\right)^2}<\delta_1$$ So we choose $$0<r<\min\left(\frac{\delta_1\left(\frac{\epsilon}2\right)}{\sqrt2},\frac12\right)$$ So that $$\begin{align}\left|I_1\right|&\le\int\int_{R_1}\left[xy(1-x)(1-y)\right]^n\left|f(x,y)-f\left(\frac12,\frac12\right)\right|d^2A\\ &<\int\int_{R_1}\left[xy(1-x)(1-y)\right]^n\left(\frac{\epsilon}2\right)d^2A\\ &<\int\int_{R}\left[xy(1-x)(1-y)\right]^n\left(\frac{\epsilon}2\right)d^2A\\ &=\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2\left(\frac{\epsilon}2\right)\end{align}$$ If $n\ge2$ then $\frac83n=2n+\frac23n>2n+1$ so $$\begin{align}\frac{(2n+1)!}{n\cdot2^{2n}\left(n!\right)^2}&=\frac{(2n+1)}n\prod_{k=1}^n\frac{(2k-1)(2k)}{(2k)(2k)}=\frac{(2n+1)}n\left(\frac12\right)\left(\frac34\right)\prod_{k=3}^n\frac{(2k-1)}{(2k)}\\ &=\frac{(2n+1)}{\frac83n}\prod_{k=3}^n\frac{(2k-1)}{(2k)}\le1\end{align}$$ Also the continuity of $f$ over a closed set $R$ implies that it is bounded on that set so there is some $M$ such that $\left|f(x,y)\right|\le M$ for all $(x,y)\in R$, so by the triangle inequality $$\left|f(x,y)-f\left(\frac12,\frac12\right)\right|\le2M$$ For $(x,y)\in R\supset R_2$. Also $xy(1-x)(1-y)$ attains its maximum absolute value in $R_2$ at $\left(\frac12\pm r,\frac12\right)$ and at $\left(\frac12,\frac12\pm r\right)$ of $\frac1{16}\left(1-4r^2\right)$. Then $$\begin{align}\left|I_2\right|&\le\int\int_{R_2}\left[xy(1-x)(1-y)\right]^n\left|f(x,y)-f\left(\frac12,\frac12\right)\right|d^2A\\ &\le\int\int_{R_2}\frac1{2^{4n}}\left(1-4r^2\right)^n(2M)d^2A\\ &=\int\int_{R_2}\left[\frac{(2n+1)!}{n\cdot2^{2n}\left(n!\right)^2}\right]^2\left[\frac{n\cdot\left(n!\right)^2}{(2n+1)!}\right]^2\left(1-4r^2\right)^n(2M)d^2A\\ &\le\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2n^2\left(1-4r^2\right)^n(2M)\int\int_{R_2}d^2A\\ &\le\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2n^2\left(1-4r^2\right)^n(2M)\int\int_{R}d^2A\\ &=\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2(2M)n^2\left(1-4r^2\right)^n\end{align}$$ We may establish via L'Hopital's Rule that for $0<r<\frac12$ $$\lim_{n\rightarrow\infty}g(n,r)=\lim_{n\rightarrow\infty}n^2\left(1-4r^2\right)^n=0$$ Which means that for any $\epsilon_2>0$ there is an $N_2\left(\epsilon_2,r\right)$ such that $g(n,r)<\epsilon_2$ whenever $n>N_2\left(\epsilon_2,r\right)$ so we choose $$n>\max\left(N_2\left(\frac{\epsilon}{4M},r\right),2\right)$$ Which yields the desired $$\left|I_2\right|\le\left[\frac{\left(n!\right)^2}{(2n+1)!}\right]^2\left(\frac{\epsilon}2\right)$$ Given such an $n$ we can now put together $$\begin{align}\left|\left[\frac{(2n+1)!}{\left(n!\right)^2}\right]^2I-f\left(\frac12,\frac12\right)\right|&=\left|\left[\frac{(2n+1)!}{\left(n!\right)^2}\right]^2\left(I_0+I_1+I_2\right)-f\left(\frac12,\frac12\right)\right|\\ &\le\left|\left[\frac{(2n+1)!}{\left(n!\right)^2}\right]^2\left(\left|I_1\right|+\left|I_2\right|\right)\right|\\ &<\left|\frac{\epsilon}2+\frac{\epsilon}2\right|=\epsilon\end{align}$$ Let's see... we divided by $4M$ earlier, but if $4M=0$ then the limit is trivially $f\left(1/2,1/2\right)=0$. So this limit gives us another realization of the Dirac $\delta$ function.
