<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/lagrangian-polynomial-interpolation-with-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/lagrangian-polynomial-interpolation-r/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>Part 2 of 4 in the series <a rel="nofollow" href="http://www.aaronschlegel.com/series/numerical-analysis/" title="Numerical Analysis" target="_blank">Numerical Analysis</a></div>
<p>Polynomial interpolation is the method of determining a polynomial that fits a set of given points. There are several approaches to polynomial interpolation, of which one of the most well known is the Lagrangian method. This post will introduce the Lagrangian method of interpolating polynomials and how to perform the procedure in R.</p>
<h6>Lagrangian Polynomial Interpolation</h6>
<p>The Lagrangian method of polynomial interpolation uses Lagrangian polynomials to fit a polynomial to a given set of data points. The Lagrange interpolating polynomial is given by the following theorem:</p>
<p>For a set of data points <img alt="(x_0, y_0), (x_1, y_1), \cdots, (x_n, y_n)" title="(x_0, y_0), (x_1, y_1), \cdots, (x_n, y_n)" src="https://s0.wp.com/latex.php?latex=%28x_0%2C+y_0%29%2C+%28x_1%2C+y_1%29%2C+%5Ccdots%2C+%28x_n%2C+y_n%29&bg=ffffff&fg=000&s=0" class="latex"> with no duplicate <img alt="x" title="x" src="https://s0.wp.com/latex.php?latex=x&bg=ffffff&fg=000&s=0" class="latex"> and there exists a function <img alt="f" title="f" src="https://s0.wp.com/latex.php?latex=f&bg=ffffff&fg=000&s=0" class="latex"> which evaluates to these points, then there is a unique polynomial <img alt="P(x)" title="P(x)" src="https://s0.wp.com/latex.php?latex=P%28x%29&bg=ffffff&fg=000&s=0" class="latex"> with degree <img alt="\leq n" title="\leq n" src="https://s0.wp.com/latex.php?latex=%5Cleq+n&bg=ffffff&fg=000&s=0" class="latex"> also exists. The polynomial is given by:</p>
<img alt=" P(x) = f(x_o)L_{n,0}(x) + \cdots + f(x_n)L_{n,n}(x) = \sum^n_{k=0} f(x_k) L_{n,k}(x) " title=" P(x) = f(x_o)L_{n,0}(x) + \cdots + f(x_n)L_{n,n}(x) = \sum^n_{k=0} f(x_k) L_{n,k}(x) " src="https://s0.wp.com/latex.php?latex=+P%28x%29+%3D+f%28x_o%29L_%7Bn%2C0%7D%28x%29+%2B+%5Ccdots+%2B+f%28x_n%29L_%7Bn%2Cn%7D%28x%29+%3D+%5Csum%5En_%7Bk%3D0%7D+f%28x_k%29+L_%7Bn%2Ck%7D%28x%29+&bg=ffffff&fg=000&s=0" class="latex">
<p>Where each <img alt="k" title="k" src="https://s0.wp.com/latex.php?latex=k&bg=ffffff&fg=000&s=0" class="latex"> in <img alt="k = 0, 1, \cdots, n" title="k = 0, 1, \cdots, n" src="https://s0.wp.com/latex.php?latex=k+%3D+0%2C+1%2C+%5Ccdots%2C+n&bg=ffffff&fg=000&s=0" class="latex"> is:</p>
<img alt=" L_{n,k} = \frac{(x - x_0)(x - x_1) \cdots (x - x_{k-1})(x - x_{k+1}) \cdots (x - x_n)}{(x_k - x_0)(x_k - x_1) \cdots (x_k - x_{k-1})(x_k - x_{k+1}) \cdots (x_k - x_n)} = \prod_{\stackrel{i \neq k}{i=0}}^{n} \frac{(x - x_i)}{(x_k - x_i)}" title=" L_{n,k} = \frac{(x - x_0)(x - x_1) \cdots (x - x_{k-1})(x - x_{k+1}) \cdots (x - x_n)}{(x_k - x_0)(x_k - x_1) \cdots (x_k - x_{k-1})(x_k - x_{k+1}) \cdots (x_k - x_n)} = \prod_{\stackrel{i \neq k}{i=0}}^{n} \frac{(x - x_i)}{(x_k - x_i)}" src="https://s0.wp.com/latex.php?latex=+L_%7Bn%2Ck%7D+%3D+%5Cfrac%7B%28x+-+x_0%29%28x+-+x_1%29+%5Ccdots+%28x+-+x_%7Bk-1%7D%29%28x+-+x_%7Bk%2B1%7D%29+%5Ccdots+%28x+-+x_n%29%7D%7B%28x_k+-+x_0%29%28x_k+-+x_1%29+%5Ccdots+%28x_k+-+x_%7Bk-1%7D%29%28x_k+-+x_%7Bk%2B1%7D%29+%5Ccdots+%28x_k+-+x_n%29%7D+%3D+%5Cprod_%7B%5Cstackrel%7Bi+%5Cneq+k%7D%7Bi%3D0%7D%7D%5E%7Bn%7D+%5Cfrac%7B%28x+-+x_i%29%7D%7B%28x_k+-+x_i%29%7D&bg=ffffff&fg=000&s=0" class="latex">
<h6>Lagrangian Polynomial Interpolation Example</h6>
<p>Consider the following set of data points:</p>
<table>
<thead>
<tr>
<th>x</th>
<th>y</th>
</tr>
</thead>
<tbody>
<tr>
<td>0</td>
<td>7</td>
</tr>
<tr>
<td>2</td>
<td>11</td>
</tr>
<tr>
<td>3</td>
<td>28</td>
</tr>
<tr>
<td>4</td>
<td>63</td>
</tr>
</tbody>
</table>
<p>The graph of the points looks like:</p>
<pre>
x &lt;- c(0, 2, 3, 4)
y &lt;- c(7, 11, 28, 63)

dat &lt;- data.frame(cbind(x, y))

ggplot(dat, aes(x=x, y=y)) + 
  geom_point(size=5, col='blue')
</pre>
<p><a rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22552" target="_blank"><img sizes="(max-width: 1060px) 100vw, 1060px" src="https://i1.wp.com/www.aaronschlegel.com/wp-content/uploads/plotted_points.png?resize=1060%2C757" alt="plotted x-y points" class="wp-image-22552" data-recalc-dims="1"></a></p>
<p>We wish to find a polynomial that passes through these points. Note that much of the intermediate algebra will be skipped for the sake of brevity.</p>
<p>Start with <img alt="L_1(x)" title="L_1(x)" src="https://s0.wp.com/latex.php?latex=L_1%28x%29&bg=ffffff&fg=000&s=0" class="latex">:</p>
<p><img alt=" \large{L_1(x) = \frac{(x-2)(x-3)(x-4)}{(0-2)(0-3)(0-4)} = -\frac{1}{24}(x-2)(x-3)(x-4)} " title=" \large{L_1(x) = \frac{(x-2)(x-3)(x-4)}{(0-2)(0-3)(0-4)} = -\frac{1}{24}(x-2)(x-3)(x-4)} " src="https://s0.wp.com/latex.php?latex=+%5Clarge%7BL_1%28x%29+%3D+%5Cfrac%7B%28x-2%29%28x-3%29%28x-4%29%7D%7B%280-2%29%280-3%29%280-4%29%7D+%3D+-%5Cfrac%7B1%7D%7B24%7D%28x-2%29%28x-3%29%28x-4%29%7D+&bg=ffffff&fg=000&s=0" class="latex"><br>
<img alt="L_2(x)" title="L_2(x)" src="https://s0.wp.com/latex.php?latex=L_2%28x%29&bg=ffffff&fg=000&s=0" class="latex"> is then:<br>
<img alt=" \large{L_2(x) = \frac{(x-0)(x-3)(x-4)}{(2-0)(2-3)(2-4)} = \frac{1}{4}x(x-3)(x-4)} " title=" \large{L_2(x) = \frac{(x-0)(x-3)(x-4)}{(2-0)(2-3)(2-4)} = \frac{1}{4}x(x-3)(x-4)} " src="https://s0.wp.com/latex.php?latex=+%5Clarge%7BL_2%28x%29+%3D+%5Cfrac%7B%28x-0%29%28x-3%29%28x-4%29%7D%7B%282-0%29%282-3%29%282-4%29%7D+%3D+%5Cfrac%7B1%7D%7B4%7Dx%28x-3%29%28x-4%29%7D+&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>Then finally, <img alt="L_3(x)" title="L_3(x)" src="https://s0.wp.com/latex.php?latex=L_3%28x%29&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="L_4(x)" title="L_4(x)" src="https://s0.wp.com/latex.php?latex=L_4%28x%29&bg=ffffff&fg=000&s=0" class="latex">:</p>
<p><img alt=" \large{L_3(x) = \frac{(x-0)(x-2)(x-4)}{(3-0)(3-2)(3-4)} = -\frac{1}{3}x(x-2)(x-4)} " title=" \large{L_3(x) = \frac{(x-0)(x-2)(x-4)}{(3-0)(3-2)(3-4)} = -\frac{1}{3}x(x-2)(x-4)} " src="https://s0.wp.com/latex.php?latex=+%5Clarge%7BL_3%28x%29+%3D+%5Cfrac%7B%28x-0%29%28x-2%29%28x-4%29%7D%7B%283-0%29%283-2%29%283-4%29%7D+%3D+-%5Cfrac%7B1%7D%7B3%7Dx%28x-2%29%28x-4%29%7D+&bg=ffffff&fg=000&s=0" class="latex"><br>
<img alt=" \large{L_4(x) = \frac{(x-0)(x-2)(x-3)}{(4-0)(4-2)(4-3)} = \frac{1}{8}x(x-2)(x-3)} " title=" \large{L_4(x) = \frac{(x-0)(x-2)(x-3)}{(4-0)(4-2)(4-3)} = \frac{1}{8}x(x-2)(x-3)} " src="https://s0.wp.com/latex.php?latex=+%5Clarge%7BL_4%28x%29+%3D+%5Cfrac%7B%28x-0%29%28x-2%29%28x-3%29%7D%7B%284-0%29%284-2%29%284-3%29%7D+%3D+%5Cfrac%7B1%7D%7B8%7Dx%28x-2%29%28x-3%29%7D+&bg=ffffff&fg=000&s=0" class="latex"></p>
<p>These polynomials are then multiplied by they corresponding <img alt="y" title="y" src="https://s0.wp.com/latex.php?latex=y&bg=ffffff&fg=000&s=0" class="latex"> value and added together, resulting in the following polynomial:</p>
<img alt=" 7\left(-\frac{1}{24}(x-2)(x-3)(x-4)\right) + 11\left(\frac{1}{4}x(x-3)(x-4)\right) - 28\left(-\frac{1}{3}x(x-2)(x-4)\right) + 63\left(\frac{1}{8}x(x-2)(x-3)\right) " title=" 7\left(-\frac{1}{24}(x-2)(x-3)(x-4)\right) + 11\left(\frac{1}{4}x(x-3)(x-4)\right) - 28\left(-\frac{1}{3}x(x-2)(x-4)\right) + 63\left(\frac{1}{8}x(x-2)(x-3)\right) " src="https://s0.wp.com/latex.php?latex=+7%5Cleft%28-%5Cfrac%7B1%7D%7B24%7D%28x-2%29%28x-3%29%28x-4%29%5Cright%29+%2B+11%5Cleft%28%5Cfrac%7B1%7D%7B4%7Dx%28x-3%29%28x-4%29%5Cright%29+-+28%5Cleft%28-%5Cfrac%7B1%7D%7B3%7Dx%28x-2%29%28x-4%29%5Cright%29+%2B+63%5Cleft%28%5Cfrac%7B1%7D%7B8%7Dx%28x-2%29%28x-3%29%5Cright%29+&bg=ffffff&fg=000&s=0" class="latex">
<p>Simplifying this equation yields:</p>
<img alt=" \Large{x^3 - 2x + 7} " title=" \Large{x^3 - 2x + 7} " src="https://s0.wp.com/latex.php?latex=+%5CLarge%7Bx%5E3+-+2x+%2B+7%7D+&bg=ffffff&fg=000&s=0" class="latex">
<h6>Polynomial Interpolation in R</h6>
<p>We can implement a function in R to perform Lagrangian interpolation. The <a rel="nofollow" href="https://cran.r-project.org/web/packages/rSymPy/index.html" target="_blank">rSymPy package</a>, which provides an R interface for the SymPy CAS (computer algebra system), will be used as we are interested in the symbolic form of the resulting polynomial.</p>
<pre>
library(rSymPy)
</pre>
<p>Our implementation of the Lagrangian polynomial interpolation method, <code>lagrange.poly()</code>, proceeds as follows:</p>
<pre>
lagrange.poly &lt;- function(x, y) {

  l &lt;- list() # List to store Lagrangian polynomials L_{1,2,3,4}
  k &lt;- 1
  
  for (i in x) {
    # Set the numerator and denominator of the Lagrangian polynomials to 1 and build them up
    num &lt;- 1
    denom &lt;- 1

    # Remove the current x value from the iterated list
    p &lt;- x[! x %in% i]

    # For the remaining points, construct the Lagrangian polynomial by successively 
    # appending each x value
    for (j in p) {
      num &lt;- paste(num, &quot;*&quot;, &quot;(&quot;, 'x', &quot; - &quot;, as.character(j), &quot;)&quot;, sep = &quot;, collapse = &quot;)
      denom &lt;- paste(denom, &quot;*&quot;, &quot;(&quot;, as.character(i),&quot; - &quot;, as.character(j), &quot;)&quot;, sep = &quot;, collapse = &quot;)
    }
  
    # Set each Lagrangian polynomial in rSymPy to simplify later.
    l[k] &lt;- paste(&quot;(&quot;, num, &quot;)&quot;, &quot;/&quot;, &quot;(&quot;, denom, &quot;)&quot;, sep = &quot;, collapse = &quot;)
    k &lt;- k + 1
  }
  
  # Similar to before, we construct the final Lagrangian polynomial by successively building 
  # up the equation by iterating through the polynomials L_{1,2,3,4} and the y values 
  # corresponding to the x values.
  eq &lt;- 0
  
  for (i in 1:length(y)) {
    eq &lt;- paste(eq, '+', as.character(y[i]), &quot;*&quot;, l[[i]], sep = &quot;, collapse = &quot;)
  }
  
  # Define x variable for rSymPy to simplify
  x &lt;- Var('x')
  
  # Simplify the result with rSymPy and return the polynomial
  return(sympy(paste(&quot;simplify(&quot;, eq, &quot;)&quot;)))
}
</pre>
<p>With the <img alt="x" title="x" src="https://s0.wp.com/latex.php?latex=x&bg=ffffff&fg=000&s=0" class="latex"> and <img alt="y" title="y" src="https://s0.wp.com/latex.php?latex=y&bg=ffffff&fg=000&s=0" class="latex"> values used previously, perform Lagrangian polynomial interpolation to see if our function returns the same polynomial we found earlier.</p>
<pre>
lagrange.poly(x, y)
## [1] &quot;7 - 2*x + x**3&quot;
</pre>
<p>The polynomial <img alt="x^3 - 2x + 7" title="x^3 - 2x + 7" src="https://s0.wp.com/latex.php?latex=x%5E3+-+2x+%2B+7&bg=ffffff&fg=000&s=0" class="latex"> resulting from our <code>lagrange.poly()</code> function matches what we found previously.</p>
<p>Lagrangian polynomial interpolation can also be performed in R using the <code>poly.calc()</code> function from the <a rel="nofollow" href="https://cran.r-project.org/web/packages/polynom/index.html" target="_blank">polynom package</a>.</p>
<pre>
library(polynom)
</pre>
<p>Use the <code>poly.calc()</code> function to interpolate the data points using the Lagrangian method.</p>
<pre>
poly.calc(x, y)
## 7 - 2*x + x^3
</pre>
<p>We can plot this polynomial against the original data points to confirm visually that the polynomial does indeed pass through the points.</p>
<pre>
f &lt;- function(x) {
  return(x^3 - 2 * x + 7)
}

ggplot(dat, aes(x=x, y=y)) + 
  geom_point(size=5, col='blue') + 
  stat_function(fun = f, size=1.25, alpha=0.4)
</pre>
<p><a rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22553" target="_blank"><img sizes="(max-width: 1060px) 100vw, 1060px" src="https://i0.wp.com/www.aaronschlegel.com/wp-content/uploads/plotted_points_w_func.png?resize=1060%2C757" alt="plotted points with polynomial resulting from polynomial interpolation method" class="wp-image-22553" data-recalc-dims="1"></a></p>
<h6>References</h6>
<p>Burden, R. L., &amp; Faires, J. D. (2011). Numerical analysis (9th ed.). Boston, MA: Brooks/Cole, Cengage Learning.</p>
<p>Cheney, E. W., &amp; Kincaid, D. (2013). Numerical mathematics and computing (6th ed.). Boston, MA: Brooks/Cole, Cengage Learning.</p>
<p>The post <a rel="nofollow" href="http://www.aaronschlegel.com/lagrangian-polynomial-interpolation-r/" target="_blank">Lagrangian Polynomial Interpolation with R</a> appeared first on <a rel="nofollow" href="http://www.aaronschlegel.com/" target="_blank">Aaron Schlegel</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.aaronschlegel.com/lagrangian-polynomial-interpolation-r/"> R – Aaron Schlegel</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/ptvNBWLsXSU" height="1"></a></div>