Polynomials and fitting with NumPy

What you will be able to do

After this lesson you can

• store, evaluate, differentiate and solve a polynomial with NumPy,
• read complex roots and pick the real ones, and
• fit a straight line to measured data with np.polyfit and judge the fit.

A polynomial is a row of coefficients

You know that $x^3 - 6x^2 + 11x - 6$ is made of the coefficients $1, -6, 11, -6$. NumPy stores a polynomial in just this way: as an array of the coefficients, from the highest power down. The degree is one less than the length of the array. A missing power must be written with a zero: $x^3 - 2$ is [1, 0, 0, -2].

Predict first. What are the values of the polynomial at 1, 2 and 3? Then run it.

The poly1d object

Roots that are not real

To keep the real roots only, we test whether the imaginary part is close to zero: roots[np.isclose(roots.imag, 0)].real. For this polynomial the answer is an empty array, and the list is empty, so there are no real roots. Use np.isclose, and not == 0, because of the tiny errors in the last digits.

Fitting a line to measurements

Data from an experiment has scatter: no point is exactly on the line. A spring is stretched with five different forces, and the extension is measured. Hooke's law says that the force is proportional to the extension. How can we find the best straight line, and the spring constant?

Predict first. The force of 4.9 N gives an extension of 0.10 m. What is the rough spring constant? Then run it.

• np.polyfit(extension, force, 1) finds the straight line of degree 1 that fits the points best, and gives two numbers: the slope first, then the intercept. (The highest power comes first, as always.) The method is called least squares: it chooses the line for which the sum of the squares of the vertical distances from the points is the smallest.
• The slope is 48.0 N/m: this is the spring constant. The intercept is 0.16 N: for a perfect spring it should be 0, and a small number like this one means that the measurements have a small error.
• The residuals are the differences between the data and the line: $-0.02, 0.02, -0.04, 0.10, -0.06$. They are small, and they are not all on one side. Look at residuals for a pattern: a curved pattern would say that a straight line is the wrong model.
• With the line, we predict the force at 0.07 m: 3.52 N.
• $R^2$ (R squared) is a number up to 1 that tells how much of the scatter of the data the line explains. Here it is 0.9983, so the line explains 99.8 percent. A value near 1 says the fit is good, but it does not prove that the model is right. Always look at the graph too.

With a higher degree, np.polyfit(x, y, 2) fits a parabola, and so on. Do not use a high degree to make the line pass through every point. With 5 points, a polynomial of degree 4 passes through all of them exactly, but it wiggles between them and means nothing. A good model is as simple as the physics allows.

"Best" here has an exact meaning: the smallest sum of squared misses. That is what least squares means.

Common mistakes

• Giving the coefficients from the lowest power. NumPy wants the highest power first.
• Leaving out a zero coefficient. $x^3 - 2$ is [1, 0, 0, -2].
• Comparing a complex root with zero by ==. Use np.isclose.
• Reading the output of polyfit in the wrong order. The slope first, then the intercept.
• Using a high degree to chase every point. The model learns the noise.
• Judging a fit by $R^2$ alone. Look at the graph and the residuals.

Practice

Start with two quick questions. Each New numbers gives you a fresh one, so you can practise as long as you like.

Then the labs and the problems:

1. Warm-up. Open the lab Polynomial tools. You write functions for the real roots and the turning points.
2. Build. Open the lab Fit a line, the problem Roots and factors and the problem Fit a spring.

Remember these

• A polynomial is an array of coefficients, from the highest power down.
• np.polyval, np.roots, np.poly, np.polyder and np.poly1d do the work.
• Roots can be complex. Keep the real ones with np.isclose(roots.imag, 0).
• np.polyfit(x, y, 1) gives the slope and then the intercept of the best line.
• Look at the residuals and the graph, and keep the model simple.

Check yourself

Try each one on your own first. If you are stuck, open a hint. Then check with Show the answer.

1. What are the coefficients of $2x^3 - x + 4$ in NumPy?

2. What does np.roots([1, 0, 1]) give, and why?

3. For np.polyfit(x, y, 1) giving [3.0, 2.0], what is the line?

4. Why not fit a polynomial of degree 9 to 10 measurements?