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")

• Vogel612: range doesn't handle floats that well (at least in python 3.5) ... so the replacement you propose isn't going to work (for that instance) there's a workaround with a list comprehension to have the range run over integers and just divide the elements so as to get the same floats. Also the equation you have in the answer isn't quite correct. the \$x_i\$ is actually just \$i\$...
• Peilonrayz: @Vogel612 Oh yeah, you're working with floats... Silly me, I think the main reasoning behind my point still stands, but I'll correct the domain. And yeah I'll fix the equation, thought you were passing a list when I was writing it. :)
• Peilonrayz: @MathiasEttinger I'll change it later, just in-case you, or anyone else, thinks of any other changes
• Vogel612: I've tried to run the code you put up here, but I keep getting a ValueError saying that the first dimension of the domain and the result of calling poly_calcmismatch. That's the case for both lebesgue_function and as_supporting_point_poly. I assume the numpy-reduce functions reduce all the dimensions into a scalar :/
• Peilonrayz: @Vogel612 I reverted to my original code, as np.multiply.reduce returns a generator if you pass one in... Which is odd. Sorry I didn't know about that. The above code now works, I've double checked.

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[\$.