Return to Answer

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Also also also, you could construct a linear system where the unknowns are the coefficients of your polynomial.

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Also also also, you could construct a linear system where the unknowns are the coefficients of your polynomial. This way, you could respect a prescribed monomial support (like if you want a circle).

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Also also also, you could construct a linear system where the unknowns are the coefficients of your polynomials.

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Also also also, you could construct a linear system where the unknowns are the coefficients of your polynomial.

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Also also also, you could construct a linear system where the unknowns are the coefficients of your polynomials.