Plotting different parameterized polynoms

Question

For a university assignment I had to plot different polynomial functions depending on one single parameter. That parameter gave the number of supporting points to interpolate a given function in the domain \$\lbrack -1.0, 1.0 \rbrack \$.

The supporting points were to be calculated in two different fashions. Either they were to be equidistant or be Chebyshev nodes. The given definitions were:

$$ x_i = \frac{2i}{n} - 1 , \quad x_i = \cos \frac{(2i + 1)\pi}{2(n + 1)} $$

The plots are to be handed in in a pdf. The polynomial functions I had to calculate were given as:

\$\Phi_n(x) = \underset{i \neq j}{\underset{i = 0}{\overset{n}{\Pi}}} (x - x_i) \$ and the slightly more complicated \$\lambda(x) = \underset{i = 0}{\overset{n}{\Sigma}} \lvert l_{i,n}(x) \rvert \$. Here \$l_{i,n}(x)\$ denotes a Lagrange polynomial. I'll just stop torturing you with math definitions, (because I'm reasonably sure I'm able to copy a formula from a script into code).

Note that \$\Phi_n\$ is called "Supporting point polynomial" and \$\lambda\$ is called "Lebesgue function" in the assignment.

So without further ado, here's my code.
Note that maintainabiltiy for future use is not a concern, so if you want you can mention docstrings and variable names, but those points don't really help me :)

import numpy as np import matplotlib.pyplot as plt def equidistant_points(count): points = [] for i in np.arange(0, count): points.append((2 * i / count) - 1) return points def tschebyscheff_points(count): points = [] for i in np.arange(0, count): points.append(np.cos(((2 * i + 1) * np.pi) / (2 * (count + 1)))) return points def as_supporting_point_poly(points, x): poly = 1 for point in points: poly = poly * (x - point) return poly def lagrange_poly(points, j, x): poly = 1 for i in range(0, len(points)): if (i != j): poly = poly * ((x - points[i]) / (points[j] - points[i])) return poly def lebesgue_function(points, x): leb = 0 for i in range(0, len(points)): leb = leb + np.fabs(lagrange_poly(points, i, x)) return leb def plot_and_save(n, x, poly_calc, name): equi = plt.plot(x, poly_calc(equidistant_points(n), x), 'r-', label='Äquidistante Stützstellen') tsch = plt.plot(x, poly_calc(tschebyscheff_points(n), x), 'g-', label='Tschebyscheff-Stützstellen') plt.xlabel("x") plt.ylabel("y") plt.title(name + " mit n = " + str(n)) plt.grid(True) plt.legend(loc='upper right', shadow=False) plt.savefig("Aufg1"+ poly_calc.__name__ + str(n) + ".png") plt.show() if __name__== '__main__': domain = np.arange(-1.0, 1.0001, 0.0001) plot_and_save(8, domain, as_supporting_point_poly, "Stützstellenpolynome") plot_and_save(20, domain, as_supporting_point_poly, "Stützstellenpolynome") plot_and_save(8, domain, lebesgue_function, "Lebesgue-Funktion") plot_and_save(20, domain, lebesgue_function, "Lebesgue-Funktion") 

I'm especially interested in ways to make the calculation of the supporting points in equidistant_points and tschebyscheff_points cleaner.

Answer 1

You're not really using numpy at the moment. If we use the simple 'translation table' below, your code would work if you just replaced the NumPy function with the Python equivalent:

• np.arange ->range, assuming your domain is just integers,
• np.fabs ->abs,
• np.cos ->math.cos, and,
• np.pi ->math.pi.

Instead you want to take advantage of NumPy. Take tschebyscheff_points, you have the equation:

$$\cos(\frac{\pi(2i + 1)}{2(\text{count} + 1)})$$

But your Python code is:

def tschebyscheff_points(count): points = [] for i in np.arange(0, count): points.append(np.cos(((2 * i + 1) * np.pi) / (2 * (count + 1)))) return points 

Yes it contains the equation, but with numpy you can just write the equation:

def tschebyscheff_points(count): return np.cos(((2 * np.arange(count) + 1) * np.pi) / (2 * (count + 1))) 

This significantly improves both performance, and readability. As you only need to read the equation.

---

I'd also change your code to use comprehensions. lebesgue_function should use sum as writing the addition yourself is WET. And in as_supporting_point_poly and lagrange_poly you should factor out the multiplication into a product function.

---

I'm not too good with NumPy and matplotlib so I can't really help with improving the display of the data. But the above should get you to the following code. Note, that I have two lagrange_poly as I don't know if a pure Python function is better than the NumPy equivalent.

import numpy as np import matplotlib.pyplot as plt from operator import mul from functools import reduce def product(it, initial=1): return reduce(mul, it, initial) def equidistant_points(count): return np.arange(count) * 2 / count - 1 def tschebyscheff_points(count): return np.cos(((2 * np.arange(count) + 1) * np.pi) / (2 * (count + 1))) def as_supporting_point_poly(points, x): return product(x - point for point in points) def lagrange_poly(points, j, x): point_j = points[j] p = np.delete(points, j) return product((x - p) / (point_j - p)) def lagrange_poly(points, j, x): return product( (x - point) / (points[j] - point) for i, point in enumerate(points) if i != j ) def lebesgue_function(points, x): return sum( np.fabs(lagrange_poly(points, i, x)) for i in range(len(points)) ) def plot_and_save(n, x, poly_calc, name): equi = plt.plot(x, poly_calc(equidistant_points(n), x), 'r-', label='Äquidistante Stützstellen') tsch = plt.plot(x, poly_calc(tschebyscheff_points(n), x), 'g-', label='Tschebyscheff-Stützstellen') plt.xlabel("x") plt.ylabel("y") plt.title(name + " mit n = " + str(n)) plt.grid(True) plt.legend(loc='upper right', shadow=False) plt.savefig("Aufg1"+ poly_calc.__name__ + str(n) + ".png") plt.show() if __name__== '__main__': domain = np.arange(-1.0, 1.0001, 0.0001) plot_and_save(8, domain, as_supporting_point_poly, "Stützstellenpolynome") plot_and_save(20, domain, as_supporting_point_poly, "Stützstellenpolynome") plot_and_save(8, domain, lebesgue_function, "Lebesgue-Funktion") plot_and_save(20, domain, lebesgue_function, "Lebesgue-Funktion") 

Answer 2

Here's my reworking of the calculation functions. Graphs look similar (didn't check details) when using the arange version of equidistant, but I think the linspace version is more balanced.

The main thing when trying to use numpy is keeping track of dimensions. In this case points is an array of 8-20 values, x an array of 2001.

def equidistant_points(count): # count values from -1 to 1-eps # slightly different end point handling #return np.arange(count)/count*2 - 1 return np.linspace(-1,1,count) def tschebyscheff_points(count): i = np.arange(count) return (np.cos(((2 * i + 1) * np.pi) / (2 * (count + 1)))) def as_supporting_point_poly(points, x): return np.product(x - points[:,None], axis=0) def lagrange_poly(points, j, x): # print('lag', points.shape, j, x.shape) pts = points[np.arange(points.shape[0])!=j] # skip the i==j point return np.product((x - pts[:,None])/(points[j]-pts[:,None]), axis=0) def lebesgue_function(points, x): # sum over i of prod over j for j!=i leb = 0 for i in range(points.shape[0]): leb += np.fabs(lagrange_poly(points, i, x)) return leb 

The lebesgue_function probably can be changed to avoid the loop (though looping over 8-20 items is better than 2000). But I haven't worked out the dimensional details. Avoiding that i==j divide by 0 is the main sticking point.

===================

I have figured out a vectorized (non-iterative) version of the lebesque:

def lebesgue_function(points, x): # sum over i of prod over j for j!=i # loop version arr = np.array([lagrange_poly(points, i, x) for i in range(points.shape[0])]) leb = np.abs(arr).sum(axis=0) return leb def lebesgue_function1(points, x): # full array version xp = x - points[:,None] # 8x2001 pp = points - points[:,None] # 8x8, diag 0s with np.errstate(invalid='ignore', divide='ignore'): xpp = xp[:,None,:]/pp[:,:,None] # 8x8x2001 n = np.arange(points.shape[0]) xpp[n,n,:] = 1 # so nan, inf don't affect prod leb = np.abs(xpp.prod(axis=0)).sum(axis=0) return leb 

with x (2000,), and small points (e.g. 10), the looped version is faster. points has to be in the 40 range to be faster. I had to play with the errstate to ignore the divide by 0 errors (which put nan and inf in xpp). Getting the prod axis right also took a bit of experimentation.

Answer 3

Besides @Peilonrayz excellent advice to "use numpy", I'd like to point out something that might be an issue with your supporting points.

From the equations, it seems they are defined for \$i \in [0;n]\$ but in your code, you only support \$i \in [0;n[\$ since it is the default behaviour of np.arange.

So I'd change your 2 functions defining supporting points to use

np.arange(count+1) 

instead of np.arange(0, count).

This will result in supporting points defined on \$[-1;1]\$ rather than \$[-1;1[\$.