Sage, ore_algebra

In [1]:
import sys
sys.path.append("/home/marc/cur/ore_algebra/src/")
In [5]:
from ore_algebra import *
In [11]:
Pols.<z> = PolynomialRing(QQ)
In [12]:
(z + 1)*(z - 1)
Out[12]:
z^2 - 1
In [13]:
z.degree()
Out[13]:
1
In [14]:
DiffOps.<Dz> = OreAlgebra(Pols)
In [6]:
Dz*z
Out[6]:
z*Dz + 1
In [7]:
dop = (z*Dz - 1)*(z*Dz  - 5)
dop
Out[7]:
z^2*Dz^2 - 5*z*Dz + 5
In [8]:
dop.polynomial_solutions()
Out[8]:
[(z,), (z^5,)]

Toy Examples

In [15]:
(Dz - 1).numerical_solution(ini=[1], path=[0, 1])
Out[15]:
[2.7182818284590452 +/- 3.66e-17]
In [16]:
(Dz - 1).numerical_solution(ini=[1], path=[0, 1], eps=1e-100)
Out[16]:
[2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642743 +/- 4.70e-102]
In [10]:
(Dz - 1).numerical_solution(ini=[1], path=[0, i*pi])
Out[10]:
[-1.0000000000000000 +/- 8.68e-19] + [+/- 1.16e-37]*I
In [11]:
dop = Dz*z*Dz; dop
Out[11]:
z*Dz^2 + Dz
In [12]:
dop(log(z))
Out[12]:
0
In [13]:
dop.numerical_solution(ini=[0,1], path=[1,2])
Out[13]:
[0.69314718055994531 +/- 1.23e-18]
In [15]:
dop.numerical_solution(ini=[0,1], path=[1, i, -1])
Out[15]:
[+/- 4.79e-18] + [3.1415926535897932 +/- 4.35e-17]*I
In [16]:
dop.numerical_solution(ini=[0,1], path=[1, -i, -1])
Out[16]:
[+/- 4.79e-18] + [-3.1415926535897932 +/- 4.35e-17]*I
In [17]:
dop.numerical_solution(ini=[0,1], path=[1, i, -1, -i, 1, i, -1])
Out[17]:
[+/- 2.91e-16] + [9.424777960769380 +/- 5.76e-16]*I

Face-centered cubic lattices

In [18]:
from ore_algebra.analytic.examples import fcc
fcc.dop4
Out[18]:
(9*z^10 + 186*z^9 + 1393*z^8 + 4608*z^7 + 6156*z^6 - 256*z^5 - 7488*z^4 - 4608*z^3)*Dz^4 + (126*z^9 + 2304*z^8 + 15322*z^7 + 46152*z^6 + 61416*z^5 + 15584*z^4 - 39168*z^3 - 27648*z^2)*Dz^3 + (486*z^8 + 7716*z^7 + 44592*z^6 + 119388*z^5 + 151716*z^4 + 66480*z^3 - 31488*z^2 - 32256*z)*Dz^2 + (540*z^7 + 7248*z^6 + 35268*z^5 + 80808*z^4 + 91596*z^3 + 44592*z^2 + 2688*z - 4608)*Dz + 108*z^6 + 1176*z^5 + 4584*z^4 + 8424*z^3 + 7584*z^2 + 3072*z
In [19]:
fcc.dop4.numerical_solution([0, 0, 0, 1], [0, 1], 1e-60)
Out[19]:
[1.1058437979212047601829954708858510744395462366387528583649984 +/- 2.78e-62] + [+/- 2.24e-74]*I
In [ ]:
fcc.dop5
In [ ]:
fcc.dop5.numerical_solution([0, 0, 0, 0, 1, 0], [0, 1/5+i/2, 1], 1e-60)
In [ ]:
fcc.dop6
In [ ]:
ini = [0, 0, 0, 0, 0, 1, 0, 0]
fcc.dop6.numerical_solution(ini, [0, 3/2 + i, 1], 1e-60)

Small-Step Walks

In [22]:
from ore_algebra.analytic.examples import ssw
dop = ssw.dop[4,1,1]; dop
Out[22]:
(64*t^6 + 40*t^5 - 30*t^4 - 5*t^3 + t^2)*Dt^3 + (576*t^5 + 200*t^4 - 252*t^3 - 33*t^2 + 5*t)*Dt^2 + (1152*t^4 + 88*t^3 - 468*t^2 - 48*t + 4)*Dt + 384*t^3 - 72*t^2 - 144*t - 12
In [24]:
terms = oeis("A151331").first_terms()
terms[:5]
Out[24]:
(1, 3, 18, 105, 684)
In [32]:
loc = dop.local_basis_monomials(0); loc
ini = [(terms[m.degree(t)] if m.is_polynomial(t) else 0) for m in loc]
ini
Out[32]:
[0, 0, 1]
In [33]:
import itertools

rts = dop.leading_coefficient().roots(QQbar, multiplicities=False)
rts = [rt for rt in rts if not rt.is_zero()]
rts = sorted(rts, key=abs)
rts
Out[33]:
[0.12500000000000000?, -0.25000000000000000?, 0.50000000000000000?, -1]
In [42]:
loc = dop.local_basis_expansions(rts[0])
In [44]:
loc[0]
Out[44]:
log(t - 1/8) - 80/9*(t - 1/8)*log(t - 1/8) + 688/9*(t - 1/8)^2*log(t - 1/8) - 208/27*(t - 1/8)^2 - 469760/729*(t - 1/8)^3*log(t - 1/8) + 250112/2187*(t - 1/8)^3 + 35202880/6561*(t - 1/8)^4*log(t - 1/8) - 24664480/19683*(t - 1/8)^4 - 290776064/6561*(t - 1/8)^5*log(t - 1/8) + 401448448/32805*(t - 1/8)^5
In [45]:
loc[1]
Out[45]:
1 - 48*(t - 1/8)^2 + 55040/81*(t - 1/8)^3 - 5343040/729*(t - 1/8)^4 + 51897344/729*(t - 1/8)^5
In [46]:
loc[2]
Out[46]:
t - 1/8 - 14*(t - 1/8)^2 + 12064/81*(t - 1/8)^3 - 1041128/729*(t - 1/8)^4 + 9473152/729*(t - 1/8)^5
In [52]:
con = dop.numerical_transition_matrix([0, rts[0]])*vector(ini)
con
Out[52]:
([-0.84882636315677512 +/- 4.54e-18] + [+/- 3.27e-31]*I, [-1.6914020435595060 +/- 5.35e-17] + [2.6666666666666667 +/- 5.65e-17]*I, [16.168933032184458 +/- 3.21e-16] + [-23.703703703703704 +/- 5.02e-16]*I)
In [65]:
B.<n> = AsymptoticRing('QQ^n * n^QQ * log(n)^QQ', QQ)
B.coefficients_of_generating_function(lambda t: log(t-1/8), (1/8,), precision=2)
Out[65]:
-8^n*n^(-1) + O(8^n*n^(-2)*log(n))
In [67]:
N(8/(3*pi))
Out[67]:
0.848826363156775