Search code examples
pythonnumpymathematical-optimizationpolynomials

Tolerances, linalg.solv, polynom solve


I have following problem: I try to solve the equilation by using linalg.solv and it seems to work. But if i try to check it by inserting the aquired coefficients and one of the required points i get a difference of about 30% to the original data. Have i done a mistake, which i dont get? Or do i have do use a different methode to get more accurate datasets. If yes, which one?

Further if i use different values from which i entered while calculating the coefficients, i get strangly high results

data = np.genfromtxt("data1.csv",dtype=float,delimiter=";")

to = data[0,2:]
tc = data[1:,0]
y = data[1:,2:]



a = np.array([[1, to[0], tc[0], to[0]**2, to[0]*tc[0], tc[0]**2, to[0]**3, tc[0]*to[0]**2, to[0]*tc[0]**2, tc[0]**3],
            [1, to[1], tc[1], to[1]**2, to[1]*tc[1], tc[1]**2, to[1]**3, tc[1]*to[1]**2, to[1]*tc[1]**2, tc[1]**3],
            [1, to[2], tc[2], to[2]**2, to[2]*tc[2], tc[2]**2, to[2]**3, tc[2]*to[2]**2, to[2]*tc[2]**2, tc[2]**3],
            [1, to[3], tc[3], to[3]**2, to[3]*tc[3], tc[3]**2, to[3]**3, tc[3]*to[3]**2, to[3]*tc[3]**2, tc[3]**3],
            [1, to[4], tc[4], to[4]**2, to[4]*tc[4], tc[4]**2, to[4]**3, tc[4]*to[4]**2, to[4]*tc[4]**2, tc[4]**3],
            [1, to[5], tc[5], to[5]**2, to[5]*tc[5], tc[5]**2, to[5]**3, tc[5]*to[5]**2, to[5]*tc[5]**2, tc[5]**3],
            [1, to[6], tc[6], to[6]**2, to[6]*tc[6], tc[6]**2, to[6]**3, tc[6]*to[6]**2, to[6]*tc[6]**2, tc[6]**3],
            [1, to[7], tc[7], to[7]**2, to[7]*tc[7], tc[7]**2, to[7]**3, tc[7]*to[7]**2, to[7]*tc[7]**2, tc[7]**3],
            [1, to[8], tc[8], to[8]**2, to[8]*tc[8], tc[8]**2, to[8]**3, tc[8]*to[8]**2, to[8]*tc[8]**2, tc[8]**3],
            [1, to[9], tc[9], to[9]**2, to[9]*tc[9], tc[9]**2, to[9]**3, tc[9]*to[9]**2, to[9]*tc[9]**2, tc[9]**3]])

b = np.array([y[0,0],y[1,1],y[2,2],y[3,3],y[4,4],y[5,5],y[6,6],y[7,7],y[8,8],y[9,9]])
c = np.linalg.solve(a, b)



ges_to = 10
ges_tc = 35


ges_y  =  c[0] + c[1]*ges_to + c[2]*ges_tc + c[3]*ges_to**2 + c[4]*ges_to*ges_tc + c[5]*ges_tc**2 + c[6]*ges_to**3 \
          + c[7]*ges_tc*ges_to**2 + c[8]*ges_to*ges_tc**2 + c[9]*ges_tc**3

Here are the values I use to calculate

('to:', array([ 15.,  10.,   5.,   0.,  -5., -10., -15., -20., -25., -30., -35.]))
('tc:', array([ 30.,  35.,  40.,  45.,  50.,  55.,  60.,  65.,  70.,  80.,  90.]))

('b', array([ 24.,  31.,  35.,  36.,  35.,  33.,  30.,  25.,  21.,  18.]))

('y:', array([[ 24.,  26.,  27.,  27.,  26.,  25.,  23.,  20.,  18.,  15.,  13.],
       [ 30.,  31.,  31.,  30.,  29.,  27.,  24.,  21.,  18.,  16.,  14.],
       [ 35.,  35.,  35.,  33.,  31.,  29.,  26.,  22.,  19.,  16.,  15.],
       [ 40.,  40.,  38.,  36.,  33.,  30.,  27.,  23.,  20.,  16.,  15.],
       [ 45.,  44.,  41.,  39.,  35.,  32.,  28.,  24.,  20.,  17.,  16.],
       [ 49.,  47.,  44.,  41.,  37.,  33.,  29.,  25.,  20.,  17.,  16.],
       [ 53.,  51.,  47.,  43.,  39.,  34.,  30.,  25.,  21.,  17.,  16.],
       [ 57.,  54.,  50.,  45.,  40.,  35.,  30.,  25.,  21.,  17.,  16.],
       [ 61.,  57.,  52.,  47.,  41.,  36.,  31.,  26.,  21.,  17.,  16.],
       [ 64.,  60.,  54.,  59.,  53.,  37.,  32.,  27.,  22.,  18.,  19.],
       [ 67.,  63.,  56.,  61.,  55.,  59.,  34.,  29.,  24.,  18.,  19.]]))

('ges_y:', 49.0625)

Solution

  • The floating point arithmetic is leading you astray. If you look at the determinant of that matrix a, it's something incredibly small like 1.551864434916621e-51. If you compute the determinate with the entries as integers (and avoid floating point arithmetic weirdness) you'll see it's actually 0, and the rank of your matrix is 5. So it's singular, and in general equations like ax = b may not have any solution.

    Another quick thing you can do to see this is np.dot(a, np.linalg.inv(a)) is nowhere close to the identity matrix. Similarly, np.dot(a, c) is nowhere close to b.

    There may or may not be an actual solution to ax = b, but np.linalg.lstsq(a,b) will get you an approximate solution in either case, if that's sufficient for you rneeds.