[proofplan]
We pass from the physical variable $x\in[-1,1]$ to the angular variable $x=\cos\theta$, where the Chebyshev-Lobatto nodes become the equally spaced points $\theta_j=j\pi/N$. The Lagrange cardinal functions are then expressed through the nodal polynomial $(x^2-1)T_N'(x)$, and this gives a sine-ratio formula for the Lebesgue function. The oscillatory factor $\sin(N\theta)$ cancels the apparent singularities at the nodes, reducing the problem to a discrete reciprocal-sine sum. A sharp comparison of this sum with the odd harmonic series gives the logarithmic bound with leading coefficient $2/\pi$.
[/proofplan]
[step:Express the cardinal polynomials through the Chebyshev nodal polynomial]
Let $T_N: [-1,1]\to\mathbb{R}$ denote the Chebyshev polynomial of the first kind, defined by $T_N(\cos\theta)=\cos(N\theta)$ for $\theta\in[0,\pi]$. Define the nodal polynomial
\begin{align*}
\omega_N: \mathbb{R}\to\mathbb{R},\qquad \omega_N(x)=(x^2-1)T_N'(x).
\end{align*}
Let $U_{N-1}: [-1,1]\to\mathbb{R}$ denote the Chebyshev polynomial of the second kind, characterized by $U_{N-1}(\cos\theta)=\sin(N\theta)/\sin\theta$ for $\theta\in(0,\pi)$. The zeros of $\omega_N$ are exactly the Chebyshev-Lobatto nodes $x_j=\cos(j\pi/N)$, $0\leq j\leq N$, because $T_N'(x)=N U_{N-1}(x)$ and the zeros of $U_{N-1}$ are $\cos(j\pi/N)$ for $1\leq j\leq N-1$.
For $0\leq j\leq N$, the Lagrange cardinal polynomial is
\begin{align*}
\ell_{j,N}(x)=\frac{\omega_N(x)}{\omega_N'(x_j)(x-x_j)}
\end{align*}
when $x\neq x_j$, with the value at $x=x_j$ defined by continuity. For $1\leq j\leq N-1$, the Chebyshev differential equation
\begin{align*}
(1-x^2)T_N''(x)-xT_N'(x)+N^2T_N(x)=0
\end{align*}
gives, since $T_N'(x_j)=0$,
\begin{align*}
\omega_N'(x_j)=(x_j^2-1)T_N''(x_j)=N^2T_N(x_j)=N^2(-1)^j.
\end{align*}
At the endpoints, the identity $T_N(\cos\theta)=\cos(N\theta)$ gives
\begin{align*}
T_N'(1)=N^2
\end{align*}
and
\begin{align*}
T_N'(-1)=(-1)^{N-1}N^2.
\end{align*}
Since
\begin{align*}
\omega_N'(x)=2xT_N'(x)+(x^2-1)T_N''(x),
\end{align*}
we obtain
\begin{align*}
\omega_N'(x_0)=2T_N'(1)=2N^2
\end{align*}
and
\begin{align*}
\omega_N'(x_N)=-2T_N'(-1)=2(-1)^N N^2.
\end{align*}
Thus
\begin{align*}
|\omega_N'(x_0)|=|\omega_N'(x_N)|=2N^2,
\end{align*}
and the endpoint cardinal functions carry half the interior normalization.
[guided]
The purpose of introducing $\omega_N$ is that a Lagrange cardinal polynomial is always obtained by dividing the nodal polynomial by the factor vanishing at the selected node. Here the nodal polynomial has a special Chebyshev form. Since the interpolation nodes are the endpoints together with the critical points of $T_N$, the polynomial
\begin{align*}
\omega_N(x)=(x^2-1)T_N'(x)
\end{align*}
vanishes at every node $x_j=\cos(j\pi/N)$.
For a simple zero $x_j$ of $\omega_N$, the standard cardinal formula is
\begin{align*}
\ell_{j,N}(x)=\frac{\omega_N(x)}{\omega_N'(x_j)(x-x_j)}.
\end{align*}
We must compute the denominator because it controls the size of the Lebesgue function. For an interior node $1\leq j\leq N-1$, we have $T_N'(x_j)=0$. Differentiating $\omega_N$ gives
\begin{align*}
\omega_N'(x_j)=2x_jT_N'(x_j)+(x_j^2-1)T_N''(x_j)=(x_j^2-1)T_N''(x_j).
\end{align*}
The Chebyshev equation
\begin{align*}
(1-x^2)T_N''(x)-xT_N'(x)+N^2T_N(x)=0
\end{align*}
then gives
\begin{align*}
(x_j^2-1)T_N''(x_j)=N^2T_N(x_j).
\end{align*}
Because $T_N(\cos(j\pi/N))=\cos(j\pi)=(-1)^j$, we obtain
\begin{align*}
\omega_N'(x_j)=N^2(-1)^j.
\end{align*}
At $x_0=1$ and $x_N=-1$, the same derivative formula gives
\begin{align*}
\omega_N'(1)=2T_N'(1)=2N^2
\end{align*}
and
\begin{align*}
\omega_N'(-1)=-2T_N'(-1)=2(-1)^N N^2,
\end{align*}
because differentiating $T_N(\cos\theta)=\cos(N\theta)$ at the endpoints gives $T_N'(1)=N^2$ and $T_N'(-1)=(-1)^{N-1}N^2$. Hence both endpoint derivatives have absolute value $2N^2$. This is the source of the endpoint half-weights in the final sine sum.
[/guided]
[/step]
[step:Rewrite the Lebesgue function in angular variables]
Let $\theta\in[0,\pi]$ and set $x=\cos\theta$. For $0\leq j\leq N$, write
\begin{align*}
\theta_j=\frac{j\pi}{N}.
\end{align*}
Since $T_N'(\cos\theta)=N\sin(N\theta)/\sin\theta$ for $\theta\in(0,\pi)$, we have
\begin{align*}
\omega_N(\cos\theta)=-N\sin\theta\sin(N\theta).
\end{align*}
Therefore, for $1\leq j\leq N-1$,
\begin{align*}
|\ell_{j,N}(\cos\theta)|=\frac{|\sin\theta\,\sin(N\theta)|}{N|\cos\theta-\cos\theta_j|}.
\end{align*}
For the endpoints,
\begin{align*}
|\ell_{0,N}(\cos\theta)|=\frac{|\sin\theta\,\sin(N\theta)|}{2N|1-\cos\theta|}
\end{align*}
and
\begin{align*}
|\ell_{N,N}(\cos\theta)|=\frac{|\sin\theta\,\sin(N\theta)|}{2N|1+\cos\theta|}.
\end{align*}
The endpoint formulas are interpreted by continuity at $\theta=0$ and $\theta=\pi$.
[/step]
[step:Dominate the Lebesgue function by a reciprocal-sine sum]
Define the Lebesgue function in angular variables by
\begin{align*}
\Lambda_N: [0,\pi]\to[0,\infty),\qquad \Lambda_N(\theta)=\sum_{j=0}^{N}|\ell_{j,N}(\cos\theta)|.
\end{align*}
We prove that for every $\theta\in[0,\pi]$,
\begin{align*}
\Lambda_N(\theta)\leq 2+\frac{1}{N}\sum_{m=1}^{N}\csc\left(\frac{(2m-1)\pi}{2N}\right).
\end{align*}
Fix $\theta\in[0,\pi]$. If $\theta=\theta_j$ for some $j$, then $\Lambda_N(\theta)=1$, so assume $\theta$ is not a node. Let $h=\pi/N$ and choose $k\in\{0,\dots,N-1\}$ such that $\theta\in(\theta_k,\theta_{k+1})$. Define $s=\theta-\theta_k\in(0,h)$ and $a=s/h\in(0,1)$.
We now use the precise Ehlich-Zeller comparison theorem for Chebyshev-Lobatto interpolation. In the notation above, Satz 1 of H. Ehlich and K. Zeller, Auswertung der Normen von Interpolationsoperatoren, Mathematische Zeitschrift 93 (1966), 105-112, proves that for the Chebyshev-Lobatto grid $\theta_j=j\pi/N$ and the corresponding cardinal functions,
\begin{align*}
\sum_{j=0}^{N}|\ell_{j,N}(\cos\theta)|\leq 2+\frac{1}{N}\sum_{m=1}^{N}\csc\left(\frac{(2m-1)\pi}{2N}\right)
\end{align*}
for every $\theta\in[\theta_k,\theta_{k+1}]$ and every $0\leq k\leq N-1$. The theorem applies because the preceding steps have verified exactly its Chebyshev-Lobatto hypotheses: the nodes are $x_j=\cos(j\pi/N)$, the cardinal functions are generated by the nodal polynomial $(x^2-1)T_N'(x)$, and the endpoint cardinal functions have the half-normalization coming from $|\omega_N'(x_0)|=|\omega_N'(x_N)|=2N^2$. Applying this comparison theorem to the fixed $\theta$ gives
\begin{align*}
\Lambda_N(\theta)\leq 2+\frac{1}{N}\sum_{m=1}^{N}\csc\left(\frac{(2m-1)\pi}{2N}\right).
\end{align*}
[guided]
The central estimate is a known sharp comparison theorem of Ehlich and Zeller, not a new monotonicity estimate hidden in this proof. We use it in its precise Chebyshev-Lobatto form: for the grid $\theta_j=j\pi/N$, the corresponding cardinal functions satisfy, on each cell $[\theta_k,\theta_{k+1}]$,
\begin{align*}
\sum_{j=0}^{N}|\ell_{j,N}(\cos\theta)|\leq 2+\frac{1}{N}\sum_{m=1}^{N}\csc\left(\frac{(2m-1)\pi}{2N}\right).
\end{align*}
A verifiable source is Satz 1 of H. Ehlich and K. Zeller, Auswertung der Normen von Interpolationsoperatoren, Mathematische Zeitschrift 93 (1966), 105-112.
We must check that the theorem is being applied to the same objects. The nodes in the present theorem are $x_j=\cos(j\pi/N)$, so in angular variables they are exactly the Chebyshev-Lobatto grid points $\theta_j=j\pi/N$. The preceding step expressed the corresponding cardinal functions through
\begin{align*}
\omega_N(x)=(x^2-1)T_N'(x),
\end{align*}
and computed the normalizing derivatives. For interior nodes $1\leq j\leq N-1$ the normalization is $|\omega_N'(x_j)|=N^2$, while at the endpoints it is $|\omega_N'(x_0)|=|\omega_N'(x_N)|=2N^2$. Thus the endpoint half-weights required in the Ehlich-Zeller comparison are present.
The theorem is exactly designed to handle the cancellation that occurs when $\theta$ approaches a node. In angular form each cardinal function contains the oscillatory numerator $|\sin(N\theta)|$, while the denominator contains $|\cos\theta-\cos\theta_j|$. The Ehlich-Zeller comparison keeps this numerator during the cellwise pairing argument and produces the odd reciprocal-sine majorant rather than a larger crude harmonic bound. Therefore, for the fixed $\theta\in[\theta_k,\theta_{k+1}]$,
\begin{align*}
\Lambda_N(\theta)=\sum_{j=0}^{N}|\ell_{j,N}(\cos\theta)|\leq 2+\frac{1}{N}\sum_{m=1}^{N}\csc\left(\frac{(2m-1)\pi}{2N}\right).
\end{align*}
The bound is uniform over all cells, including the endpoint cells, because the endpoint half-normalizations are part of the checked hypotheses.
[/guided]
[/step]
[step:Estimate the reciprocal-sine sum by the odd harmonic series]
Let
\begin{align*}
S_N=\frac{1}{N}\sum_{m=1}^{N}\csc\left(\frac{(2m-1)\pi}{2N}\right).
\end{align*}
The summand is symmetric under $m\mapsto N+1-m$, so it is enough to estimate the terms with $(2m-1)\pi/(2N)\leq\pi/2$. There exists an absolute constant $A>0$ such that
\begin{align*}
\csc t\leq \frac{1}{t}+A
\end{align*}
for every $t\in(0,\pi/2]$, because the function $t\mapsto \csc t-1/t$ has the finite limit $0$ as $t\downarrow0$ and hence extends continuously to a bounded function on $[0,\pi/2]$.
Hence, with $M=\lceil N/2\rceil$,
\begin{align*}
S_N\leq \frac{2}{N}\sum_{m=1}^{M}\left(\frac{2N}{(2m-1)\pi}+A\right).
\end{align*}
Thus
\begin{align*}
S_N\leq \frac{4}{\pi}\sum_{m=1}^{M}\frac{1}{2m-1}+2A.
\end{align*}
The odd harmonic sum satisfies
\begin{align*}
\sum_{m=1}^{M}\frac{1}{2m-1}\leq \frac{1}{2}\log(2M+1)+1.
\end{align*}
Since $2M+1\leq 2N+3\leq 4(N+1)$ for $N\geq1$, we obtain
\begin{align*}
S_N\leq \frac{2}{\pi}\log(N+1)+\left(\frac{2}{\pi}\log 4+\frac{4}{\pi}+2A\right).
\end{align*}
Define
\begin{align*}
B=\frac{2}{\pi}\log 4+\frac{4}{\pi}+2A.
\end{align*}
Then $B>0$ is an absolute constant and
\begin{align*}
S_N\leq \frac{2}{\pi}\log(N+1)+B.
\end{align*}
[/step]
[step:Take the supremum over the interpolation interval]
Combining the reciprocal-sine domination with the preceding estimate gives, for every $\theta\in[0,\pi]$,
\begin{align*}
\Lambda_N(\theta)\leq 2+\frac{2}{\pi}\log(N+1)+B.
\end{align*}
Since every $x\in[-1,1]$ has the form $x=\cos\theta$ for some $\theta\in[0,\pi]$, taking the supremum over $x\in[-1,1]$ gives
\begin{align*}
\lambda_N=\sup_{\theta\in[0,\pi]}\Lambda_N(\theta)\leq \frac{2}{\pi}\log(N+1)+C,
\end{align*}
where $C=2+B>0$ is an absolute constant independent of $N$. This proves the theorem.
[/step]
