Quadratic interpolation

The table below shows a list of data points.

n	xn	yn
0	1.04	2.71
1	1.51	2.55
2	1.83	2.82

A common task in applied math is interpolation: given a set of data points, construct a curve passing through those data points and use the curve to predict values of the curve at other points. In the current example we want to interpolate a curve through the three data points shown above. Since we have three data points, the simplest curve that passes through all three of those points will be a quadratic polynomial. Our job is to compute that polynomial and use it to predict the value of y at other points.

The first method to compute the polynomial starts with the assumption that the polynomial takes the form

p(x) = a x² + b x + c

Substituting the three data points into this curve gives a set of three equations in three unknowns.

These equations can be rewritten as a matrix-vector system.

If you know how to compute the inverse of a matrix, you can then solve this system for the coefficients of the polynomial.

Newton polynomials

Here is a second method to solve this problem. The method is based a slightly different assumption about the form of the polynomial we are trying to solve for, the so-called Newton form of the polynomial:

p(x) = a (x-x₀) (x - x₁) + b (x-x₀) + c

The algebra needed to solve for the coefficients of this polynomial is a lot simpler:

y₀ = p(x₀) = c

y₁ = p(x₁) = b (x₁ - x₀) + y₀

In most applications of this method we manipulate this expression a bit more to put it into a more convenient form for computation:

Doing arithmetic with Python

Now that we have reduced the computation of the polynomial to a set of formulas, the last thing left to do is to do the necessary arithmetic to compute p(x) for some desired x.

Below you will find the code for a Python program that can do the necessary calculations for us. Make a new project named Quadratic in PyCharm and add a file quadratic.py to the project. Paste the code below into that file.

x0 = 1.04
y0 = 2.71
x1 = 1.51
y1 = 2.55
x2 = 1.83
y2 = 2.82
diff10 = (y1-y0)/(x1-x0)
diff21 = (y2-y1)/(x2-x1)
c = y0
b = diff10
a = (diff21-diff10)/(x2-x0)
x = 1.7
p = a*(x-x1)*(x-x0)+b*(x-x0)+c
print('At x = 1.7, y is approximately =',p)

When you run the program it will print an approximation for p(1.7).

The example program demonstrates some of the basics of working with variables and arithmetic in Python. To use variables in a Python program we simple write a statement that assigns a value to the variable.

To do arithmetic with the values stored in the variables we use the operators +, -, *, /, and ** for addition, subtraction, multiplication, division, and exponentiation. For grouping we use the parenthesis notation you are used from mathematics.

Adding comments

Once a program grows beyond just a few lines, you may want to insert comments into the program to help explain what the various parts of the program are doing. To type a comment in Python, type the hash symbol # followed by the text of the comment. Python will ignore text of the comment - comments are used simply to help document what various parts of the program are doing.

Here is the program above again with comments inserted at appropriate locations to help document what is going on in the program.

# This is a program to construct a quadratic polynomial
# from a set of data points. We use the Newton method to
# construct the polynomial

# The initial set of data points.
x0 = 1.04
y0 = 2.71
x1 = 1.51
y1 = 2.55
x2 = 1.83
y2 = 2.82

# Constructing the coefficients of the polynomial
diff10 = (y1-y0)/(x1-x0)
diff21 = (y2-y1)/(x2-x1)
c = y0
b = diff10
a = (diff21-diff10)/(x2-x0)

# This is the value of x at which we want to evaluate
x = 1.7

# Compute the value of the polynomial at x using the
# Newton form of the polynomial
p = a*(x-x1)*(x-x0)+b*(x-x0)+c

# Print the result
print('At x = 1.7, y is approximately =',p)

Programming assignment

As an alternative to fitting a quadratic polynomial to the original set of three data points you can try doing two separate linear interpolations. For the first calculation, construct the equation of the line passing through the data points (x₀ , y₀) and (x₁ , y₁) and use that linear function to estimate the value that y would have at x = 1.7. For the second computation do the same thing, but instead construct a line that passes through (x₁ , y₁) and (x₂ , y₂).

Write a Python program that does both of these linear interpolations, and have your program print the two estimates for y(1.7).

To submit your work for grading, please email me the Python file that contains the code for your program as an attachment to an email message.