Solving a mathematical equation recursively in Python

Question

I want to solve an equation which I am supposed to solve it recursively, I uploaded the picture of formula (Sorry! I did not know how to write mathematical formulas here!) I wrote the code in Python as below:

import math
alambda = 1.0
rho = 0.8
c = 1.0
b = rho * c / alambda
P0 = (1 - (alambda*b))
P1 = (1-(alambda*b))*(math.exp(alambda*b) - 1)

def a(n):
    a_n = math.exp(-alambda*b) * ((alambda*b)**n) / math.factorial(n)
    return a_n

def P(n):
    P(n) = (P0+P1)*a(n) + sigma(n)

def sigma(n):
    j = 2
    result = 0
    while j <= n+1:
        result = result + P(j)*a(n+1-j)
        j += 1
    return result

It is obvious that I could not finish P function. So please help me with this. when n=1 I should extract P2, when n=2 I should extract P3. By the way, P0 and P1 are as written in line 6 and 7. When I call P(5) I want to see P(0), P(1), P(2), P(3), P(4), P(5), P(6) at the output.

Answer 1

You need to reorganize the formula so that you don't have to calculate P(3) to calculate P(2). This is pretty easy to do, by bringing the last term of the summation, P(n+1)a(0), to the left side of the equation and dividing through by a(0). Then you have a formula for P(n+1) in terms of P(m) where m <= n, which is solvable by recursion.

As Bruce mentions, it's best to cache your intermediate results for P(n) by keeping them in a dict so that a) you don't have to recalculate P(2) etc everytime you need it, and b) after you get the value of P(n), you can just print the dict to see all the values of P(m) where m <= n.

import math
a_lambda = 1.0
rho = 0.8
c = 1.0
b = rho * c / a_lambda
p0 = (1 - (a_lambda*b))
p1 = (1-(a_lambda*b))*(math.exp(a_lambda*b) - 1)
p_dict = {0: p0, 1: p1}

def a(n):
    return math.exp(-a_lambda*b) * ((a_lambda*b)**n) / math.factorial(n)

def get_nth_p(n, p_dict):
    # return pre-calculated value if p(n) is already known
    if n in p_dict:
        return p_dict[n]

    # Calculate p(n) using modified formula
    p_n = ((get_nth_p(n-1, p_dict)
            - (get_nth_p(0, p_dict) + get_nth_p(1, p_dict)) * a(n - 1)
            - sum(get_nth_p(j, p_dict) * a(n + 1 - j) for j in xrange(2, n)))
          / a(0))

    # Save computed value into the dict
    p_dict[n] = p_n
    return p_n

get_nth_p(6, p_dict)
print p_dict

Edit 2

Some cosmetic updates to the code - shortening the name and making p_dict a mutable default argument (something I try to use only sparingly) really makes the code much more readable:

import math

# Customary to distinguish variables that are unchanging by making them ALLCAP
A_LAMBDA = 1.0
RHO = 0.8
C = 1.0
B = RHO * C / A_LAMBDA
P0 = (1 - (A_LAMBDA*B))
P1 = (1-(A_LAMBDA*B))*(math.exp(A_LAMBDA*B) - 1)

p_value_cache = {0: P0, 1: P1}

def a(n):
    return math.exp(-A_LAMBDA*B) * ((A_LAMBDA*B)**n) / math.factorial(n)

def p(n, p_dict=p_value_cache):
    # return pre-calculated value if p(n) is already known
    if n in p_dict:
        return p_dict[n]

    # Calculate p(n) using modified formula
    p_n = ((p(n-1)
            - (p(0) + p(1)) * a(n - 1)
            - sum(p(j) * a(n + 1 - j) for j in xrange(2, n)))
          / a(0))

    # Save computed value into the dict
    p_dict[n] = p_n
    return p_n

p(6)
print p_value_cache

Answer 2

You could fix if that way:

import math
alambda = 1.0
rho = 0.8
c = 1.0
b = rho * c / alambda


def a(n):
  # you might want to cache a as well
  a_n = math.exp(-alambda*b) * ((alambda*b)**n) / math.factorial(n)
  return a_n


PCache={0:(1 - (alambda*b)),1:(1-(alambda*b))*(math.exp(alambda*b) - 1)}

def P(n):
 if n in PCache:
   return PCache[n]
 ret= (P(0)+P(1))*a(n) + sigma(n)
 PCache[n]=ret
 return ret

def sigma(n):
  # caching this seems smart as well
  j = 2
  result = 0
  while j <= n+1:
    result = result + P(j)*a(n+1-j)
    j += 1
  return result

void displayP(n):
  P(n) # fill cache :-)
  for x in range(n):
    print ("%u -> %d\n" % (x,PCache[x]))

Instead of managing the cache manually, you might want to use a memoize decorator (see http://www.python-course.eu/python3_memoization.php )

Notes:

• not tested, but you should get the idea behind it
• your recurrence won't work P(n) depends on P(n+1) on your equation... This will never end
• It looks like I misunderstood P0 and P1 as being Both constants (big numbers) and results (small numbers, indices)... Notation is not the best choice I guess...