recurrent equations solver

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
# Just for fun ;)
# Usage - if you have recurrent equation u(n) + a_k-1 u(n-1) + ... + a_0 u(n-k) = 0
# and condition that u(0) = b_0, u(1) = b_1, ... , u(k-1) = b_k-1
# then you just put the coefficients 1, a_k-1, a_k-2, ... , a_0 to initial_recurrent_coeffs list
# and put b_0, b_1, ... , b_k-1 to first_values list
# For me, I have an equation 1*u(n) + 2*u(n-1) + 1*u(n-2) = 0, and u(0) = -1, u(1) = -1,
# so I put 1, 2, 1 into initial_recurrent_coeffs and -1, -1 into first_values
initial_recurrent_coeffs = [1, 2, 1]
# Add first k-1 values of recurrent sequence
first_values = [-1, -1]
import sympy
from sympy import Symbol
from sympy import symbols
from sympy import poly
from sympy import Idx, Indexed, Matrix
from sympy.solvers import solve_linear_system
def solve_recurrent(coeffs):
from sympy.abc import x
ex = 0
for i in xrange(len(coeffs)):
ex += coeffs[i]*x**(len(coeffs) - 1 - i)
return ex
def roots_counts(polynom):
roots = polynom.all_roots()
counts = {}
for i in xrange(len(roots)):
if counts.has_key(roots[i]):
counts[roots[i]] += 1
else:
counts[roots[i]] = 1
return counts
def solve_linear(counts, first_values):
pow_values = counts.values()
# Warning! I didn't check that str(pow_values[i]) is right choice for cc making
cc = [list(symbols('C'+str(i)+':'+str(i+1)+':'+str(pow_values[i]))) for i in xrange(len(pow_values))]
counts_list = list(counts)
print 'counts list: ', counts_list
under_matrix = [[] for i in xrange(sum(counts.values()))]
for power_of_roots, ex in enumerate(under_matrix):
for i_idx, root_value in enumerate(counts):
for j_idx in xrange(counts[root_value]):
under_matrix[power_of_roots].append(power_of_roots**j_idx * root_value**power_of_roots)
for i, item in enumerate(under_matrix):
item.append(first_values[i])
linear_system = Matrix(under_matrix)
all_variables = [elem for item in cc for elem in item]
solution = solve_linear_system(linear_system, *all_variables)
# Build recurrent function
result_expression = 0
x = Symbol('x')
print cc
root_values = list(counts.keys())
for root_idx, root_coeffs in enumerate(cc):
subresult_expression = 0
for coeff_idx, coeff_symbol in enumerate(root_coeffs):
subresult_expression += solution[cc[root_idx][coeff_idx]] * x**coeff_idx
subresult_expression *= root_values[root_idx]**x
result_expression += subresult_expression
print result_expression
print sympy.simplify(result_expression)
return sympy.lambdify(x, result_expression)
def naive_recursive_equation(n):
if n < 0: return 0
if n == 0: return -2
if n == 1: return 1
return 3*naive_recursive_equation(n-1) - 2*naive_recursive_equation(n-2)
def naive_recursive_equation_b(n):
if n < 0: return 0
if n == 0: return -1
if n == 1: return -1
return -2*naive_recursive_equation_b(n-1) - naive_recursive_equation_b(n-2)
# Make polynom from recurrent coefficients
characteristic_polynom = sympy.poly_from_expr(solve_recurrent(initial_recurrent_coeffs))[0]
# find all roots of polynom with powers of each root
counts = roots_counts(characteristic_polynom)
# solve recurrent equations and give all Cij, that exactly describes solution
recurrent_generator = solve_linear(counts, first_values)
for i in xrange(20):
print recurrent_generator(i), naive_recursive_equation_b(i)
# Nice job, my friend 8-)