Как подготовить систему уравнений для python

Насколько я понимаю, вы должны настроить матрицу вашей системы, чтобы позаботиться о RyA и других переменных, которые в настоящее время находятся в правой части. Вы можете сделать это вручную (в этом случае этот вопрос выходит за рамки этого сайта, это чисто математическое упражнение) или использовать, например. sympy вместо np.linalg.solve(), который может решить алгебраическую часть задачи за вас:

from sympy import Matrix, symbols, solve

x1, x2, f1 = symbols('x1 x2 f1')
X = Matrix([0, x1, x2])
B = Matrix([f1, 50, 60])
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

solve(M * X - B, [x1, x2, f1])

# {f1: 40, x2: 100/3, x1: -30}

Обратите внимание, что sympy может медленнее решать числовые линейные системы, чем numpy.linalg, поэтому вы можете использовать sympy для выполнения алгебраической части работы, вычисления матрицы и правой части, а затем использовать numpy.linalg.solve для ее решения.

import numpy as np
from sympy import expand
def symbolic_to_matrix(F, variables):

    """
    F is a symbolic vector function that is a left hand side of equation F = 0
    variables is a list of variables (sympy.Symbol's) which F depends on.

    Assuming that there exists numeric matrix A such that equation F = 0
    is equivalent to linear equation Ax = b, this function returns 
    tuple (A, b)
    """
    A = []
    b = []
    for row in F:
        coeffs = expand(row).as_coefficients_dict()
        A.append([float(coeffs[x]) for x in variables])
        b.append(-float(coeffs[1]))
    return np.array(A), np.array(b)

A, b = symbolic_to_matrix(M * X - B, [x1, x2, f1])
# A
# array([[ 2.,  3., -1.],
#       [ 5.,  6.,  0.],
#       [ 8.,  9.,  0.]])
# b
# array([ -0.,  50.,  60.])

np.linalg.solve(A, b)
# array([-30.        ,  33.33333333,  40.        ])
# the same answer 

вы можете представить систему линейных уравнений в виде списка списков, а затем определить очень простую функцию для ее решения, как показано ниже:

def solve(equations):
     """
     the constants of a system of linear equations are stored in a list for each equation in the system
     for example the system below:
          2x+9y-3z+7w+8=0
          7x-2y+6z-1w-10=0
          -8x-3y+2z+5w+4=0
          0x+2y+z+w+0=0
     is expressed as the list:
          [[2,9,-3,7,8],[7,-2,6,-1,-10],[-8,-3,2,5,4],[0,2,1,1,0]]
     """
     for i in equations:
          if len(i)<>(len(equations)+1):
               raise ValueError("your equation system has not a valid format")
     lists=[] # I failed to name it meaningfully
     for eq in range(len(equations)):
          #print "equations 1", equations
          #find an equation whose first element is not zero and call it index
          index=-1
          for i in range(len(equations)):
               if equations[i][0]<>0:
                    index=i;
                    break;
          if index==-1:
               raise ValueError("your equation system can not be solved")
          #print "index "+str(eq)+": ",index
          #for the equation[index] calc the lists next item  as follows
          lists.append([-1.0*i/equations[index][0] for i in equations[index][1:]])
          #print "list"+str(eq)+": ", lists[-1]
          #remove equation[index] and modify the others
          equations.pop(index)
          for i in equations:
               for j in range(len(lists[-1])):
                    i[j+1]+=i[0]*lists[-1][j]
               i.pop(0)



 lists.reverse()
 answers=[lists[0][0]]
 for i in range(1,len(lists)):
      tmpans=lists[i][-1]
      for j in range(len(lists[i])-1):
           tmpans+=lists[i][j]*answers[-1-j]
      answers.append(tmpans)
 answers.reverse()
 return answers