Why is this matrix not working in my program?

[etotheipi]

https://projecteuler.net/problem=101

import numpy as np
for j in range (1,11):
    M = np.empty([j, j])
    for x in range(1,j+1):
        for y in range(1,j+1):
            M[y,x] = y**(j-x)
    Minv = np.linalg.inv(M)

The $j^{\mathrm{th}}$ estimate $\mathrm{OP}(j,n)$ which fits $j$ data points is a $(j-1)^{\mathrm{th}}$ degree polynomial $u^{(j)}(n) = a_{j-1} n^{j-1} + \dots + a_{0}$ which agrees with $u_n$ on the set $n \in \{1, \dots, j \}$, so$$\mathbf{a} := \begin{pmatrix} a_{j-1} \\ \vdots \\ a_0 \end{pmatrix} = \begin{pmatrix} 1 & 1 & \cdots & 1 \\ 2^{j-1} & 2^{j-2} & \cdots & 1 \\ \vdots & \vdots & & \vdots \\ j^{j-1} & j^{j-2} & \cdots & 1\end{pmatrix}^{-1} \begin{pmatrix} u_1 \\ \vdots \\ u_j \end{pmatrix} := \mathbf{M}^{-1} \mathbf{u}$$I need to make the matrix $\mathbf{M}$, why does it not work?

M[y,x] = y**(j-x)
IndexError: index 1 is out of bounds for axis 0 with size 1

 

[jedishrfu]

For the M[x,y] assignment do you mean M[x-1,y-1] ?

 

[etotheipi]

oh yeah that's it thanks

 

[etotheipi]

It does not give me the correct answer

import numpy as np
sum = 0
def p(k):
    s = 0
    for i in range(0,11):
        s += ((-1)**i)*(k**i)
    return s
for j in range (1,11):
    M = np.empty(shape=(j,j))
    for x in range(1,j+1):
        for y in range(1,j+1):
            M[y-1,x-1] = y**(j-x)
    Minv = np.linalg.inv(M)
    U = np.empty(shape=(j,1))
    for i in range(1, j+1):
        U[i-1] = p(i)
    A = np.matmul(Minv, U)
    z = 1
    FIT_found = False
    while(FIT_found == False):
        estimate = 0
        for r in range(0,j):
            estimate += A[r]*z**r
        if(estimate != p(z)):
            sum += estimate
            FIT_found = True
        else:
            z = z + 1      
print(sum)

 

[pbuk]

It does not give me the correct answer

I think the idea of doing Project Euler problems is to work on your code until it does

When posting code here it helps if you add the language so syntax highlighting works:

import numpy as np
sum = 0
def p(k):
    s = 0
    for i in range(0,11):
        s += ((-1)**i)*(k**i)
    return s
for j in range (1,11):
    M = np.empty(shape=(j,j))
    for x in range(1,j+1):
        for y in range(1,j+1):
            M[y-1,x-1] = y**(j-x)
    Minv = np.linalg.inv(M)
    U = np.empty(shape=(j,1))
    for i in range(1, j+1):
        U[i-1] = p(i)
    A = np.matmul(Minv, U)
    z = 1
    FIT_found = False
    while(FIT_found == False):
        estimate = 0
        for r in range(0,j):
            estimate += A[r]*z**r
        if(estimate != p(z)):
            sum += estimate
            FIT_found = True
        else:
            z = z + 1
print(sum)

 

[pbuk]

I always get confused by Python's range() function: note that range(1, 11) gives you [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]: this is the first thing I'd check.

And then i'd check that my j, x, and y values are doing what I want in the inner loop (either with an IDE or a simple print statement).

Edit: first of all I'd refactor it so all my loops start at 0, it makes it much easier to keep track.

 

[etotheipi]

thanks okay I'm going to finish this tomorrow because I have homework and my stupid code is giving me -670.99996864

 

[pbuk]

thanks okay I'm going to finish this tomorrow because I have homework and my stupid code is giving me -670.99996864

yeah, I don't think that's right.

 

[jedishrfu]

I always get confused by Python's range() function: note that range(1, 11) gives you [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]: this is the first thing I'd check.

And then i'd check that my j, x, and y values are doing what I want in the inner loop (either with an IDE or a simple print statement).

Edit: first of all I'd refactor it so all my loops start at 0, it makes it much easier to keep track.

That’s how I checked the code. I called that seldom used secret programmers function:

print(‘j = %d, x = %d, y = %d’%(j,x,y))

just before his M assignment.

 

[pbuk]

I'd look for a simpler method of polynomial interpolation if I were you.