Return to Answer

By the Euler–Maclaurin formula (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula), $$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\ =-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$ for each $f\in\mathcal F$.

However, one cannot get a constant factor $O(n^{-4})$ good for all $f\in\mathcal F$ at once. Specifically, for $n\ge2$, $$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n) =\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2}$$ as $n\to\infty$, where $$f_n(x):=\frac1{\pi^2}\, \Big(\sin (\pi x)-\frac{\sin (\pi (2 n+1) x)}{(2 n+1)^2}\Big),$$ so that $\int_0^1(f_n'')^2=1$ and $f_n(0)=0=f_n(1)$.

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula $$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $f\in\mathcal F$.)

By the Euler–Maclaurin formula (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula), $$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\ =-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$ for each $f\in\mathcal F$.

However, one cannot get a constant factor $O(n^{-4})$ good for all $f\in\mathcal F$ at once. Specifically, for $n\ge2$, $$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n) =\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2} \tag{1}\label{1}$$ as $n\to\infty$, where $$f_n(x):=\frac1{\pi^2}\, \Big(\sin (\pi x)-\frac{\sin (\pi (2 n+1) x)}{(2 n+1)^2}\Big),$$ so that $\int_0^1(f_n'')^2=1$ and $f_n(0)=0=f_n(1)$.

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula $$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $f\in\mathcal F$.)

On the other hand, even without the conditions $f(1)=f'(1)=0$ on $f\in\mathcal F$ (added later by the OP), one can see that $$\sup_{f\in\mathcal F}d_n(f)\le \frac1{72 n^5}+\frac1{12 n^2}\sim\frac1{12 n^2} \tag{2}\label{2}$$ as $n\to\infty$, so that the lower bound on $\sup_{f\in\mathcal F}d_n(f)$ in \eqref{1} is sharp up to a universal positive real constant factor.

To get \eqref{2}, use the Euler–Maclaurin formula again, but this time with $p=2$, which together with the condition $\int_0^1(f'')^2\le1$ yields $$d_n(f)\le-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+\frac1{12 n^2} \tag{3}\label{3}$$ for all $f\in\mathcal F$. Using the condition $\int_0^1(f'')^2\le1$ again, now together with the condition $f'(0)=0$, we get $|f'(1)|\le\int_0^1|f''|\le1$, whence $$-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1) \le-\frac1{2n}\,f^2(1)+\frac1{6n^3}\,|f(1)|\le\frac1{72 n^5},$$ so that \eqref{2} follows from \eqref{3}.

By the Euler–Maclaurin formula (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula), $$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\ =-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$ for each $f\in\mathcal F$.

However, one cannot get a constant factor $O(n^{-4})$ good for all $f\in\mathcal F$ at once. Specifically,
  $$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n) =\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2}$$ as $n\to\infty$, where $$f_n(x):=\frac1{\pi^2}\, \Big(\sin (\pi x)-\frac{\sin (\pi (2 n+1) x)}{(2 n+1)^2}\Big),$$ so that $\int_0^1(f_n'')^2=1$ and $f_n(0)=0=f_n(1)$.

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula $$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $\mathcal F$.)

By the Euler–Maclaurin formula (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula), $$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\ =-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$ for each $f\in\mathcal F$.

However, one cannot get a constant factor $O(n^{-4})$ good for all $f\in\mathcal F$ at once. Specifically, for $n\ge2$, $$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n) =\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2}$$ as $n\to\infty$, where $$f_n(x):=\frac1{\pi^2}\, \Big(\sin (\pi x)-\frac{\sin (\pi (2 n+1) x)}{(2 n+1)^2}\Big),$$ so that $\int_0^1(f_n'')^2=1$ and $f_n(0)=0=f_n(1)$.

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula $$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $f\in\mathcal F$.)

By the Euler–Maclaurin formula (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula), $$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\ =-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$ for each $f\in\mathcal F$.

However, one cannot get a constant factor $O(n^{-4})$ good for all $f$ at once. Specifically,
$$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n) =\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2}$$ as $n\to\infty$, where $$f_n(x):=\frac1{\pi^2}\, \Big(\sin (\pi x)-\frac{\sin (\pi (2 n+1) x)}{(2 n+1)^2}\Big).$$

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula $$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $\mathcal F$.)

By the Euler–Maclaurin formula (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula), $$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\ =-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$ for each $f\in\mathcal F$.

However, one cannot get a constant factor $O(n^{-4})$ good for all $f\in\mathcal F$ at once. Specifically,
$$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n) =\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2}$$ as $n\to\infty$, where $$f_n(x):=\frac1{\pi^2}\, \Big(\sin (\pi x)-\frac{\sin (\pi (2 n+1) x)}{(2 n+1)^2}\Big),$$ so that $\int_0^1(f_n'')^2=1$ and $f_n(0)=0=f_n(1)$.

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula $$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $\mathcal F$.)