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