0

Is there a simply way of solving large arrays of polynomial equations in Python?

I have to solve 16 equations that look like:

W1*(P1^0) + W2*(P2^0) + ... + W8*(P8^0)     = 2/1
W1*(P1^1) + W2*(P2^1) + ... + W8*(P8^1)     = 0
W1*(P1^2) + W2*(P2^2) + ... + W8*(P8^2)     = 2/3
...
W1*(P1^14) + W2*(P2^14) + ... + W8*(P8^14)  = 2/15
W1*(P1^15) + W2*(P2^15) + ... + W8*(P8^15)  = 0

My idea was to divide this into 2 arrays, J (left part) and I (right side). With a for loop I managed to make an array of 2/(n+1) if n is an even number, with n starting at 0. For the left part I generated an array with indexing Wi and Pi.

I've tried several different solvers, but i'm not managing to succeed with any.

Also, I'm willing to learn a better way of solving this problem if there are any ideas.

from sympy import *
from sympy.solvers.solveset import linsolve
init_printing(use_unicode=True)
   
#definição dos
w1, w2, w3, w4, w5, w6, w7, w8 = symbols('w1:9')
p1, p2, p3, p4, p5, p6, p7, p8 = symbols('p1:9')

#vetor de coeficientes do polinómio I de n=15
I=[]
for i in range(16):
    if i%2==0:
        I.append(2/(i+1))
    else:
        I.append('0')

#vetor de coeficientes do polinómio J de n=15
J=[]

W = IndexedBase('W')
P = IndexedBase('P')

i = symbols('i', cls=Idx)
s=0
J=[]
for n in range(16):
    for i in range(1,9):
        s += W[i]*P[i]**n
    J.append(s)
    s=0

lst_PW = []
for i in range(1,9):
    lst_PW.append(W[i])
    lst_PW.append(P[i])

solve(J-I,lst_PW)
Luis Corte-Real
  • 55
  • 5
  • 1
    I'm not familiar with sympy so maybe this is dumb, but why can't you invert the matrix of coefficients? Assuming det != 0, unique solution exists etc – gph Dec 07 '20 at 00:00
  • Maybe this question is more appropriate for math.stackexchange.com ? Not asking for a programmatic solution, but for a mathematical approach to tackle these equations? – JohanC Dec 07 '20 at 11:23

1 Answers1

0

15th degree polynomials look too heavy to be solved by sympy or any related software.

Apart from that, the code makes a few mistakes:

  • Instead of the string '0', the numberic value 0 should be appended to I.
  • In the given code, the variables W[i] aren't related to the variables w1, w2, w3, ...
  • In Python, you can't use J - I with two lists. For element-wise subtraction, a loop is needed.
  • Writing 2/(i+1) results in a float value, which hampers sympy's attemps for exact symbolic solutions. Writing it as 2/S(i+1) results in a sympy fraction. See Python numbers vs sympy numbers.

If also Python's list comprehensions would be used, the code could look as follows:

from sympy import init_printing, symbols, S, Eq, solve

init_printing(use_unicode=True)

w1, w2, w3, w4, w5, w6, w7, w8 = W = symbols('w1:9')
p1, p2, p3, p4, p5, p6, p7, p8 = P = symbols('p1:9')

I = [2 / S(i + 1) if i % 2 == 0 else 0 for i in range(16)]
J = [sum([W[i] * P[i] ** n for i in range(8)]) for n in range(16)]

solve([Eq(i, j) for i, j in zip(I, J)], list(P) + list(W))

But, unfortunately, sympy can't find a solution.

JohanC
  • 71,591
  • 8
  • 33
  • 66
  • I agree that it is unlikely that a general purpose computational polynomial solver could find a solution to this. The system is linear in some variables but still there might be as many as `8**16 == 281474976710656` solutions *if* the solution set is finite. The system does have a lot of structure though so it might be possible to make progress by reasoning rather than computation: this might be a better question for math.stackexchange. It might also be possible to find some solutions numerically. – Oscar Benjamin Dec 07 '20 at 00:47
  • @OscarBenjamin @JohanC That's a nice way of rewriting it. I am thinking that if the unknowns were just `W`, and if `P` were known, then the resulting matrix would be of [Vandermonde type](https://en.wikipedia.org/wiki/Vandermonde_matrix), and so [the linear system](https://i.imgur.com/bmfolSj.png) would be easy to solve. – wsdookadr Dec 07 '20 at 09:28
  • And if it can be reformulated as Vandermonde, then one can [solve it like so](https://math.stackexchange.com/a/384167/68328). But I think the constant term column is also supposed to be a list of consecutive increasing powers which it's not. Maybe by way of some variable changes it could be brought to be a list of powers though ? – wsdookadr Dec 07 '20 at 09:43
  • I forgot to mention something, which can be quite relevante. The variables Pi and Wi are the values from Gauss Quadrature Points and Weights. – Luis Corte-Real Dec 07 '20 at 10:53
  • Doesn't that mean that Pi and Wi can be calculated as the roots of Legendre polynomials? – Oscar Benjamin Dec 07 '20 at 12:31