Multivariate Newton Interpolation

For $m,n\in \mathbb{N}$, $m\ge 1$ and a given function $f : $\mathbb{R}$^m\longrightarrow $\mathbb{R}$$, the polynomial interpolation problem (PIP) is to determine a unisolvent node set $P_{m,n}\subseteq {\mathbb{R}}^m$ of $N(m,n):=\mathrm{\mid }{P}_{m,n}\mathrm{\mid }=\left(\genfrac{}{}{0px}{}{m+n}{n}\right)$ points and the uniquely defined polynomial $Q_{m,n,f}\in {\mathrm{\Pi }}_{m,n}$ in $m$ variables of degree $\mathrm{d}\mathrm{e}\mathrm{g}(Q_{m,n,f})\le n\in \mathbb{N}$ that fits $f$ on $P_{m,n}$, i.e., $Q_{m,n,f}(p)=f(p)$, $\mathrm{\forall }p\in P_{m,n}$. For $m=1$ the solution to the PIP is well known. In higher dimensions, however, no closed framework was available. We here present a generalization of the classic Newton interpolation from one-dimensional to arbitrary-dimensional spaces. Further we formulate an algorithm, termed PIP-SOLVER, based on a multivariate divided difference scheme that computes the solution $Q_{m,n,f}$ in $\mathcal{O}(N(m,n{)}^2)$ time using $\mathcal{O}(mN(m,n))$ memory. Further, we introduce unisolvent Newton-Chebyshev nodes and show that these nodes avoid Runge’s phenomenon in the sense that arbitrary periodic Sobolev functions $f \in H^k(\Omega,$\mathbb{R}$) \subsetneq C^0(\Omega,$\mathbb{R}$)$,$\Omega =[-1,1]^m$ofregularity$k >m/2$ can be uniformly approximated, i.e., $ \lim_{n\rightarrow \infty}||,f -Q_{m,n,f} ,||_{C^0(\Omega)}= 0$. Numerical experiments demonstrate the computational performance and approximation accuracy of the PIP-SOLVER in practice. We expect the presented results to be relevant for many applications, including numerical solvers, quadrature, non-linear optimization, polynomial regression, adaptive sampling, Bayesian inference, and spectral analysis.