Hermite interpolation by polynomials

Directly to the input form

• By interpolation one means a process to replace a function, given by a finite set of data, by a "simple" function, given as a (in the range of interest) everywhere evaluable formula, with the aim of evaluation or otherwise mathematical manipulation (say integration, differentiation etc). Here we have a special case: the "simple" function is chosen as a polynomial of given maximum degree and the dimension of the independent variable is one.
• We have data x_i, y_i, y_i⁽¹⁾, i=0,...,n, (with the meaning of y_i = f(x_i), y⁽¹⁾_i = f'(x_i), f not necessarily known), and want to compute a polynomial p of maximum degree 2n+1 such that
  p(x_i) = y_i, i=0,...,n .
  p'(x_i) = y_i⁽¹⁾, i=0,...,n .
  The unique solvability of this problem requires only that x_i, i=0,...,n, are pairwise different.
• This problem can be solved in a variety of ways, depending on the choice of the basis of the linear space of polynomials of maximum degree 2n+1. The most obvious choice, the monomials, is discouraged because of the inherent bad conditioning of the resulting Vandermonde system. For numerical representation the Newton basis is much better. In order to use this in a simple way, the nodes x_i are considered as ''double nodes''. In order to describe this and to use divided differences, we introduce new artificial variables ξ , η and set
  ξ_2i = ξ_2i+1 = x_i, i=0,...,n .
  Then the Newton basis is
  N_i(x) = Π _j=0,..,i-1 (x-ξ_j), i=0,...,2n+1
  with the empty product considered as the constant 1 .
  The function values are copied into the η variables:
  η_2i = η_2i+1 = y_i, i=0,...,n .
  The computation of the coefficients is done using the scheme of "divided differences" with a special treatment of the second column (index 1)
  c_i,0 = η_i, i=0,...,2n+1
  c_2i,1 = y_i⁽¹⁾, i=0,...,n
  c_2i+1,1 = (η_2i+2 - η_2i)/(ξ_2i+2 - ξ_2i), i=0,...,n-1
  c_i,j+1 = (c_i+1,j-c_i,j)/(x_j+i-x_i), j=1,..,2n, i=0,..,2n-j
  with the final form
  p(x) = ∑_{i=0 to 2n+1} c_0,i N_i(x)
  
• This scheme is obtained from the ordinary Newton interpolation scheme by coalescing two consecutive (different) nodes, sometimes also called the ''confluent'' case. The computational effort for the coefficients is O(n²) and O(n) for a single polynomial value, which can be done especially efficient using the common factors of the N_i.
• If y_i = f(x_i), y_i⁽¹⁾ = f'(x_i) and f is 2n+2 times continuously differentiable, then the interpolation error is
  f⁽²ⁿ⁺²⁾(x^~)/((2n+2)!)* ((x-x₀)*...*(x-x_n))²
  where x^~ is an unknown value between x and the x_i and hence
  <= M_2n+2/((2n+2)!)*(b-a)²ⁿ⁺²
  for an arbitrary distribution of the abscissae in [a,b], where M_2n+2 is an upper bound for the absolute value of the 2n+2-th derivative of f on [a,b]. If this derivative has no zeroes on [a,b] and m_2n+2>0 is a lower bound for its absolute value on [a,b], then there is at least one x in [a,b] where
  |f(x)-p(x)| >= m_2n+2/((2n+2)!)*2*((b-a)/4)²ⁿ⁺²
  whatever (optimal) distribution of abscissae one chooses. This shows that this method always gives good results if the interval span is small, but not in the case that there exists a (possibly complex) singularity of f whose distance from [a,b] is smaller than 1/(b-a) , because then the growth of the derivatives (in n) is stronger than this.

Input

• You either might interpolate a discrete data set of your own (getting then a plot ) or might choose a demonstration using predefined or self defined functions for the generation of the y and y⁽¹⁾ data.
• There are 5 predefined cases of functions with different growth properties of their derivatives. You also might define a function of your choice providing a piece of code following FORTRAN conventions.
• You also specify the degree 2n+1 by choosing n which is restricted to 1 <= n <= 50 !
• Then either you give a list of your (x,y,y⁽¹⁾) -data or
• specify a demonstration case:
• In this case you also have the choice how to distribute the abscissae in this interval: either equidistant (a common case in practice), hence
  a+i*(b-a)/n
  or in an optimal way (as far as it is independent of f), namely
  (a+b)/2+((b-a)/2)*cos(((2*i+1)/(2*n+2))*π),
  the so called Chebyshev abscissae, which are obtained by
  min_{x0,...,xn in [a,b]} max_{x in [a,b]} | Π _{j=0 to n}(x-x_j) | .
  Next you specify n, the function case, the interval . There is a special choice, namely setting y⁽¹⁾ =0 for all i with the Chebyshev abscissae and only (x,f) data. This is Fejers method and one can show that this defines an interpolation process by polynomials which is convergent in the max-norm for every continuous function on [a,b]. You will of course see that the convergence of this process is quite slow (why?)
• you also have the choice to see the scheme of divided differences and the evaluation formula of the polynomial.

Output

• In case of own discrete data a plot, showing these data and the interpolating polynomial
• If you have generated the data artificially, then you get a plot of f and one of the error curve f - p .
• In both cases you also have the choice to see the scheme of divided differences which gives you information about the growth of derivatives of f as long as these are no spoiled by roundoff. Remember that in each column of this scheme roundoff will be in order of ''function value size times machine precision didived by the product of the abscissa differences involved'', in case of equidistant data and f values in the order of one hence something like
  2^jε/h^j
  in column j, with ε the machine precision (2^-53 here).
• In both cases you also have the choice the see the formula which evaluates the polynomial explicitly.

Questions?!

• For which functions decreases the error with increasing degree and why?
• For which functions this appears just opposite (divergence) and why?
• What appears if you take the Chebyshev abscissae ?
• What takes place of you try a nonsmooth function?

To the input form