How to prepare a system of equations for python

Question

This VIDEO is a really nice demonstration for how to do it in a simplified case - in case of three equations and three variables.

Let's say I wan to solve the following system

for variables f1, x1 and x2. Since this is a rather small system I could easily do it manually. But this is a simplified example - in reality my system consists of 100 variables and 100 equations.

So my question is how to separate variables in order to solve this system? How to gather all the variables in one vector and rewrite the system so I can solve it? In the end all I want is numerical values for f1, x1 and x2.

ps.: I just made this system by inserting random numbers. I'm not sure the system cae be solved but... you get the idea. (adjust the numerical values in that case).

Answer 1

you can represent a system of linear equations in a list of lists and then define a very simple function to solve it as demonstrated below:

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

Answer 2

As far as I understand, you have to adjust the matrix of your system to take care of RyA and other variables that currently in the right hand side. You can do it manually (in which case this question is out of the scope of this site, it is purely mathematical excercise) or use e.g. sympy instead of np.linalg.solve() which can do the algebra part of the problem for you:

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}

Note that sympy can be slower in solving numeric linear systems than numpy.linalg, so you may want to use sympy to do algebraic part of the work, calculate the matrix and the right hand side, and then use numpy.linalg.solve to solve it.

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