I’m trying to create a Python algorithm for calculating the linear systems by Gauss. The method is as follows. We draw up a matrix of coefficients, including the constant term. Then give a matrix to triangular form. For this, first along the first column (with index 0), each element is divided by the diagonal (a0,0) (Example – 3.8), calculating an index m, and then recalculate the line 2 of the formula: each element (without intercept member from the last column) minus the product of the element above it (from the zero line) and the index m for the second row. Separately, work it out with a free member column (here, the algorithm is not important).

Following similar actions must be done for the third row elements (but considering that in the first iteration of the elements of the second line transformed above algorithm, and m the coefficient will be considered in the second column: respectively dividing all its elements by the diagonal element of the 2nd row a1 1) (example 1.3).
Question: I calculated the column vector m: m = ([1,000,
1,684,
0,632])

And now we have to work with the second row of the matrix. And here is the difficulty with indexation. Firstly, I can not enumerate the values ​​of m, the type of which float. Second, I ask the wrong indexing elements of the second line (in fact – after zero is the first row)

import numpy as np
matrix = np.array ([[3.8, 6.7, -1.2, 5.2],
          [6.4, 1.3, -2.7, 3.8]
          [2.4, -4.5, 3.5, -0.6]])
def gaussFunc (matrix):
  # Calculation len1 (3) and len2 (4) - are not given here
  # Calculation code m of the zero column:
  for row in range (len1):
    for column in range (len2-3):
      m = matrix [row] [column] / matrix [column] [column]
  elem = row-column values ​​# 1 at zero line put in
                # Variable elem
  for i in range (len (m) -1): # go cycle over a range of three values ​​of m minus
                 # The last third, of len error for float
    while row & lt; (Len1-1): # until the first or second line (in len2 them all
                 # 3):
       while column & lt; (Len2-1): # until the column first, second or third
                      # (Minus column free
                     # Member):
         # Recalculated coefficients of the second (first numpy) lines:
         # The current element - m on this line * in the top element
         # Column (from line 0):
         a = matrix [row] [column] - m [i] * matrix [elem] [column]

Answer 1, Authority 100%

At the end of a link to the jupyter notebook with a more or less complete solver Slough. Plus trick, how to count on pyton almost as fast as Fortran:)

Initial response

If you do not pay attention to a possible division by zero, the ghost triangular form can be written very simply:

def gaussFunc (matrix):
  # Function changes through the matrix side effects
  # If you want to keep the old matrix, copy it np.copy
  for nrow, row in enumerate (matrix):
    # Nrow equal line number
    # Row comprises a row of the matrix itself
    divider = row [nrow] # diagonal element
    # Divided by a diagonal element.
    row / = divider
    # Now we have to subtract the reduced line of all of the underlying lines
    for lower_row in matrix [nrow + 1:]:
      factor = lower_row [nrow] # line element in the column nrow
      lower_row - = factor * row # subtract to get the zero in the column nrow
  # All the rows of the matrix changed, in principle, can not return
  return matrix

The result for your example:

array ([[1., 1.76315789, -0.31578947, 1.36842105],
    [-0. 1., 0.06800211, 0.49657354]
    [0, 0, 1, 0.09309401]])

What is the ambush. During the calculations may be zero to the diagonal.

matrix = np.array ([[1, 0, 0, 1],
          [0, 0, 1, 2],
          [0, 1, 0, 3]])

As far as I remember, before dividing the diagonal element first view all the rows starting with the current, down. Select a string with a maximum value in the current column and rearranged from the current one. After that continue.

Result check.

The Make_identity function takes the matrix obtained by the Gauss method and brings it to one. For this line, it is shifted from the bottom up and deducted from the overlying strings to reset the corresponding speakers.

def make_identity (Matrix):
  # Switching rows in reverse order
  For NROW IN RANGE (LEN (Matrix) -1,0, -1):
    ROW = MATRIX [NROW]
    For Upper_ROW in Matrix [: nrow]:
      Factor = Upper_Row [NROU]
      Upper_ROW - = Factor * Row
  Return Matrix.

Result for your example: make_identity (gaussfunc (NP.Copy (Matrix))

array ([[1., 0., 0., 0., 0.53344344],
    [-0. , 1., 0., 0.49024295],
    [0., 0., 1., 0.09309401]])

Cut out the last column, we get the root string: roots = make_identity (gaussfunc (np.copy (Matrix)) [:, 3]

array ([0.53344344, 0.49024295, 0.09309401])

Multiply the found roots on the original matrix and compare with the last column of the original task: NP.Matmul (Matrix [: 3], roots.t) - Matrix [: 3]

The result of the calculation Array ([0.00000000E + 00, -4.44089210E-16, -2.22044605e-16])

Consequently, the roots are found correctly. What and congratulations.

Update

Gauss method with choosing the main element. Plus Added zero treatment to the diagonal.

Def GausspivotFunc (Matrix):
  For NROW IN RANGE (LEN (MATRIX)):
    # NROW is equal to line number
    # NP.Argmax Returns the line number with a maximum element in a reduced matrix
    # which starts with the nrow row. Therefore, you need to add nrow to the result
    Pivot = Nrow + NP.Argmax (ABS (Matrix [NROU: NROW]))
    IF Pivot! = Nrow:
      # Swap
      # Matrix [Nrow], Matrix [Pivot] = Matrix [Pivot], Matrix [NROU] - does not work.
      # You need to rearrange the rows as written below
      Matrix [[Nrow, Pivot]] = Matrix [[Pivot, Nrow]]]
    ROW = MATRIX [NROW]
    Divider = Row [NROW] # diagonal element
    IF ABS (Divider) & LT; 1E-10:
      # Almost zero on the diagonal. It does not make sense to continue, the result of the invoice is unstable
      Raise Valueerror (F "Matrix is ​​incomprehensible. The maximum element in the column {nrow}: {Divider: .3g}")
    # Delim on a diagonal element.
    ROW / = Divider
    # Now you need to subtract the line of all underlying lines
    For Lower_ROW in Matrix [NROW + 1:]:
      Factor = Lower_Row [NROU] # Row element in the NROW column
      Lower_ROW - = Factor * Row # Remove to get zero in the nrow column
  # lead to a diagonal view
  MAKE_IDENTIY (MATRIX)
  Return Matrix.

In this feature, the main focus is how to rearrange the lines. Simple Formula Matrix [Nrow], Matrix [Pivot] = Matrix [Pivot], Matrix [NROW] does not work. With this assignment to the right there are pointers on the string, and the left – addresses where you need copy values. That is, when first assigned to the nrow line, the row Pivot will be copied, and in the pivot line – the contents of the line nrow – – But it is already rewritten from the line Pivot ! That is, actually the row Pivot will be copied to itself. And instead of rearrangement of two lines there will be a copy of one line.

Matrix [[NROW, PIVOT]] = Matrix [[Pivot, Nrow]] works. And with explicit copying also works: Matrix [NRW], Matrix [Pivot] = Matrix [Pivot], NP.Copy (Matrix [NROU])

Update 2

JUPYTER Notebook with a clerk code Slava

In addition to the actual solver, a comparison with the Si-shnaster numpy.linalg.solve and trick how to speed up the speed of the riser’s account on Pajton is 20 times, which is almost as fast as a purely sulder.

Strictly speaking, the riser in numpy Do not even si-shny, and Forthrovsky LAPACK GESV

Answer 2

Trying to view and rearrange the strings – what’s wrong?

Def Gaussfunc (Matrix):
  For NROW, ROW in Enumerate (Matrix):
    MAX_ELEM = ABS (ROW [NROW])
    Max_row = Nrow.
    For i in Range (NROW + 1, LEN (Matrix)):
      IF ABS (Matrix [i] [Nrow]) & gt; ABS (MAX_ELEM):
        max_elem = Matrix [i] [nrow]
        max_row = i # fix the string
    if max_row! = nrow:
      Matrix [NROW], MATRIX [max_row] = Matrix [max_row], Matrix [NROU]