Lab 07: Polynomials

Lab 07: Polynomials#

Examples (by instructor)#

In this lab, you’ll learn to work with polynomials in NumPy: creating and evaluating polynomial objects, fitting data with least-squares polynomial regression, and assessing fit quality with residuals, MAE, MSE, and \(R^2\).

Example 1: Polynomial Basics — np.poly1d, np.polyval#

A degree-\(n\) polynomial is:

\[p(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0\]

NumPy stores coefficients in descending order (highest power first).

Function

Purpose

np.poly1d(coeffs)

Create a callable polynomial object

np.polyval(coeffs, x)

Evaluate polynomial at \(x\)

np.polyfit(x, y, deg)

Least-squares fit — returns coefficients

Define \(p(x) = 2x^3 - 3x^2 + x - 5\) and:

  1. Print its poly1d representation

  2. Evaluate it at \(x = 2\) using both poly1d and polyval

  3. Print its derivative and its roots

Example 2: Root Finding — Chemical Equilibrium#

In chemical engineering, we often need to find where a polynomial equals zero. A classic example is solving an equilibrium expression that reduces to a polynomial equation.

Problem: The equilibrium conversion \(x\) for the gas-phase reaction \(A \rightleftharpoons 2B\) at constant \(T\) and \(P\) satisfies:

\[K_p = \frac{4x^2}{1 - x^2}\]

Rearranging for \(K_p = 2.0\):

\[K_p(1 - x^2) = 4x^2 \implies K_p = (4 + K_p)\,x^2 \implies (4 + K_p)\,x^2 - K_p = 0\]

This is a degree-2 polynomial in \(x\). We will:

  1. Express it as np.poly1d and find its roots via p.roots

  2. Keep only the physically meaningful root (\(0 < x < 1\))

Example 3: Fitting \(C_p(T)\) for CO\(_2\) with np.polyfit#

Below is NIST tabulated heat capacity data for CO\(_2\). Use np.polyfit to fit degree 2 polynomial, then:

  1. Plot fit against the data

  2. Compute and print MAE, MSE, and \(R^2\) for each degree

\[ \text{MAE} = \frac{1}{N}\sum|e_i|, \quad \text{MSE} = \frac{1}{N}\sum e_i^2, \quad R^2 = 1 - \frac{\sum e_i^2}{\sum(y_i - \bar{y})^2} \]

\(T\) (K)

300

400

500

600

700

800

900

1000

1100

1200

\(C_p\) (J/mol/K)

37.13

41.33

44.60

47.33

49.65

51.61

53.26

54.31

55.37

56.21

 
import numpy as np
import matplotlib.pyplot as plt

# NIST tabulated Cp data for CO2
T_data  = np.array([300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200], dtype=float)
Cp_data = np.array([37.13, 41.33, 44.60, 47.33, 49.65,
                    51.61, 53.26, 54.31, 55.37, 56.21])

Warm-Up: Syntax Practice#

Short exercises to get comfortable with the key NumPy polynomial functions before the main problems.

Exercise 1 — Define a polynomial.

Define \(q(x) = x^2 - 5x + 6\) using np.poly1d and print it. Then confirm that the coefficients stored inside match what you passed in.

import numpy as np

# Define q(x) = x^2 - 5x + 6
q = np.poly1d([___,  ___, ___])   # fill in the three coefficients

print(q)
print("Coefficients:", q.coeffs)  # should print [1, -5, 6]

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[2], line 6
      3 # Define q(x) = x^2 - 5x + 6
      4 q = np.poly1d([___,  ___, ___])   # fill in the three coefficients
----> 6 print(q)
      7 print("Coefficients:", q.coeffs)  # should print [1, -5, 6]

File ~/miniconda3/envs/chemcomp/lib/python3.10/site-packages/numpy/lib/polynomial.py:1284, in poly1d.__str__(self)
   1282 for k, coeff in enumerate(coeffs):
   1283     if not iscomplex(coeff):
-> 1284         coefstr = fmt_float(real(coeff))
   1285     elif real(coeff) == 0:
   1286         coefstr = '%sj' % fmt_float(imag(coeff))

File ~/miniconda3/envs/chemcomp/lib/python3.10/site-packages/numpy/lib/polynomial.py:1277, in poly1d.__str__.<locals>.fmt_float(q)
   1276 def fmt_float(q):
-> 1277     s = '%.4g' % q
   1278     if s.endswith('.0000'):
   1279         s = s[:-5]

TypeError: must be real number, not numpy.str_

Exercise 2 — Evaluate at a point.

Using the same \(q(x) = x^2 - 5x + 6\), evaluate \(q(2)\) and \(q(3)\) two ways: via the poly1d object and via np.polyval. Confirm both give the same answer, and check manually that the values make sense (hint: \(q(2) = 0\) and \(q(3) = 0\)).

import numpy as np

coeffs = [1, -5, 6]
q = np.poly1d(coeffs)

for x_val in [2, 3]:
    via_poly1d = ___          # call q like a function
    via_polyval = ___         # use np.polyval
    print(f"q({x_val}): poly1d = {via_poly1d},  polyval = {via_polyval}")

Exercise 3 — Find roots.

Find the roots of \(q(x) = x^2 - 5x + 6\) using p.roots. You should get two real roots. Print them and verify they match what you’d get by factoring: \(q(x) = (x-2)(x-3)\).

import numpy as np

q = np.poly1d([1, -5, 6])

roots = ___          # get the roots of q
print("Roots:", roots)

# Verify by evaluating q at each root — should both be ~0
for r in roots:
    print(f"q({r:.1f}) = {q(r):.6f}")

Exercise 4 — Fit a polynomial to data.

Five measurements of reaction rate \(r\) at different temperatures \(T\) are given below. Fit a degree-2 polynomial using np.polyfit, then evaluate it at \(T = 350\) K.

\(T\) (K)

300

325

350

375

400

\(r\) (mol/L/s)

1.2

2.3

3.8

5.7

8.1

 
import numpy as np

T_data = np.array([300, 325, 350, 375, 400], dtype=float)
r_data = np.array([1.2, 2.3, 3.8, 5.7, 8.1])

# Fit a degree-2 polynomial
coeffs = np.polyfit(___, ___, ___)   # x, y, degree
p = np.poly1d(coeffs)

print("Fitted polynomial:")
print(p)

# Evaluate at T = 350 K
T_query = 350.0
print(f"\nr({T_query} K) = {p(___):.4f} mol/L/s")   # call p at T_query

Practice Problems (by students)#

Problem 1: Polynomial Evaluation#

Consider the polynomial \(p(x) = 4x^3 - 2x^2 + 7x - 1\).

  1. Create it as a np.poly1d object and print its representation.

  2. Evaluate \(p(3)\) using both poly1d and np.polyval. Verify manually.

  3. Find all roots of \(p(x)\). Which roots are real? Which are complex?

  4. Plot \(p(x)\) over \(x \in [-2, 4]\) and mark the real root(s) with a red dot.

Hint: p.roots returns all roots (may be complex); filter with np.isreal() to find real ones.

Problem 2: Fitting Viscosity Data — Degrees of Freedom & Overfitting#

The dynamic viscosity of liquid water varies with temperature (11 data points):

\(T\) (°C)

0

10

20

30

40

50

60

70

80

90

100

\(\mu\) (mPa·s)

1.793

1.307

1.002

0.798

0.653

0.547

0.467

0.404

0.355

0.315

0.282

(a) Fill in the degrees of freedom table below. Recall: \(\text{DOF} = N - (n+1)\) where \(N\) = number of data points and \(n\) = polynomial degree.

Degree \(n\)

Parameters \(n+1\)

DOF = \(11 - (n+1)\)

Overdetermined?

1

3

5

10

(b) Fit polynomials of degree 1, 3, 5, and 10, plot all fits over \(T \in [-10, 120]\) °C, and print MAE and \(R^2\) for each. Then answer: which degrees overfit, and how can you tell?

Contents