HiveBrain v1.2.0
Get Started
← Back to all entries
patternpythonMinor

Plotting different parameterized polynoms

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
plottingdifferentparameterizedpolynoms

Problem

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 lebesgu

Solution

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

Code Snippets

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 tschebyscheff_points(count):
    return np.cos(((2 * np.arange(count) + 1) * np.pi) / (2 * (count + 1)))
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")

Context

StackExchange Code Review Q#148451, answer score: 7

Revisions (0)

No revisions yet.