CODEWITHSUNDEEP
pythonmathnumpyscipy
Simd
Simd

Reputation: 21194

How to solve recurrence relations in Python

I am trying to write code to give numerical answers to a recurrence relation. The relation itself is simple and is defined as follows. The variable x is an integer

This is also in this code.

from __future__ import division
def p(i):
    if (i == 0):
        return p(2)/2
    if (i >= x):
        return 1
    return p(i-1)/2+p(i+2)/2


x = 4
#We would like to print p(0)  for example.

This of course doesn't actually let you compute p(0). How can you do this in python?

Is it possible to set up a system of simultaneous equations which numpy.linalg.solve can then solve?

Upvotes: 7

Views: 12505

Answers (4)

Bai Li
Bai Li

Reputation: 1138

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).

Upvotes: 5

TooTone
TooTone

Reputation: 8126

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.

Upvotes: 7

peterx
peterx

Reputation: 196

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.

Upvotes: 2

Adam Smith
Adam Smith

Reputation: 54183

I'm confused because your code seems like it should do just that.

def p(i):
    x = 4 # your constant should be defined in-function

    if (i == 0):
        return p(2)/2
    elif (i >= x):
        return 1
    return p(i-1)/2+p(i+2)/2

The big problem here is your recursion. For p(1) it does:

p(0)/2                    +                    p(3)/2
p(2)/2                    +           p(2)/2 + p(4)/2
p(1)/2                    +              p(1)/2 + 1/2
# each side of the operator is now the same as the original problem!
# that's a sure sign of infinite recursion.

What do you EXPECT to be the output?

Upvotes: 0

Related Questions