In [1]:
from sympy import *
from sympy.abc import n, i, N, x, lamda, phi, z, j, r, k, a, alpha
from commons import *
from matrix_functions import *
from sequences import *
import functions_catalog
init_printing()
In [2]:
m=8
In [3]:
C = define(let=Symbol(r'\mathcal{{C}}_{{ {} }}'.format(m)),
be=Matrix(m, m, lambda n,k: (k+1)*binomial(2*n-k, n-k)/(n+1) if n > 0 else int(not k)))
C
Out[3]:
In [4]:
eigendata = spectrum(C)
eigendata
Out[4]:
In [5]:
data, eigenvals, multiplicities = eigendata.rhs
In [6]:
Phi_poly = Phi_poly_ctor(deg=m-1)
Phi_poly
Out[6]:
In [7]:
Phi_polynomials = component_polynomials(eigendata, early_eigenvals_subs=True)
Phi_polynomials
Out[7]:
In [8]:
cmatrices = component_matrices(C, Phi_polynomials)
list(cmatrices.values())
Out[8]:
In [9]:
f_power, g_power, G_power = functions_catalog.power(eigendata, Phi_polynomials)
In [13]:
C_power = G_power(C)
C_power
Out[13]:
In [15]:
define(C_power.lhs, C_power.rhs.applyfunc(factor)) # factored
Out[15]:
In [16]:
assert C_power.rhs == (C.rhs**r).applyfunc(simplify)
In [17]:
inspect(C_power.rhs)
Out[17]:
In [18]:
production_matrix(C_power.rhs).applyfunc(factor)
Out[18]:
In [19]:
f_inverse, g_inverse, G_inverse = functions_catalog.inverse(eigendata, Phi_polynomials)
In [20]:
C_inverse = G_inverse(C)
C_inverse, G_inverse(C_inverse)
Out[20]:
In [21]:
inspect(C_inverse.rhs)
Out[21]:
In [22]:
production_matrix(C_inverse.rhs)
Out[22]:
In [23]:
assert C_inverse.rhs*C.rhs == Matrix(m, m, identity_matrix())
assert C_inverse.rhs == C_power.rhs.subs({r:-1})
In [24]:
f_sqrt, g_sqrt, G_sqrt = functions_catalog.square_root(eigendata, Phi_polynomials)
In [25]:
C_sqrt = G_sqrt(C)
C_sqrt
Out[25]:
In [26]:
inspect(C_sqrt.rhs)
Out[26]:
In [27]:
production_matrix(C_sqrt.rhs)
Out[27]:
In [28]:
assert C.rhs**(S(1)/2) == C_sqrt.rhs
assert C_sqrt.rhs == C_power.rhs.subs({r:S(1)/2})
In [29]:
f_exp, g_exp, G_exp = functions_catalog.exp(eigendata, Phi_polynomials)
In [30]:
C_exp = G_exp(C)
C_exp
Out[30]:
In [35]:
define(C_exp.lhs, C_exp.rhs.applyfunc(factor))
Out[35]:
In [37]:
C_exp1 = define(let=Subs(C_exp.lhs, alpha, 1), be=C_exp.rhs.subs({alpha:1}))
C_exp1
Out[37]:
In [38]:
inspect(C_exp.rhs)
Out[38]:
In [58]:
inspect(C_exp1.rhs)
Out[58]:
In [42]:
eigendata_Cexpt = spectrum(C_exp1)
eigendata_Cexpt
Out[42]:
In [43]:
Phi_polynomials_Cexpt = component_polynomials(eigendata_Cexpt, early_eigenvals_subs=True)
Phi_polynomials_Cexpt
Out[43]:
In [44]:
f_log, g_log, G_log = functions_catalog.log(eigendata, Phi_polynomials)
In [45]:
C_log = G_log(C)
C_log
Out[45]:
In [46]:
inspect(C_log.rhs[1:,:-1])
Out[46]:
In [48]:
production_matrix(C_log.rhs[1:,:-1])
Out[48]:
In [49]:
g_log_Cexpt = Hermite_interpolation_polynomial(f_log, eigendata_Cexpt, Phi_polynomials_Cexpt)
g_log_Cexpt
Out[49]:
In [50]:
g_log_Cexpt = g_log_Cexpt.subs(eigendata_Cexpt.rhs[1])
g_log_Cexpt
Out[50]:
In [52]:
with lift_to_matrix_function(g_log_Cexpt) as G_log_Cexpt:
CC = G_log_Cexpt(C_exp1)
CC
Out[52]:
In [53]:
f_sin, g_sin, G_sin = functions_catalog.sin(eigendata, Phi_polynomials)
In [54]:
C_sin = G_sin(C)
C_sin
Out[54]:
In [ ]:
inspect(C_sin.rhs) # takes long to evaluate
In [55]:
f_cos, g_cos, G_cos = functions_catalog.cos(eigendata, Phi_polynomials)
In [56]:
C_cos = G_cos(C)
C_cos
Out[56]:
In [ ]:
inspect(C_sin.rhs) # takes long to evaluate
In [57]:
assert (C_sin.rhs**2 + C_cos.rhs**2).applyfunc(trigsimp) == Matrix(m,m, identity_matrix())
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.