In this chapter, we will present some comprehensive and longer examples together with a brief introduction to the theoretical background and the examples' complete implementation. Here, we want to show you how the concepts defined in this book are used in practice.
The following topics are covered in this chapter:
First, we will demonstrate the power of the Python constructs presented so far by designing a class for polynomials.
Note that this class differs conceptually from the class numpy.poly1d.
We will give some theoretical background, which will lead us to a list of requirements, and then we will give the code, with some comments.
A polynomial 
is defined by its degree, representation, and 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 
.
Alternatively, the polynomial can be written in:
and 
points, 
:
and 
points, 
:
with the cardinal functions:
There are infinitely many representations, but we restrict ourselves here to these three typical ones.
A polynomial can be determined from interpolation conditions:
with the given distinct values 
and arbitrary values 
as input. In the Lagrange formulation, the interpolation polynomial is directly available, as its coefficients are the interpolation data. The coefficients for the interpolation polynomial in Newton representation can be obtained...
We can now formulate some programming tasks:
.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
import matplotlib.pyplot as mp 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...
Let's give some usage examples.
First, we create a polynomial instance from the given interpolation points:
p = PolyNomial(points=[(1,0),(2,3),(3,8)])
The polynomial's coefficients with respect to the monomial basis are available as an attribute of p:
p.coeff # returns array([ 1., 0., -1.]) (rounded)
This corresponds to the polynomial 
. The default plot of the polynomial, obtained by p.plot((-3.5,3.5)), results in the following figure (Figure 19.1):
Finally, we compute the zeros of the polynomial, which in this case are two real numbers:
pz = p.zeros() # returns array([-1.+0.j, 1.+0.j])
The result can be verified by evaluating the polynomial at these points:
p(pz) # returns array([0.+0.j, 0.+0.j])
The class NewtonPolyNomial 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 command super (see Section 8.5: Subclasses and inheritance):
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()))
We used divided differences for computing...
An interesting application of eigenvectors is for clustering data. Using the eigenvectors of a matrix derived from a distance matrix, unlabeled 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 
elements (for example, the pairwise distance between data points) is an 
symmetric matrix. Given such an 
distance matrix 
with distance values 
, we can create the Laplacian matrix of the data points as follows:
Here, 
is the identity matrix and 
is the diagonal matrix containing the row sums of 
:
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 of simple two-class clustering. The following code creates some 2D data points and clusters them...
In this section, we will consider the mathematical task of numerically solving a system of ordinary equations for given initial values:
.
The solution to this problem is a function 
. A numerical method computes approximations, 
at discrete communications points, 
, within the interval of interest 
. 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
"""
self.rhs...Most of what we explained in this book was bundled into the three longer examples in this chapter. These examples 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 relationship 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.
Ex. 1: Implement a method __add__ that constructs a new polynomial 
by adding two given polynomials 
and 
. In monomial form, polynomials are added by just adding the coefficients, whereas in Newton form, the coefficients depend on the abscissas 
of the interpolation points. Before adding the coefficients of both polynomials, the polynomial 
has to get new interpolation points with the property that their abscissas 
coincide with those of 
, 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 
as parameters and returns a new polynomial that interpolates self.points and 
.
Ex. 4: Write a class called LagrangePolynomial...
© 2021 Packt Publishing Limited All Rights Reserved
