Interpolation by rational functions

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 and reproducing these data, with the aim of evaluation or otherwise mathematical manipulation (say integration, differentiation etc). Here we have a somewhat tricky case: the "simple" function is chosen as a fraction of two polynomials of given maximum degree and the dimension of the independent variable is one.
• We have data x_i, y_i, i=0,...,n, n=l+m, (with the meaning of y_i = f(x_i), f not necessarily known), and want to compute a rational function of numerator degree at most l and denominator degree at most m (with the cases l=0 or m=0 provided) such that
  R(x_i) = y_i, i=0,...,n .
  The solvability of this problem is not universally guaranteed. We discuss this below.
• This problem can be solved in a variety of ways, depending on the choice of the representation of the function R(x)=p_l(x)/q_m(x).
• For strictly monotonic data the representation by a sum of a polynomial of degree l-m and a Thiele continued fraction multiplied by a polynomial of degree l-m+1 is a choice, requiring however l >= m, but for arbitrary input data this has often severe trouble since there can occur undefined subexpressions. Hence here we reduce the problem to the solution of a homogeneous linear system of equations which is always solvable.
• Let
  p_l(x) = ∑_{i=0 to l} a_ixⁱ
  and
  q_m(x) = ∑_{i=0 to m}b_ixⁱ
  Then
  R(x_i) = y_i i=0,...,l+m
  is equivalent to the linear system with l+m+1 rows and l+m+2 columns (unknowns)
  

  [	1	x0	...	x0l	-y0	...	-y0 x0m	]	[a0]	=	0
  [	1	x1	...	x1l	-y1	...	-y1 x1m	]	[a1]	=	0
  [	..	..	...	...	...	...	...	]	[...]	=	0
  [	..	..	...	...	...	...	...	]	[...]	=	0
  [	1	xn-1	...	xn-1l	-yn-1	...	-yn-1 xn-1m	]	[bm-1]	=	0
  [	1	xn	...	xnl	-yn	...	-yn xnm	]	[bm]	=	0

  The interpolation problem is solvable if and only if the system matrix has rank l+m+1 and the two polynomials p_l and q_m defined by a normalized solution have no common factor. This is a system with a generalized Vandermonde type of matrix and one may have good reason to suspect a severe illconditioning of this solution approach. We try to overcome this here by linearly transforming the x-range to [-1,1] and using the chebyshev polynomials of the first kind
  T₀(x) = 1 , T₁(x) = x , T_i(x) = 2 x T_i-1(x) - T_i-2(x) , i=2,...
  as basis of the polynomial space.
• We solve the homogeneous system with an automatic normalization of the solution by computing the singular values decomposition of the matrix and taking the column of the right singular vector matrix which corresponds to the smallest singular value (=0 of course) as the coefficient vector.
• We report the singular values such that you might judge the rank condition yourself. Indeed test cases can be constructed with several singular values zero.
• We also report the zeroes of the polynomials p_l and q_m, computed by the method of Hirano (which at least in theory can never fail). It is easy to construct cases there these polynomials have nearly a common factor.
• Even worse, the interpolation problem can be uniquely solvable but yielding a solution which exhibits poles inside the data interval. This problem cannot be overcome. In order to avoid this a much more involved problem must be solved: an approximation problem with an infinite number of linear constraints.
• In many cases the rational interpolation results in a much lower approximation error for the function than a polynomial approximation, e.g. if f itself exhibits poles or bounded asymptotes. Our testcases provided here demonstrate this clearly.
• If y_i = f(x_i) and f is m+l+1 times continuously differentiable then the interpolation error is
  ∏ _{i = 0 to l+m} (x-x_i)/((l+m+1)!) × (1/q_m(x)) × ∑_{j = 0 to m} (l+m+1 over j+l+1) f^(j+l+1)(θ)q_m^(m-j)(θ)
  where θ is an unknown value between x and the x_i. This shows that this method always gives good results if the interval span is small and |q_m| is sufficiently bounded away from zero, but otherwise it can not be excluded that the error grows even unbounded in the interval of consideration.

Input

• You either might interpolate a discrete data set of your own (getting then a plot showing the rational interpolant and your data) or might choose a demonstration using predefined or self defined functions for the generation of the y-data.
• You get a printout showing the singular values of the system matrix , of the coefficients of p_l and q_m in the standard monomial basis and for the original interval and the computed zeroes of numerator and denominator.
• There are 5 predefined cases of functions with different growth properties of their derivatives. You also might define a function of your choice using a piece of code following FORTRAN conventions.
• You also specify the degrees l and m which are restricted to 0 <= l,m <= 40 and l+m <= 40!
• Then either you give a list of your (x,y)-data or
• specify the interval from where the abscissae should be drawn.
• 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 solving
  min _{x0,...,xl+m in [a,b]} max_{x in [a,b]} | Π _{j=0 to l+m}(x-x_j) |.
  

Output

• The result list mentioned above.
• In case of own discrete data a plot, showing these data and the interpolating function.
• If you have generated the data artificially, then you get a plot of f and R and one of the error curve f(.) - R(.) . If the error is quite small, then you will see only one graph since one curve is plotted over the other. It may be that there is a pole inside the interval without the plot showing this: if there are almost identical zeroes of the numerator and the denominator our resolution (2000 points in the interval) might not suffice to resolve this. Hence always have a look on the table of zeroes! If a value of the interpolant is larger than 10 times the maximum interpolated value (in absolute values) then we have a gap in the plot!

Questions?!

• For which functions works the rational interpolation much better than the polynomial interpolation with the same degrees of freedom? (You can choose m=0!)
• For which functions this appears just opposite and why?
• What appears if you take the Chebyshev abscissae ? Are these much better here?
• What takes place of you increase the degrees with data which could be represented exactly with lower degrees?

To the input form