Math 441

Interpolation: The Algorithm

Newton’s Polynomial

  • Data:
    • \(x_i\), \(i = 0, 1, \dots, n\)
    • \(y_i\), \(i = 0, 1, \dots, n\)
  • The polynomial \[p(x) = a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \cdots + a_n(x-x_0)(x-x_1)\cdots(x-x_{n-1})\]
    • \(a_i\), \(i = 0, 1, \dots, n\)
    • \(x_i\), \(i = 0, 1, \dots, {\color{green}n-1}\)
    • Evaluated using Horner’s algorithm
  • The coefficients \(a_i\), \(i = 0, 1, \dots, n\)
    • Generated from \(x_i\)’s and \(y_i\)’s …
    • using divided differences

Horner’s Algorithm

\[p(z_0,z_1,\dots,z_{n-1}) = a_0 + a_1z_0 + a_2z_0z_1 + \cdots + a_nz_0z_1\cdots z_{n-1}\]

  • Start with \(a_n\).
  • Multiply by \(z_{n-1}\).
  • Add \(a_{n-1}\).
  • Multiply by \(z_{n}\).
  • Add \(a_{n}\).
  • \(\vdots\)
  • Multiply by \(z_{k}\).
  • Add \(a_{k}\).
  • \(\vdots\)
  • Multiply by \(z_{0}\) and add \(a_{0}\).
  • Start with \(a_n\).
  • For \(k\) from \(n-1\) down to 0:
    • Multiply by \(z_{k}\).
    • Add \(a_{k}\).

Horner in Julia

  • Start with \(a_n\).
  • For \(k\) from \(n-1\) down to 0:
    • Multiply by \(z_{k}\).
    • Add \(a_{k}\).
  • Julia indexing starts at 1, not 0.
    • a[1] instead of \(a_0\).
    • a[n+1] instead of \(a_n\).
  • Multiply - Add is very common
    • many CPUs do it with a single instruction
    • muladd() function
  • Inputs:
    • a :: Array of length \(n+1\)
    • z :: Array of length \(n\)
p = a[n+1]
for k in n:-1:1
    p = muladd(p, z[k], a[k])
end

Divided Differences Example

Use the Newton’s basis to find a polynomial \(p\) of degree 2 such that \(p(-1) = 11\), \(p(1) = -1\), and \(p(2) = 2\).

𝑥 𝑖 𝑓[𝑥 𝑖 ] 𝑓[𝑥 𝑖 ,𝑥 𝑗 ] 𝑓[𝑓 𝑖 ,𝑓 𝑗 ,𝑥 𝑘 ] −1 11 −6 1 −1 3 3 2 2

\(p(x) = 11 - 6(x+1) + 3(x+1)(x-1)\)

Divided Differences

𝑥 0 𝑥 1 𝑥 2 𝑥 3 𝑥 4 𝑦 0 𝑦 1 𝑦 2 𝑦 3 𝑦 4 𝑓[𝑥 0 ] 𝑎 0 𝑓[𝑥 1 ] 𝑓[𝑥 1 ] 𝑓[𝑥 2 ] 𝑓[𝑥 2 ] 𝑓[𝑥 3 ] 𝑓[𝑥 3 ] 𝑓[𝑥 4 ] 𝑓[𝑥 4 ] 𝑓[𝑥 0 𝑥 1 ]= 𝑓[𝑥 1 ]−𝑓[𝑥 0 ] 𝑥 1 −𝑥 0 𝑓[𝑥 0 𝑥 1 ]= 𝑓[𝑥 1 ]−𝑓[𝑥 0 ] 𝑥 1 −𝑥 0 𝑓[𝑥 0 𝑥 1 ] 𝑎 1 𝑓[𝑥 1 𝑥 2 ]= 𝑓[𝑥 2 ]−𝑓[𝑥 1 ] 𝑥 2 −𝑥 1 𝑓[𝑥 1 𝑥 2 ]= 𝑓[𝑥 2 ]−𝑓[𝑥 1 ] 𝑥 2 −𝑥 1 𝑓[𝑥 1 𝑥 2 ] 𝑓[𝑥 2 𝑥 3 ]= 𝑓[𝑥 3 ]−𝑓[𝑥 2 ] 𝑥 3 −𝑥 2 𝑓[𝑥 2 𝑥 3 ]= 𝑓[𝑥 3 ]−𝑓[𝑥 2 ] 𝑥 3 −𝑥 2 𝑓[𝑥 2 𝑥 3 ] 𝑓[𝑥 3 𝑥 4 ]= 𝑓[𝑥 4 ]−𝑓[𝑥 3 ] 𝑥 4 −𝑥 3 𝑓[𝑥 3 𝑥 4 ]= 𝑓[𝑥 4 ]−𝑓[𝑥 3 ] 𝑥 4 −𝑥 3 𝑓[𝑥 3 𝑥 4 ] 𝑓[𝑥 0 𝑥 1 𝑥 2 ]= 𝑓[𝑥 1 𝑥 2 ]−𝑓[𝑥 0 𝑥 1 ] 𝑥 2 −𝑥 0 𝑓[𝑥 0 𝑥 1 𝑥 2 ] 𝑎 2 𝑓[𝑥 1 𝑥 2 𝑥 3 ]= 𝑓[𝑥 2 𝑥 3 ]−𝑓[𝑥 1 𝑥 2 ] 𝑥 3 −𝑥 1 𝑓[𝑥 1 𝑥 2 𝑥 3 ] 𝑓[𝑥 2 𝑥 3 𝑥 4 ]= 𝑓[𝑥 3 𝑥 4 ]−𝑓[𝑥 2 𝑥 3 ] 𝑥 4 −𝑥 2 𝑓[𝑥 2 𝑥 3 𝑥 4 ] 𝑓[𝑥 0 𝑥 1 𝑥 2 𝑥 3 ]= 𝑓[𝑥 1 𝑥 2 𝑥 3 ]−𝑓[𝑥 0 𝑥 1 𝑥 2 ] 𝑥 3 −𝑥 0 𝑓[𝑥 0 𝑥 1 𝑥 2 𝑥 3 ] 𝑎 3 𝑓[𝑥 1 𝑥 2 𝑥 3 𝑥 4 ]= 𝑓[𝑥 2 𝑥 3 𝑥 4 ]−𝑓[𝑥 1 𝑥 2 𝑥 3 ] 𝑥 4 −𝑥 1 𝑓[𝑥 1 𝑥 2 𝑥 3 𝑥 4 ] 𝑓[𝑥 0 𝑥 1 𝑥 2 𝑥 3 𝑥 4 ]= 𝑓[𝑥 1 𝑥 2 𝑥 3 𝑥 4 ]−𝑓[𝑥 0 𝑥 1 𝑥 2 𝑥 3 ] 𝑥 4 −𝑥 0 𝑓[𝑥 0 𝑥 1 𝑥 2 𝑥 3 𝑥 4 ] 𝑎 4

The Algorithm

  • The first divided differences are the \(y_i\)’s.
  • To find the second divided differences:
    • Calculate \(f[x_kx_{k-1}]=\frac{f[x_k] - f[x_{k-1}]}{x_k - x_{k-{\color{green}1}}} = \frac{y_k - y_{k-1}}{x_k - x_{k-1}}\),
    • starting with the highest \(k\),
    • replacing \(y_k\) by the result, so that \(y_k = f[x_kx_{k-1}]\).
    • then the previous \(k\), …, down to the second \(k\).
  • To find the third divided differences:
    • Calculate \(f[k_kx_{k-1}x_{k-2}]=\frac{f[x_kx_{k-1}] - f[x_{k-1}x_{k-2}]}{x_k - x_{k-{\color{green}2}}} = \frac{y_k - y_{k-1}}{x_k - x_{x-2}}\)
    • starting with the highest \(k\),
    • replacing \(y_k\) by the result,
    • then the previous \(k\), …, down to the third \(k\).
  • And so on…

The \(i\)-th divided differences

  • The first \(i-1\) entries in \(y_i\) already are the final \(a\)’s
  • The rest of the entries contain the \((i-1)\)-st divided differences:
    • \(y_k = f[x_kx_{k-1}\cdots x_{k-i+1}]\) for \((i-1)\)-st, \(i\)-th, \((i+1)\)-st, … entries
  • Starting with the highest \(k\), and going down to the \(i\)-th entry, calculate the new \(y_k\): \[y_k = \frac{y_k - y_{k-1}}{x_k - x_{k - i + 1}}\]

Inputs:

  • x :: Array
  • y :: Array
  • Both have length \(m\).

Code:

for i in 2:m        # Start with the second diffs
    for k in m:-1:i # From m down to i
        y[k] = (y[k] - y[k-1])/(x[k] - x[k-i+1])
    end
end