Since we know that we're searching for a degree 4 polynomial, we could create $p(x) = ax^4 + bx^3 + cx^2 + dx + e$. Substituting each of the five points into $p(x)$ would give us a solvable system of 5 equations and 5 unknowns. These equations would be as follows:\n", "\n", "

\n", "\n", "

\n",
"\n",
"$$a(0)^4 + b(0)^3 + c(0)^2 + d(0) + e = -1\\\\\\ a(1)^4 + b(1)^3 + c(1)^2 + d(1) + e = 0\\\\\\ a(2)^4 + b(2)^3 + c(2)^2 + d(2) + e = -11\\\\\\ a(3)^4 + b(3)^3 + c(3)^2 + d(3) + e = 2\\\\\\ a(4)^4 + b(4)^3 + c(4)^2 + d(4) + e = 99$$\n",
"\n",
"

\n",
"\n",
"\n", "\n", "However, solving this system of equations and unknowns would take quite some time.

\n", "\n", "

\n", "\n", "

\n",
"\n",
"$$\n",
"p_{i}(x) = \\frac{\\Pi_{j \\neq i} (x - x_j)}{\\Pi_{j \\neq i} (x_i - x_j)}\n",
"$$\n",
"\n",
"

\n",
"\n",
"\n", "\n", "For clarity, let's calculate $p_1(x)$ and $p_3(x)$. Recall, we had $S = \\{(0, -1), (1, 0), (2, -11), (3, 2), (4, 99)\\}$, meaning that $(x_1, y_1) = (0, -1)$ and $(x_3, y_3) = (2, -11)$.\n", "\n", "

\n", "\n", "

\n",
"$$\n",
"p_1(x) = \\frac{\\Pi_{j \\neq 1}(x - x_j)}{\\Pi_{j \\neq 1}(0 - x_j)} = \\frac{(x-1)(x-2)(x-3)(x-4)}{(0-1)(0-2)(0-3)(0-4)} = \\frac{1}{24}(x-1)(x-2)(x-3)(x-4)$$\n",
"\n",
"

\n", "$$p_3(x) = \\frac{\\Pi_{j \\neq 3}(x - x_j)}{\\Pi_{j \\neq 3}(0 - x_j)} = \\frac{(x-0)(x-1)(x-3)(x-4)}{(2-0)(2-1)(2-3)(2-4)} = \\frac{1}{4}(x)(x-1)(x-3)(x-4)\n", "$$\n", "

\n",
"\n",
"\n", "$$p_3(x) = \\frac{\\Pi_{j \\neq 3}(x - x_j)}{\\Pi_{j \\neq 3}(0 - x_j)} = \\frac{(x-0)(x-1)(x-3)(x-4)}{(2-0)(2-1)(2-3)(2-4)} = \\frac{1}{4}(x)(x-1)(x-3)(x-4)\n", "$$\n", "

\n", "\n", "The second-to-last step of the above expansions best illustrate why we've chosen to craft our sub-polynomials in this way; if we were to evaluate $p_1(0)$, the numerator and denominator would be exactly the same. If we were to instead evaluate $p_1(1), p_1(2), p_1(3)$ or $p_1(4)$, since $x-1, x-2, x-3, x-4$ are all factors of the numerator the result would be 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're almost done. We can now say that our final polynomial $p(x)$ is constructed as follows:\n", "\n", "

\n", "\n", "

\n",
"$$\n",
"p(x) = \\sum_{i = 1}^n y_i p_i(x)\n",
"$$\n",
"

\n",
"\n",
"\n", "\n", "This is where the $y$ values of each of the given points come into play. Looking at our example more closely, we have:\n", "\n", "

\n", "\n", "

\n",
"$$\n",
"p(x) = -p_1(x) + 0p_2(x) -11p_3(x) + 2p_4(x) + 99p_5(x)\n",
"$$\n",
"

\n",
"\n",
"\n", "\n", "From the way each $p_i(x)$ was constructed, $p(0) = (-1) \\cdot 1 + 0 \\cdot 0 + (-11) \\cdot 0 + 2 \\cdot 0 + 99 \\cdot 0 = -1$, and so on and so forth, as we expected. Doing the arithmetic yields the original polynomial, $p(x) = x^4 - 13x^2 + 13x - 1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = np.linspace(-1, 6)\n", "y = x**4 - 13*x**2 + 13*x - 1\n", "plt.scatter([0, 1, 2, 3, 4], [-1, 0, -11, 2, 99], color='r')\n", "plt.plot(x, y);\n", "plt.ylim(-40, 200);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementing Lagrange Interpolation in Code" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's write a (higher-order) function in Python that will allow us to take a set of points $S$ and return the polynomial that interpolates it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def sub_polynomial(S, i):\n", " def f(x):\n", " num, den = 1, 1\n", " for p in range(len(S)):\n", " if p != i:\n", " num *= (x - S[p][0])\n", " den *= (S[i][0] - S[p][0])\n", " return num / den\n", " return f\n", "\n", "def interpolate(S):\n", " def f(x):\n", " return sum([S[i][1] * sub_polynomial(S, i)(x) for i in range(len(S))])\n", " return f" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our function `sub_polynomial` creates the $p_i(x)$ functions that we saw earlier. `interpolate` puts all of this together. Let's test out `interpolate` with the same set `S` we saw before." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "S = [(0, -1), (1, 0), (2, -11), (3, 2), (4, 99)]\n", "x = np.linspace(-2, 6, 1000)\n", "y = [interpolate(S)(i) for i in x]\n", "plt.scatter([s[0] for s in S], [s[1] for s in S], color='r')\n", "plt.plot(x, y);\n", "plt.ylim(-40, 200);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's write a function that interpolates and plots any arbitrary set." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def interpolate_and_plot(S):\n", " xs = [s[0] for s in S]\n", " ys = [s[1] for s in S]\n", " \n", " f = interpolate(S)\n", " \n", " x = np.linspace(min(xs) - 1, max(xs) + 1, 1000)\n", " y = [f(i) for i in x]\n", " \n", " plt.scatter(xs, ys, color='r');\n", " plt.plot(x, y);" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "#interpolate_and_plot([(100, 2), (102, 10), (103, -3)])\n", "interpolate_and_plot([(-5, 2), (1, 5), (2, 17)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Neat! Notice what happens when we pass in $S = \\{ (1, 3), (2, 4), (3, 5) \\}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "interpolate_and_plot([(1, 3), (2, 4), (3, 5)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we passed in 3 points, so we'd expect the degree of our polynomial to be 2. However, it turns out that these three points were _colinear_, meaning they all lie on the same line. The $x^2$ term cancelled out in Lagrange Interpolation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "\n", "## Modular Arithmetic... with Polynomials?\n", "\n", "

\n", "\n", "