Chapter 14. Comprehensive Examples

In this chapter, we present some comprehensive and longer examples together with a brief introduction to the theoretical background and their complete implementation. By this, we want to show you how the concepts defined in this book are used in practice.

Polynomials

First, we will demonstrate the power of the Python constructs presented so far by designing a class for polynomials. We will give some theoretical background, which leads us to a list of requirements, and then we will give the code, with some comments.

Note, this class differs conceptually from the class numpy.poly1d.

Theoretical background

A polynomial: p(x) = a_n x ⁿ + a_n-1 x^n-1 +…+ a₁ x + a₀ is defined by its degree, its representation, and its coefficients. The polynomial representation shown in the preceding equation is called a monomial representation. In this representation, the polynomial is written as a linear combination of monomials, xⁱ . Alternatively, the polynomial can be written in:

• Newton representation with the coefficients c_i and n points, x₀, …, x_n-1:

  p(x) = c₀ + c₁ (x - x₀) + c₂ (x - x₀)(x-x₁) + ... + c_n(x - x₀) … (x - x_n-1)

• Lagrange representation with the coefficients y_i and n+1 points, x₀, … , x_n:

                     p(x) = y₀ l₀(x) + y₁ l₁(x) + … + y_n l_n(x)

  with...

The polynomial class

Let's now design a polynomial base class based on a monomial formulation of the polynomial. The polynomial can be initialized either by giving its coefficients with respect to the monomial basis or by giving a list of interpolation points, as follows:

import scipy.linalg as sl

class PolyNomial:
    base='monomial'
    def __init__(self,**args):
        if 'points' in args:
            self.points = array(args['points'])
            self.xi = self.points[:,0]
            self.coeff = self.point_2_coeff()
            self.degree = len(self.coeff)-1
        elif 'coeff' in args:
            self.coeff = array(args['coeff'])
            self.degree = len(self.coeff)-1
            self.points = self.coeff_2_point()
        else:
            self.points = array([[0,0]])
            self.xi = array([1.])
            self.coeff = self.point_2_coeff()
            self.degree = 0

The __init__...

Newton polynomial

The NewtonPolyNomial class defines a polynomial described with respect to the Newton basis. We let it inherit some common methods from the polynomial base class, for example, polynomial.plot, polynomial.zeros, and even parts of the __init__ method, by using the super command (refer to section Subclassing and Inheritance in Chapter 8, Classes):

class NewtonPolynomial(PolyNomial):
    base = 'Newton'
    def __init__(self,**args):
        if 'coeff' in args:
            try:
                self.xi = array(args['xi'])
            except KeyError: 
                raise ValueError('Coefficients need to be given'
                'together with abscissae values xi')
        super(NewtonPolynomial, self).__init__(**args)

Once the interpolation points are given, the computation of the coefficients is performed by:

def point_2_coeff(self):
    return array(list(self.divdiff()))

Here we used divided differences for computing the Newton...

Spectral clustering

An interesting application of eigenvectors is for clustering data. Using the eigenvectors of a matrix derived from a distance matrix, unlabelled data can be separated into groups. Spectral clustering methods get their name from the use of the spectrum of this matrix. A distance matrix for n elements (for example, the pairwise distance between data points) is an n × n symmetric matrix. Given such an n × n distance matrix M with distance values m_ij , we can create the Laplacian matrix of the data points as follows:

Here, I is the identity matrix and D is the diagonal matrix containing the row sums of M,

 

The data clusters are obtained from the eigenvectors of L. In the simplest case of data points with only two classes, the first eigenvector (that is, the one corresponding to the largest eigenvalue) is often enough to separate the data.

Here is an example for simple two-class clustering. The following code creates some 2D data points and clusters them based on the first...

Solving initial value problems

In this section, we will consider the mathematical task of numerically solving a system of ordinary equations for given initial values:

y'(t) = f(t, y)      y(t₀) = y₀∈ ℝⁿ

The solution of this problem is a function y. A numerical method aims at computing good approximations, y_i≈ y(t_i) at discrete points, the communications points t_i, within the interval of interest [t₀, t_e]. We collect the data that describes the problem in a class, as follows:

class IV_Problem:
    """
    Initial value problem (IVP) class
    """
    def __init__(self, rhs, y0, interval, name='IVP'):
        """
        rhs 'right hand side' function of the ordinary differential
                                                   equation f(t,y)
        y0 array with initial values
        interval start and end value of the interval of independent
        variables often initial and end time
        name descriptive name of the problem...

Summary

Most of what we explained in this book is bundled into the three longer examples of this chapter. They mimic code development and give prototypes, which you are encouraged to alter and confront with your own ideas.

You saw that code in scientific computing can have its own flavor due to its strong relation with mathematically defined algorithms and that it is often wise to keep the relationship between code and formula visible. Python has techniques for this, as you have seen.

Exercises

Ex. 1 → Implement a method __add__, which constructs a new polynomial p+q by adding two given polynomials p and q. In monomial form, polynomials are added by just adding the coefficients, whereas in Newton form, the coefficients depend on the abscissa x_i of the interpolation points. Before adding the coefficients of both polynomials, the polynomial q has to get new interpolation points with the property that their abscissa x_i coincides with those of p and the method __changepoints__ has to be provided for that. It should change the interpolation points and return a new set of coefficients.

Ex. 2 → Write conversion methods to convert a polynomial from Newton form into monomial form and vice versa.

Ex. 3 → Write a method called add_point that takes a polynomial q and a tuple (x,y) as parameters and returns a new polynomial that interpolates self.points and (x,y).

Ex. 4 → Write a class called LagrangePolynomial that implements polynomials in Lagrange form and inherits as much as...