""" A collection of functions dealing with Chebyshev polynomials and series. The argument naming conventions observed herein are as follows: x : array of points a,b : real interval where sample points are located. This interval is affinely mapped to the interval [-1,1]. nord : degree + 1, the nord of a Chebyshev polynomial or the number of standard or modified points in an interval. coef : 1D array of coefficients, [c_0,...,c_n], or an array where the coefficients run down the first index. The other convention is that of standard and modified Chebyshev points. On the interval [-1,1] these are defined as follows. standard points: x := cos(pi*(i + .5)/nord), 0 <= i < nord. The Chebyshev polynomial values at the standard points form orthogonal vectors. npote that the end points are not included, so the standard points are useful if something funny, i.e., a singularity, happens at one of the end points. modified points: x := cos(pi*i/(nord - 1)), 0 <= i < nord. The Chebyshev polynomial values at the modified points form orthogonal vectors if the end points are weighted by .5 in forming the inner product. npote that the end points of the interval are always included, so the nord must be > 1. ---------------------------------------------------------------------------- Included Functions : chebyshev_values(x, nord, a=-1, b=1) Evaluate Chebyshev polynomials at the points x. series_values(x, coef, a=-1, b=1) Evaluate a Chebyshev series at x. series_derivative(coef, a=-1, b=1) Derivative series of a Chebyshev series. series_integral(coef, a=-1, b=1) Integral series of a Chebyshev series. series_to_power_series(x, coef, a=-1, b=1) Convert a Chebyshev series to a power series about x. standard_points(nord, a=-1, b=1) Standard Chebyshev points in [a,b]. standard_series(nord) Matrix to convert function samples to a Chebyshev series. standard_values(x, nord, nder=0, a=-1, b=1) Matrix to interpolate function samples. standard_derivative(nord, a=-1, b=1) Matrix to convert function samples to derivative samples. standard_integral(nord, a=-1, b=1) Matrix to convert function samples to integral samples. modified_points(nord, a=-1, b=1) Modified Chebyshev points in [a,b] modified_series(nord) Matrix to convert function samples at modified points to a Chebyshev series. modified_values(x, nord, nder=0, a=-1, b=1) Matrix to interpolate function samples. modified_derivative(nord, a=-1, b=1) Matrix to convert function samples to derivative samples. modified_integral(nord, a=-1, b=1) Copyright December 16, 2006 by Charles R. Harris Email: charlesr.harris@gmail.com License: BSD """ import numpy as np # some private helper functions def _map_from_interval(x,a,b) : """F(x), where F: [a,b] -> [-1,1].""" return (x - (b + a)/2.0)*(2.0/(b - a)) def _map_to_interval(x,a,b) : """F(x), where F: [-1,1] -> [a,b].""" return x*((b - a)/2.0) + (b + a)/2.0 # main functions def chebyshev_values(x, nord, a=-1.0, b=1.0) : """Evaluate Chebyshev functions at the points x. Evaluates all Chebyshev functions of degree < nord at the points x. The function values are determined by the recursion T_0 = 1 T_1 = x T_n = 2x*T_{n-1} + T_{n-2} Before evaluation the points are shifted and scaled so that the interval [a,b] maps to [-1,1]. Arguments: x -- scalar or 1D array of points where the values will be computed. nord -- evaluate for all Chebyshev functions of degree < nord. a,b -- interval that maps to [-1,1]. Returns: Array of values. The last index corresponds to the degree of the chebyshev polynomial, the leading indices to the indices of x. The result has shape x.shape() + (nord,). Exceptions: ValueError -- if x is empty. ValueError -- if nord < 1 """ x = np.asarray(x, dtype=np.double) if nord < 1 : raise ValueError, 'nord must be >= 1' if x.size == 0 : raise ValueError, 'x array is empty' x = _map_from_interval(x, a, b) x2 = x*2.0 val = np.empty(x.shape + (nord,)) if nord == 1 : val[...,0] = 1. else : val[...,0] = 1. val[...,1] = x for i in range(2,nord) : val[...,i] = x2*val[...,i-1] - val[...,i-2] return val def series_values(x, coef, a=-1.0, b=1.0) : """Evaluate a Chebyshev series at x. The Chebyshev series is evaluated using inverse recursion - a generalized version of the Horner method for power series. The x values are scaled and shifted from the interval [a,b] to the interval [-1, 1] and the chebyshev series defined by coef is evaluated at the resulting points. Arguments: x -- scalar or array of points at which to evaluate the series. coef -- array, coef[i] contains the i'th coefficients. a,b -- interval that maps to [-1,1] Returns: the Chebyshev series evaluated at the point(s) x. Exceptions: ValueError -- if x is empty. ValueError -- if coef is empty. """ x = np.asarray(x, dtype=np.double) c = np.array(coef, dtype=np.double, copy=0, ndmin=1) if x.size == 0 : raise ValueError, 'x array is empty' if c.size == 0 : raise ValueError, 'coef array is empty' x = _map_from_interval(x, a, b) n = c.shape[0] s = x.shape + c.shape[1:] x = x.reshape(x.size,-1) c = c.reshape(n,-1) if n == 1 : c1 = np.zeros_like(c[0]) c0 = c[0] else : c1 = c[-1] c0 = c[-2] x2 = x*2.0 for i in range(3, n+1) : t = x2*c1 + c0 c0 = c[-i] - c1 c1 = t return (x*c1 + c0).reshape(s) def series_derivative(coef, a=-1.0, b=1.0) : """Derivative series of a Chebyshev series. Computes the Chebyshev series expansion of the derivative of a function defined by a Chebyshev series defined in the interval [a,b]. This is based on the recursion D'_n = T_{n-1} + D'_{n-2}, n > 2 D'_2 = T_1 D'_1 = T_0/2 where D_n = T_n/2n, n > 0. The expansion is assumed to hold for a function defined in the interval [a,b]. Arguments: coef -- Array, coef[i] contains the i'th coefficients. a,b -- interval that maps to [-1,1] Returns: Array of coefficients. Exceptions: ValueError -- if coef is empty or has rank > 2. """ c = np.array(coef, dtype=np.double, ndmin=1) if c.size == 0 : raise ValueError, "coef array is empty" n = c.shape[0] if n == 1 : dc = np.zeros(c.shape) else : dc = np.zeros((n-1,) + c.shape[1:]) c *= 4.0/(b - a) for n in range(n-1, 2, -1) : dc[n-1] += c[n]*n dc[n-3] += dc[n-1] for n in range(n-1, 0, -1) : dc[n-1] += c[n]*n dc[0] *= 0.5 return dc def series_integral(coef, a=-1.0, b=1.0) : """Integral series of a Chebyshev series. Computes the Chebyshev series expansion of the indefinite integral of a function defined by a Chebyshev series in the interval [a,b]. The lower bound of the integral is a, so that the series evaluates to zero at that point. The algorithm relies on the recursion \int T_{n} = D_{n+1} - D_{n-1}, n > 1 \int T_1 = D_2 \int T_0 = 2*D_1, where D_n = T_n/2n, n > 0. Arguments: coef -- Array, coef[i] contains the coefficients of order i. a,b -- interval that maps to [-1,1] Returns: Array of coefficients. Exceptions: ValueError -- if coef is empty. """ c = np.array(coef, dtype=np.double, ndmin=1) if c.size == 0 : raise ValueError, 'coef array is empty' n = len(c) ic = np.zeros((n+1,) + c.shape[1:]) c *= (b - a)*0.5 for n in range(n-1,0,-1) : ic[n-1] -= c[n] ic[n+1] += c[n] ic[n+1] *= 0.5/(n+1) ic[1] = ic[1]*0.5 + c[0] # make series evaluate to zero at -1 ic[0] = ic[1::2].sum(0) - ic[2::2].sum(0) return ic def series_to_power_series(x, coef, a=-1.0, b=1.0) : """Convert Chebyshev series to power series about x. Computes the power series expansion about the point x of a function defined by a Chebyshev series in the interval [a,b]. The computation proceeds by repeated series differentiation and evaluation at x to determine the coefficients in the Taylors series. npumerically, using the power series is a bad idea, as much cancellation takes place when the series is evaluated. nponetheless, for series of low nord, the results can be useful. Arguments: x -- scalar point at which the expansion takes place a,b -- interval that maps to [-1,1] coef -- vector of Chebyshev coefficients. Returns: coef.shape[0] + x.shape + coef.shape[1:] Exceptions: ValueError -- if x is not scalar ValueError -- if coef is empty or has rank > 2 """ x = float(x) c = np.array(coef, dtype=np.double, ndmin=1) if c.size == 0 : raise ValueError, 'coef array is empty' n = len(c) ps = np.zeros(c.shape) ps[0] = series_values(x, c, a, b) for i in range(1,n) : c = series_derivative(c, a, b)*(1.0/i) ps[i] = series_values(x, c, a, b) return ps def standard_points(nord, a=-1.0, b=1.0) : """Standard Chebyshev points in [a,b]. These are the points cos(pi*(i + .5)/nord), 0 <= i < nord in the interval [-1,1] shifted and scaled to the interval [a,b]. The Chebyshev functions form orthogonal vectors when evaluated at these points. Arguments: a,b -- interval that maps to [-1,1] nord -- integer, number of standard sample points (nord >= 1) Returns: 1D array of standard sample points norded from a to b Exceptions: ValueError -- if nord < 1 """ if nord < 1 : raise ValueError, 'nord must be >= 1' phi = (np.arange(nord) + .5)*(np.pi/nord) # reverse a,b to get points in increasing order pts = _map_to_interval(np.cos(phi),b,a) return pts def standard_series(nord) : """Matrix to convert function samples to a Chebyshev series. The matrix inverse of a square matrix whose columns are Chebyshev functions evaluated at the standard Chebyshev points in the interval [-1,1]. A function sampled at the standard points in any interval can be expanded in a Chebyshev series by multiplying the resulting vector of values by this matrix. It is the matrix of the DCT-2 transform. Arguments: nord -- integer, number of standard sample points (nord >= 1) Returns: 2D square array of the inverse Exceptions: ValueError -- if nord < 1 """ if nord < 1 : raise ValueError, 'nord must be >= 1' x = standard_points(nord) inv = np.transpose(chebyshev_values(x, nord))*(2.0/nord) inv[0] *= .5 return inv def standard_values(x, nord, nder=0, a=-1.0, b=1.0) : """Matrix to interpolate function samples. Returns a matrix on the interval [a,b] of the given nord. If a function on [a,b] is sampled at nord standard points then multiplication of the sample values by this matrix returns the values of the function at the points x after the operater D^nder is applied, where D is the derivative operator. npegative values of nder imply integration |nder| times, a zero value is the identity, and positive values of nder imply |nder| derivatives. Arguments: x -- points at which the values are computed. nord -- integer, number of standard sample points (nord >= 1) nder -- integer, number of derivatives to take. a,b -- interval that maps to [-1,1] Returns: Matrix that converts function values at the standard points to the values of D^nder(f) at the points x. Exceptions: ValueError -- nord must be >= 1 """ nops = int(nder) coef = standard_series(nord) if nops > 0 : for i in range(nops) : coef = series_derivative(coef, a, b) elif nops < 0 : for i in range(-nops) : coef = series_integral(coef, a, b) return series_values(x, coef, a, b) def standard_derivative(nord, a=-1.0, b=1.0) : """Matrix to convert function samples to derivative samples. Returns the derivative matrix on the interval [a,b] of the given nord. If a function on [a,b] is sampled at nord standard points then multiplication of the sample values by this matrix returns the derivative of the function at the standard points. Can be used to convert DEs to matrix equations. Arguments: nord -- integer, number of standard sample points (nord >= 1) a,b -- interval that maps to [-1,1] Returns: Matrix to convert function values at the standard points to the derivative values at the standard points. Exceptions: ValueError -- nord must be >= 1 ValueError -- der must be an integer """ x = standard_points(nord, a, b) return standard_values(x, nord, 1, a, b) def standard_integral(nord, a=-1.0, b=1.0) : """Matrix to convert function samples to integral samples. Returns the integral matrix on the interval [a,b] of the given nord. If a function on [a,b] is sampled at nord standard points then multiplication of the sample values by this matrix returns the integral of the function from a to the standard points. Can be used to convert integral equations to matrix equations. Arguments: nord -- integer, number of standard sample points (nord >= 1) a,b -- interval that maps to [-1,1] Returns: Matrix to convert function values at the standard points to the integral values at the standard points. The integral starts at a. Exceptions: ValueError -- nord must be >= 1 """ x = standard_points(nord, a, b) return standard_values(x, nord, -1, a, b) def modified_points(nord, a=-1.0, b=1.0) : """Modified Chebyshev points in [a,b]. These are the points cos(pi*i/(nord - 1)), 0 <= i < nord, shifted and scaled to the interval [a,b]. The Chebyshev functions form orthogonal vectors when evaluated at these points if the inner product is weighted by .5 at the end points. Arguments: a,b -- interval that maps to [-1,1] nord -- integer, number of modified sample points (nord >= 2) Returns: vector of sample points spaced from a to b Exceptions: ValueError -- if nord < 2. """ if nord < 2 : raise ValueError, 'nord must be >= 2' phi = np.arange(nord)*(np.pi/(nord - 1)) # reverse a,b to get points in increasing order pts = _map_to_interval(np.cos(phi), b, a) return pts def modified_series(nord) : """Matrix to convert function samples to a Chebyshev series. Returns the matrix inverse of a square matrix whose columns are Chebyshev functions evaluated at the modified Chebyshev points in the interval [-1,1]. A function sampled at the modified points in any interval can be expanded in a Chebyshev series by multiplying the resulting vector of values by this matrix. It is the matrix of the DCT-1 transform. Arguments: nord -- integer, number of samples. The degree of the fit is nord - 1. Returns: 2D square array of the inverse Exceptions: ValueError -- if nord < 2 """ if nord < 2 : raise ValueError, 'nord must be >= 2' x = modified_points(nord) inv = np.transpose(chebyshev_values(x, nord))*(2.0/(nord - 1)) inv[ 0,:] *= .5 inv[-1,:] *= .5 inv[:, 0] *= .5 inv[:,-1] *= .5 return inv def modified_values(x, nord, nder=0, a=-1.0, b=1.0) : """Matrix to interpolate function samples. Returns a matrix on the interval [a,b] of the given nord. If a function on [a,b] is sampled at nord modified points then multiplication of the sample values by this matrix returns the values of the function at the points x after the operater D^nder is applied, where D is the derivative operator. npegative values of nder imply integration |nder| times, a zero value is the identity, and positive values of nder imply |nder| derivatives. Arguments: x -- points at which the values are computed. nord -- integer, number of modified sample points (nord >= 1) nder -- integer, number of derivatives to take. a,b -- interval that maps to [-1,1] Returns: Matrix that converts function values at the modified points to the values of D^nder(f) at the points x. Exceptions: ValueError -- nord must be >= 2 """ nops = int(nder) coef = modified_series(nord) if nops > 0 : for i in range(nops) : coef = series_derivative(coef, a, b) elif nops < 0 : for i in range(-nops) : coef = series_integral(coef, a, b) return series_values(x, coef, a, b) def modified_derivative(nord, a=-1.0, b=1.0) : """Matrix to convert function samples to derivative samples. Returns the derivative matrix on the interval [a,b] of the given nord. If a function on [a,b] is sampled at nord modified points then multiplication of the sample values by this matrix returns the derivative of the function at the modified points. Can be used to convert DEs to matrix equations. Arguments: nord -- integer, number of modified sample points (nord >= 2) a,b -- interval that maps to [-1,1] Returns: Matrix to convert function values at the modified points to the derivative values at the modified points. Exceptions: ValueError -- nord must be >= 2 """ x = modified_points(nord, a, b) return modified_values(x, nord, 1, a, b) def modified_integral(nord, a=-1.0, b=1) : """Matrix to convert function samples to integral samples. Returns the integral function on [a,b] is sampled at nord modified points then multiplication of the sample values by this matrix returns the integral of the function from a to the sample points. Can be used to convert integral equations to matrix equations. Arguments: nord -- integer, number of modified sample points (nord >= 1) a,b -- interval that maps to [-1,1] Returns: Matrix to convert function values at the modified points to the integral values at the modified points. The integral starts at a. Exceptions: ValueError -- nord must be >= 1 """ x = modified_points(nord, a, b) return modified_values(x, nord, -1, a, b)