How to solve recurrence relations in Python

You're right this can be solved using linear algebra. What I've done below is a simple hard-coded translation. Your equations for p(0) to p(3) are coded up by rearranging them so that the right hand side is =0. For p(4) and p(5) which appear in the recurrence relations as base cases, there is an =1 on the right hand side.

• -p(0) + p(2)/2 = 0

• p(i-1)/2 - p(i) + p(i+2)/2 = 0 for i > 0 and i < x

• p(i) = 1 if i >= x

Here is the program hardcoded for n=4

import numpy
a=numpy.array([[-1,   0, 0.5,  0,   0,   0], # 0
               [0.5, -1,   0,0.5,   0,   0], # 1
               [0,  0.5,  -1,  0, 0.5,   0], # 2
               [0,    0, 0.5, -1,   0, 0.5], # 3
               [0,    0,   0,  0,   1,   0], # 4
               [0,    0,   0,  0,   0,   1], # 5
              ])
b=numpy.array([0,0,0,0,1,1])
# solve ax=b
x = numpy.linalg.solve(a, b)
print x

Edit, here is the code which constructs the matrix programmatically, only tested for n=4!

n = 4

# construct a
diag = [-1]*n + [1]*2
lowdiag = [0.5]*(n-1) + [0]*2
updiag = [0.5]*n
a=numpy.diag(diag) + numpy.diag(lowdiag, -1) + numpy.diag(updiag, 2)

# solve ax=b
b=numpy.array([0]*n + [1]*2)
x = numpy.linalg.solve(a, b)

print a
print x[:n]

This outputs

[[-1.   0.   0.5  0.   0.   0. ]
 [ 0.5 -1.   0.   0.5  0.   0. ]
 [ 0.   0.5 -1.   0.   0.5  0. ]
 [ 0.   0.   0.5 -1.   0.   0.5]
 [ 0.   0.   0.   0.   1.   0. ]
 [ 0.   0.   0.   0.   0.   1. ]]
[ 0.41666667  0.66666667  0.83333333  0.91666667]

which matches the solution in your comment under your question.

The issue here is that you end up in an infinite recursion regardless of where you start, because the recursion isn't explicit, but rather ends up yielding systems of linear equations to solve. If this were a problem you had to solve using Python, I would use Python to calculate the coefficients of this system of equations and use Cramer's rule to solve it.

Edit: Specifically, your unknowns are p(0), ..., p(x-1). One coefficient row vector right off the bat is (1, 0, -1/2, 0, ..., 0) (from p(0)-p(2)/2=0), and all the others are of the form (..., -1/2, 1, 0, -1/2, ...). There are x-1 of these (one for each of p(1), ..., p(x-1)) so the system either has a unique solution or none at all. Intuitively, it seems like there should always be a unique solution.

The two last equations would be unique since they would feature p(x) and p(x+1), so those terms would be ommitted; the column vector for the RHS of Cramer's rule would then be (0, 0, ..., 0, 1/2, 1/2), I believe.

Numpy has matrix support.

This is not an answer to the posted question, but this page is the top Google hit for "solve recurrence relation in Python" so I will write an answer.

If you have a linear recurrence and you want to find the recursive formula, you can use Sympy's find_linear_recurrence function. For example, suppose you have the following sequence: 0, 1, 3, 10, 33, 109, 360, 1189, 3927, 12970. Then the following code produces the recurrence relation:

import sympy
from sympy.abc import n
L = [0, 1, 3, 10, 33, 109, 360, 1189, 3927, 12970]
print(sympy.sequence(L, (n, 1, len(L))).find_linear_recurrence(len(L)))

The output is:

[3, 1]

So you know A(n) = 3*A(n-1) + A(n-2).