2

Consider this difference equation:

    

The solution is

    

I am trying to solve it numerically in python, to explicate issues that arise with floating point computations.

I wrote a function that computes xn+1

def diff(n):
    c = 1
    b = -1/5.0
    a =  0
    for i in xrange(n):
        a = 14/5.0*b+3/5.0*c
        b, c = a, b
    return a

but I don't know how to solve this numerically and then to explain why python can not provide the xn = (-1/5)n solution.

I can see that for larger n, the return value of my function diverges from the true value.

Bob__
  • 12,361
  • 3
  • 28
  • 42
Tinatim
  • 143
  • 6
  • 2
    Please note that MathJax isn't available in SO, you should rewrite your formulas (unicode art), post them as images or try this: https://stackoverflow.com/a/47798853/4944425 . – Bob__ Jul 09 '18 at 08:53
  • 1
    Should there be an `x_{n-1}` somewhere in the original formula? Otherwise you just have `x_{n+1} = 17/5 x_n`, which doesn't match the solution you give. – Mark Dickinson Jul 10 '18 at 18:51
  • @MarkDickinson Given the python code, I guess the last term (3/5), but the OP should clarify (I didn't change the formulas, with my edit). – Bob__ Jul 10 '18 at 20:28

1 Answers1

0

The second term in the equation should be x_{n-1} . As you know floating point errors propagate and lead to errors. One quick solution is to use higher precision arithmetic like so:

from decimal import *

getcontext().prec = 100

def diff(n):
    c = Decimal(1)
    b = Decimal(-1)/Decimal(5)
    a =  Decimal(0)
    coeff1 = Decimal(14) / Decimal(5)
    coeff2 = Decimal(3) / Decimal(5)
    for i in range(n):
        a = coeff1 * b + coeff2 * c
        b, c = a, b
    return a

Which should be good for the first 100 values. This doesn't solve the issue of how to do this with floating point precision. Most likely some transformation of the problem (such as taking the exponent of both sides) might allow you to do this. However finding the right transformation is tricky.

jman
  • 685
  • 5
  • 15