[proofplan] The interpolation operator is characterized by a universal property: its value on a function is the unique polynomial of degree at most $n$ matching that function at the nodes $x_0,\dots,x_n$. We take the interpolants of $f$ and $g$, form the same linear combination of those polynomial functions, and check that this new polynomial has the node values prescribed by $\alpha f+\beta g$. Uniqueness of the interpolation polynomial then forces this polynomial to equal the interpolant of $\alpha f+\beta g$. [/proofplan]

[step:Form the candidate polynomial by taking the same linear combination of interpolants] Define polynomial functions $p,q:\mathbb{R}\to\mathbb{R}$ by \begin{align*} p:=I_{x_0,\dots,x_n}f \end{align*} and \begin{align*} q:=I_{x_0,\dots,x_n}g. \end{align*} By the definition of $I_{x_0,\dots,x_n}$, we have $p,q \in \mathcal{P}_n(\mathbb{R})$ and \begin{align*} p(x_i)=f(x_i) \end{align*} and \begin{align*} q(x_i)=g(x_i) \end{align*} for every $i \in \{0,1,\dots,n\}$. Define the polynomial function $r:\mathbb{R}\to\mathbb{R}$ by \begin{align*} r(t):=\alpha p(t)+\beta q(t). \end{align*} Since $\mathcal{P}_n(\mathbb{R})$ is a real [vector space](/page/Vector%20Space), $r=\alpha p+\beta q$ belongs to $\mathcal{P}_n(\mathbb{R})$. [/step]

[step:Verify that the candidate has the interpolation data of $\alpha f+\beta g$]For each $i \in \{0,1,\dots,n\}$, evaluating the definition of $r$ at the node $x_i$ gives \begin{align*} r(x_i)=\alpha p(x_i)+\beta q(x_i). \end{align*} Using the defining interpolation identities for $p$ and $q$, this becomes \begin{align*} r(x_i)=\alpha f(x_i)+\beta g(x_i). \end{align*} By the pointwise vector space structure on the set of functions $A\to\mathbb{R}$, the function $\alpha f+\beta g:A\to\mathbb{R}$ satisfies \begin{align*} (\alpha f+\beta g)(x_i)=\alpha f(x_i)+\beta g(x_i). \end{align*} Hence \begin{align*} r(x_i)=(\alpha f+\beta g)(x_i) \end{align*} for every $i \in \{0,1,\dots,n\}$.[/step]

[guided]The goal is to identify the interpolating polynomial of the function $\alpha f+\beta g:A\to\mathbb{R}$. Since interpolation is determined by values at the nodes, we test whether the polynomial function $r=\alpha p+\beta q$ has exactly those node values. For every $i \in \{0,1,\dots,n\}$, the function $r:\mathbb{R}\to\mathbb{R}$ was defined by $r(t)=\alpha p(t)+\beta q(t)$, so substitution of $t=x_i$ gives \begin{align*} r(x_i)=\alpha p(x_i)+\beta q(x_i). \end{align*} The functions $p$ and $q$ are the interpolants of $f$ and $g$, respectively. Therefore their defining interpolation properties give \begin{align*} p(x_i)=f(x_i) \end{align*} and \begin{align*} q(x_i)=g(x_i). \end{align*} Substituting these two equalities into the formula for $r(x_i)$ gives \begin{align*} r(x_i)=\alpha f(x_i)+\beta g(x_i). \end{align*} Finally, the linear combination $\alpha f+\beta g:A\to\mathbb{R}$ is defined pointwise, meaning that for every $a\in A$, \begin{align*} (\alpha f+\beta g)(a)=\alpha f(a)+\beta g(a). \end{align*} Because each node $x_i$ belongs to $A$, we may set $a=x_i$ and obtain \begin{align*} (\alpha f+\beta g)(x_i)=\alpha f(x_i)+\beta g(x_i). \end{align*} Thus \begin{align*} r(x_i)=(\alpha f+\beta g)(x_i) \end{align*} for every node $x_i$.[/guided]

[step:Use uniqueness of interpolation to identify the candidate with the desired interpolant] We have shown that $r\in\mathcal{P}_n(\mathbb{R})$ and that $r$ agrees with $\alpha f+\beta g$ at all nodes $x_0,\dots,x_n$. By the uniqueness part of the interpolation theorem for pairwise distinct nodes, equivalently [citetheorem:8328], the polynomial function in $\mathcal{P}_n(\mathbb{R})$ with these prescribed values is unique. Therefore \begin{align*} r=I_{x_0,\dots,x_n}(\alpha f+\beta g). \end{align*} Since $r=\alpha p+\beta q$, $p=I_{x_0,\dots,x_n}f$, and $q=I_{x_0,\dots,x_n}g$, we conclude that \begin{align*} I_{x_0,\dots,x_n}(\alpha f+\beta g)=\alpha I_{x_0,\dots,x_n}f+\beta I_{x_0,\dots,x_n}g. \end{align*} This is the desired linearity identity in $\mathcal{P}_n(\mathbb{R})$. [/step]