Determining connection coefficients with parameter initial conditions.
Requirements: ore_algebra package
# Define the rings of differential and shift operators
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)
# Start with the differential equation
diff = 5*z^2*(2*z-1)*(z-3)*Dz^4 + 2*z*(59*z^2-139*z+36)*Dz^3 + 6*(61*z^2-80*z+6)*Dz^2 + 12*(25*z-11)*Dz + 36
# Get recurrence on coefficients
rec = diff.to_S(Shift).primitive_part()
rec
(15*n^2 + 57*n + 36)*Sn^2 + (-35*n^2 - 103*n - 66)*Sn + 10*n^2 + 28*n + 18
# The following gives asymptotic expansions for a basis of solutions to the recurrence
show(rec.generalized_series_solutions(n=3))
# There is a basis with the following series expansions near the origin
print(diff.local_basis_expansions(0, order = 4))
[z^(-1), 1 - 1/2*z^2, z + 11/6*z^2, z^(6/5) + 11/6*z^(11/5) + 187/63*z^(16/5) + 1705/351*z^(21/5)]
# A power series solution starting f0 + f1*z + ... is expressed in this basis with coordinates
var('f0 f1')
ini = [0,f0,f1,0]
ini
[0, f0, f1, 0]
# This generating function potentially has singularities at z = 1/2 and z = 3
diff.leading_coefficient().factor()
(10) * (z - 3) * (z - 1/2) * z^2
# Expansions of a basis near z=1/2
for k in diff.local_basis_expansions(1/2, order = 4):
show(k)
# We can compute the change of basis matrix between the basis with series at the origin and the basis with series at 1/2
# We take a low numeric accuracy here to fit on screen
show(diff.numerical_transition_matrix([0,1/2],1e-1))
# The coefficient of the dominant singular term is f0/2 - 3f1/2
bas = diff.numerical_transition_matrix([0,1/2],1e-10) * vector(ini)
show(bas[0])
# Expansions of a basis near z=3
for k in diff.local_basis_expansions(3, order = 5):
show(k)
# We can compute the change of basis matrix between the basis with series at the origin and the basis with series at 3
# We take a low numeric accuracy here to fit on screen
# Note -- the path for continuation must move off the real line to *avoid* the singularity at z=1/2
show(diff.numerical_transition_matrix([0,I,3],1e-1))
# The coefficient of the singular term is (9/2)(f1-f0)
bas2 = diff.numerical_transition_matrix([0,I,3],1e-10) * vector(ini)
show(bas2[0])