# -*- coding: utf-8 -*- """ Created on Sat Jan 14 17:05:15 2017 @author: Mariusz """ import numpy as np import math from chebyshev import chebyshev from nest import nest from newtdd import newtdd import matplotlib.pyplot as plt def myexp(x): try: y = math.exp(x) except OverflowError: y = float('inf') return y def discretize(x, edges): ''' discretize group numeric data into bins or categories. BINS = discretize(X, EDGES) returns the indices of the bins that the elements of X fall into. EDGES is a numeric vector that contains bin edges in monotonically increasing order. An element X[i] falls into the j-th bin if EDGES[j] <= X[i] < EDGES[j+1], for 0 <= j < N-1 where N is the number of bins and len(EDGES) = N. The last bin includes the right edge such that it contains EDGES[N-2] <= X[i] <= EDGES[N-1]. For out-of-range values where X[i] < EDGES[0] or X[i] > EDGES[N-1] or np.isnan(X[i]), BINS[i] returns float(nan). ''' bins = float('nan')*np.ones(np.shape(x)) for i in range(len(x)): J = np.where(edges <= x[i])[0] if J.size == 0: continue elif (J[-1] == len(edges) and edges(J[-1]) < x[i]): continue else: bins[i] = J[-1] return bins def ln(x): ''' ln(x) function that is discretized with Chebyshev nodes and interpolated ''' assert(all(1.0e-4 <= x)) assert(all(x <= 1.0e4)) # find k such that e^k <= x < e^k+1 # note: exp(-10) ~ 4.54E-05 # exp(10) ~ 2.20E+04 powers = np.arange(-10, 10+1, 1) edges = list(myexp(z) for z in powers) k = discretize(x, edges).astype(int) x2 = np.zeros(np.shape(x)) for i in range(len(x)): x2[i] = x[i]/edges[k[i]] ''' now x2 is in the fundamental domain [1,e] build the approximating polynomial based on N Chebyshev nodes from Exersise 3.3.7 we know the worst-case error formula: error <= (e-1)^N/N/2^(2N-1) Thus N = 25 guarantees error < 10^(-10) ''' N = 25 Xb = chebyshev(myexp(0), myexp(1), N) Yb = np.log(Xb) c = newtdd(Xb, Yb) y = nest(c, x2, Xb) + powers[k] y = np.reshape(y, np.shape(x)) return y, k # Test of function ln(x) if __name__ == '__main__': # 10^-4 <= x <= 10^4 x = np.linspace(1.0e-4, 1.0e4, 200) y, k = ln(x) y_real = np.log(x) error = np.abs(y-y_real) print(error) plt.plot(x,y,'r') plt.plot(x,y_real,'b') plt.show()