Hermite interpolation in Python

Question

I have this program for calculating Hermite interpolation.

Problem is, that its behave really bad. interpolation for 35 Chebyshev nodes

This is chart for 35 Chebyshev nodes. If I put more points, peak on the beginning will be higher(its about 10^7 with this amount of nodes).

I interpolated this same function using Lagrange method (green, its shifted so it can be seen) and as you can see it looks fine.

Here is the code:

def hermit_interpolate(input): #input is list of tuples [(x1,y1),(x2,y2)...] xi are Chebyshev nodes
points = [(input[0][0], input[0][1] - 0), (input[0][0], calculate_f_p_x(input[0][0]))] #"input[0][1] - 0" this is just to change type of second element
                                                                                    #calculate_f_p_x returns value of derivative
for k in range(1, len(input)): #Divided differences and derivatives in one list alternately
    points.append((input[k][0], (input[k][1] - input[k - 1][1]) / (
        input[k][0] - input[k - 1][0])))
    points.append((input[k][0], calculate_f_p_x(input[k][0])))
x, c = zip(*points)
x = list(x)
c = list(c)
n = len(points)
for i in range(2, n): #calculating factors
    for j in range(n - 1, i - 1, -1):
        c[j] = (c[j] - c[j - 1]) / (x[j] - x[j - i])

def result_polynomial(xpoint): #here is function to calculate value for given x
    val = c[0]
    factor = 1.0
    for l in range(1, n):
        factor *= (xpoint - x[l - 1])
        val += (c[l] * factor)
    return val

return result_polynomial

I can't seen what's wrong here. Thanks!

It had enormous errors also. Maybe it was something wrong with data preprocessing. Anyway, I wrote its again. This algorithm is less optimal as it use matrix instead of vector, but it works. I put proper code in answer if anyone need it — Purple, Apr 23 '15 at 16:17
Not sure, but there's a pretty good chance the errors may have been introduced by integer division, so if you have integers in your data and you're not taking special steps, the factor calculation where you divide by a difference may be truncated. Try using `from __future__ import division` right at the beginning of your code and see if the errors go away. — chthonicdaemon, Apr 28 '15 at 10:43

Purple · Accepted Answer · 2022-05-14T18:24:40.407

This code actually works:

def hermit_interpolate(input):  #input is list of tuples [(x1,y1),(x2,y2),...,(xn,yn)] xi are Chebyshev nodes
  n = len(input)
  points = numpy.zeros(shape=(2 * n + 1, 2 * n + 1))
  X, Y = zip(*input)
  X = list(X)
  Y = list(Y)

  for i in range(0, 2 * n, 2):
      points[i][0] = X[i / 2]
      points[i + 1][0] = X[i / 2]
      points[i][1] = Y[i / 2]
      points[i + 1][1] = Y[i / 2]

  for i in range(2, 2 * n + 1):
      for j in range(1 + (i - 2), 2 * n):
          if i == 2 and j % 2 == 1:
              points[j][i] = calculate_f_p_x(X[j / 2]);

          else:
              points[j][i] = (points[j][i - 1] - points[j - 1][i - 1]) / (
                  points[j][0] - points[(j - 1) - (i - 2)][0])

  def result_polynomial(xpoint):  #here is function to calculate value for given x
      val = 0
      for i in range(0, 2 * n):
          factor = 1.
          j = 0
          while j < i:
              factor *= (xpoint - X[j / 2])
              if j + 1 != i:
                  factor *= (xpoint - X[j / 2])
                  j += 1
              j += 1
          val += factor * points[i][i + 1]
      return val
return result_polynomia

Where is the return for the hermit_interpolate function? ans what is calculate_f_p_x ? — Oren, Dec 15 '15 at 16:15
This code have broken indentation. It's one function `hermit_interpolate` and it returns `result_polynomial` function. `Calculate_f_p_x` returns value of derivative. — Purple, Jan 13 '16 at 14:02

Hermite interpolation in Python

1 Answers1