In [1]:
%pylab inline
from sympy.utilities.lambdify import lambdify
from sympy import init_session
init_session()
from numpy import arange, abs
from pandas import DataFrame
from pylab import plot
In [2]:
p1 = (1/2)*((x + sqrt(x**2 - 1))**m + (x - sqrt(x**2 - 1))**m)
p2 = cos(m*acos(x))
Eq(p1, p2)
Out[2]:
In [3]:
result = []
powers = range(1,10)
points = arange(-2, 2, 0.25)
def cleverLog(x):
if x == 0: return '-∞'
else: return log(x)
for xv in points:
poly1 = p1.subs(x, xv)
poly2 = p2.subs(x, xv)
result.append([cleverLog(abs(N(poly1.subs(m, mv) - poly2.subs(m, mv), 3))) for mv in powers])
DataFrame(data=result, index=points, columns=powers)
Out[3]:
In [4]:
def drawChebyshev():
xs = pylab.linspace(-1, 1, 1000)
for i in range(3,10,2):
cheb = lambdify(x, p2.subs(m, i), "numpy")
plot(xs, cheb(xs))
title(p2.subs(m, i))
show()
drawChebyshev()
In [9]:
def interpolate(f, points):
values = f(points)
M = pylab.ma.vander(points)
invM = pylab.linalg.inv(M)
a = dot(invM, values)
return Poly(a, x).as_expr()
def absOmega(points):
pcopy = points
return pylab.vectorize(lambda x: abs(pylab.prod([x - p for p in pcopy])))
def chebspace(start, end, n):
cos = pylab.cos
pi = pylab.pi
return [(start+end)/2 + ((end-start)/2)*cos(pi*(2*k+1)/(2*n)) for k in range(0,n)]
def maximum(f, start, end):
from scipy.optimize import minimize_scalar
return f(minimize_scalar(lambda x: -f(x), bounds=(start, end), method='bounded').x)
fun = lambdify(x, tanh(x), 'numpy')
dfun = lambdify(x, tanh(x).diff(x), 'numpy')
Ifun = lambdify(x, log(cosh(x)), 'numpy')
In [6]:
def graphics():
n = 20
points = pylab.linspace(-2, 2, n)
chebPoints = chebspace(-2, 2, n)
L1 = interpolate(fun, points)
L1_f = lambdify(x, L1, 'numpy')
L2 = interpolate(fun, chebPoints)
L2_f = lambdify(x, L2, 'numpy')
dL1_f = lambdify(x, L1.diff(x), 'numpy')
dL2_f = lambdify(x, L2.diff(x), 'numpy')
IL1_f = lambdify(x, L1.integrate(x), 'numpy')
IL2_f = lambdify(x, L2.integrate(x), 'numpy')
t = pylab.linspace(-2, 2, 1000)
plot(t, fun(t))
title(r'$\tanh(x)$')
show()
plot(t, L1_f(t))
title(r'$L_1(x)$')
show()
plot(t, L2_f(t))
title(r'$L_2(x)$')
show()
plot(t, pylab.vectorize(lambda x: abs(fun(x) - L1_f(x)))(t))
title(r'$|\Delta_1|$')
show()
plot(t, pylab.vectorize(lambda x: abs(fun(x) - L2_f(x)))(t))
title(r'$|\Delta_2|$')
show()
plot(t, dfun(t))
title(r"$\tanh'(x)$")
show()
plot(t, dL1_f(t))
title(r'$dL_1$')
show()
plot(t, dL2_f(t))
title(r'$dL_2$')
show()
plot(t, pylab.vectorize(lambda x: abs(dfun(x) - dL1_f(x)))(t))
title(r'$|d\Delta_1|$')
show()
plot(t, pylab.vectorize(lambda x: abs(dfun(x) - dL2_f(x)))(t))
title(r'$|d\Delta_2|$')
show()
plot(t, Ifun(t))
title(r"$∫\tanh(x)$")
show()
plot(t, IL1_f(t))
title(r'$∫L_1$')
show()
plot(t, IL2_f(t))
title(r'$∫L_2$')
show()
plot(t, pylab.vectorize(lambda x: abs(Ifun(x) - IL1_f(x)))(t))
title(r'$|∫\Delta_1|$')
show()
plot(t, pylab.vectorize(lambda x: abs(Ifun(x) - IL2_f(x)))(t))
title(r'$|∫\Delta_2|$')
show()
graphics()
In [7]:
def errors(L, start, end):
L_f = lambdify(x, L, 'numpy')
err = maximum(lambda x: abs(fun(x) - L_f(x)), start, end)
dL = lambdify(x, diff(L, x), 'numpy')
errd = maximum(lambda x: abs(dfun(x) - dL(x)), start, end)
IL = lambdify(x, integrate(L, x), 'numpy')
errI = maximum(lambda x: abs(Ifun(x) - IL(x)), start, end)
return (err, errd, errI)
def approximate(n, start, end):
unifiedPoints = pylab.linspace(start, end, n)
chebyshevPoints = chebspace(start, end, n)
L1 = interpolate(fun, unifiedPoints)
Ω1 = absOmega(unifiedPoints)
err1, errd1, errI1 = errors(L1, start, end)
L2 = interpolate(fun, chebyshevPoints)
Ω2 = absOmega(chebyshevPoints)
err2, errd2, errI2 = errors(L2, start, end)
teorRel = maximum(Ω2, start, end)/maximum(Ω1, start, end)
rel = abs(err2/err1)
return (err1, err2, errd1, errd2, errI1, errI2, teorRel, rel)
In [10]:
results = []
for n in range(5,51,5):
results.append(approximate(n, -2, 2))
DataFrame(results, index=range(5,51,5), columns=('Δ1','Δ2', 'dΔ1', 'dΔ2', '∫Δ1', '∫Δ2', '|Ω2|/|Ω1|', 'Δ2/Δ1'))
Out[10]: