diegoaguilar
diegoaguilar

Reputation: 8376

Correct recursive python implementation of Newton 's differences interpolation method, obtaining some of the returned values inside recursion

I wrote a recursion function in python to evaluate the sequence of an interpolation method.

It's graphically explained in this image:

enter image description here

f[x]=f(x) and f[x0,x1]= f[x1]-f[x0]) / (x1 - x0) and so when f[x0,x1,...xn]=f[all_leastFirst,allbutLast] / xlast-xfirst.

This is it then, recursively.

I had got the following code:

xxs=[]
yys=[]
coeficientes = []
h = {}
r = 0.0

def a_j(xx,yy):
    global r
    if len(yy) == 1:
        h[xx[0]] = yy[0]
        return yy[0]
    else:
        r = (a_j(xx[1:],yy[1:])  - a_j(xx[:-1],yy[:-1])) / (xx-1]-xx[0])
        h[''.join(str(i) for i in xx[::-1])]=r
        coeficientes.append(r)
        return ( r )

But it was needed to get as output an array with only the numbers marked in a green circle. I was lost about how to get only those in a recursive implementation. One common pattern about them will be they ALWAYS start at X_0, so I opted about tagging them or using a dictionary might help.

Expected result would be:

[1,1.71828,1.47625,.84553]

I was obtaining:

[1, 2.71828, 7.3890599999999997, 20.085540000000002, 1.71828, 4.6707799999999997, 12.696480000000001, 1.4762499999999998, 4.0128500000000003, 0.84553333333333347]

For another run wit different parameters, if it's called by:

a_j([1,2,3,5][4,3.5,4,5.6])

Should output:

[4,-0.5,0.5,-0.1]

I was obtaining:

[4, 3.5, 4, 5.6, -0.5, 0.5, 0.5, 0.7999999999999998, 0.09999999999999994, -0.10000000000000002]

Another example:

a_j([-2,-1,0,1,2], [13,24,39,65,106])

Will output:

[13, 24, 39, 65, 106, 11, 15, 2, 26, 5, 1, 41, 7, 0, -1]

But the output should be:

[13,11,2,1.167,-0.125]

I also managed to code this iterative implementation, which is already correct:

diferencias = {}
coeficientes = []

def sublists_n(l, n):
    subs = []
    for i in range(len(l)-n+1):
        subs.extend([l[i:i+n]])
    return subs

def sublists(l):
    subs = []
    for i in range(len(l)-1,0,-1):
        subs.extend(sublists_n(l,i))
    subs.insert(0,l)
    return subs[::-1]


def diferenciasDivididas(xx,yy,x):

    combinaciones = sublists([i for i in range(len(xx))])

    for c in combinaciones:

        if len(c) == 1:
            diferencias[str(c[0])]= float(yy[c[0]])
            if c[0] == 0:
                coeficientes.append(float(yy[c[0]]))

        else:
            c1 = diferencias.get(''.join(str(i) for i in c[1:]))
            c2 = diferencias.get(''.join(str(i) for i in c[:-1]))

            d = float(( c1 - c2 ) / ( xx[c[len(c)-1]] - xx[c[0]] ))

            diferencias[''.join(str(i) for i in c)] = d

            if c[0] == 0:
                coeficientes.append(float(d))

I only wonder what was I missing?

Upvotes: 3

Views: 2496

Answers (4)

sanooj
sanooj

Reputation: 493

I have modified the script a bit.

    array=[]
    r='s'
    s=0
    def a_j(xx,yy):
        global r,s
        if r == 's':
            s=xx[0]
            r=0.0
        if len(yy) == 1:
            if xx[0]==s: array.append(yy[0])
            return float(yy[0])

        else:
            r=( a_j(xx[1:],yy[1:])  - a_j(xx[:-1],yy[:-1])) / (xx[-1]-xx[0])
            if xx[0]==s: array.append(r)
            return float(r)

    a_j([1,2,3,5],[4,3.5,4,5.6])
    print array

Output: [4, -0.5, 0.5, -0.10000000000000002]

also, the second example that you have given doesnt look correct. a_j([-2,-1,0,1,2], [13,24,39,65,106]) --> [13,11,4,7,-3]

above answer says that the 3rd element is 4.

    3rd element means --> x(-2,-1,0) -> x(-1,0)  -  x(-2,-1)/(2)
                                     -> x(0)-x(-1)/1  -  x(-1)-x(-2)/(1) /(2)
                                     ->(39-24) - (24-13)   /(2)
                                     ->15-11/(2)
                                     ->4/2 =2

Please correct me if i am wrong.

Upvotes: 1

sanooj
sanooj

Reputation: 493

Try this:

    array=[]
    h={}
    r=0

    def a_j(xx,yy):
        global r
        if len(yy) == 1:
            h[int(xx[0])]=yy[0]
            return yy[0]

        else:
            r=( a_j(xx[1:],yy[1:])  - a_j(xx[:-1],yy[:-1])) / (xx[-1]-xx[0])
            h[int(''.join(str(i) for i in xx[::-1]))]=r
            return r

    a_j([0,1,2,3], [1,2.71828,7.38906,20.08554])
    array=[h[key] for key in  sorted(h.keys())]
    print array

Output: [1, 2.71828, 7.3890599999999997, 20.085540000000002, 1.71828, 4.6707799999999997, 12.696480000000001, 1.4762499999999998, 4.0128500000000003, 0.84553333333333347]

In this code, the values are first assigned to a dict with keys as the elements of xx reversed and converted to an integer.

Upvotes: 1

sanooj
sanooj

Reputation: 493

You are getting negative values here because you have not enclosed the subtraction in parenthesis.Otherwise the code looks good.

   r = ( a_j(xx1,yy1)  - a_j(xx0,yy0)  ) / (xx[len(xx)-1]-xx[0])

http://www.mathcs.emory.edu/~valerie/courses/fall10/155/resources/op_precedence.html

Upvotes: 1

sabbahillel
sabbahillel

Reputation: 4425

Since you are doing recursion, you will append each value as it exits the function, you will wind up getting the appending done in reverse xn, ... , x3, x2, x1

Once you finish the total recursion and exit for the last time, just reverse the list, which is relatively simple by several methods and has been asked before. I leave the method you want to use up to you (or remember "Google is your friend")

Upvotes: 0

Related Questions