Leitura Complementar: Como resolver sistema de equações em python

Muitos problemas de equilíbrio químico envolvem a resolução de equações de sistemas lineares e não lineares. Os exemplos servirão como base para resolução de exercícios de equilíbrio químico, como por exemplo no cálculo da composição de equilibrio de reações químicas, resolvendo um sistema de equações não lineares. Os exemplo resolvidos motram várias maneiras de se resolver sistemas de equações em python.

Para resolução dos exemplos serão utilizados as bibliotecas Numpy e SciPy.

# Módulos necessários para resolver sistema de equações
import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import least_squares

O caso mais simples são sistemas lineares. Nestes casos, a resolução do sistema pode ser feita pela regra de Cramer: \( A = A^{-1} B \). Basta escrever as equações no formato de um sistema de equações e resolvê-lo usando a rotina np.linalg.solve(A, B). Onde A e B são duas matrizes e a instrução linalg.solve() calculará a solução do sistema de equações lineares.

Ex.) Resolva o sistema de equações abaixo:

\[ x + 2y = 1 \]

\[ 3x + 5y = 2 \]

Solução

Em primeiro lugar é preciso definir o sistema de equaçõesa ser resolvido. Isto é feito escrevendo as funções forma que sejam calculadas na forma de uma equação matricial:

\( A = A^{-1} B \), onde:

\[\begin{split} A = \begin{bmatrix}1 & 2\\ 3 & 5\\ \end{bmatrix} \end{split}\]

\[\begin{split} B = \begin{bmatrix}1 \\ 2 \\ \end{bmatrix} \end{split}\]

#Definição das matrizes

A = np.array([[1, 2], [3, 5]])

B = np.array([1, 2])

x = np.linalg.solve(A, B)

print('x = ' + format(x[0] , ' 6.3f'), 'y = ' + format(x[1] , ' 6.3f'))

x = -1.000 y =  1.000

Ex.) Resolva o sistema de equações abaixo:

\[ 4x^2 + y^2 = 1 \]

\[ x^2 + 9y^2 = 1 \]

Solução

Novamente a primeira coisa a ser feita é definir o sistema de equações a ser resolvido. Isto é feito escrevendo as funções forma que sejam calculadas suas raízes no formato do python. Neste exemplo deve-se escrever da forma:

\[ 4x^2 + y^2 - 1 = 0 \]

\[ x^2 + 9y^2 -1 = 0 \]

# Definição das funções

def equations(vars):
    x, y = vars
    eq1 = 4*x**2 + y**2 - 1
    eq2 = x**2 + 9*y**2 - 1
    return [eq1, eq2]

Existem duas maneiras de se resolver este sistema. Uma faz uso da rotina least_squares, que faz uso dos mínimos quadrados para encontrar a solução impondo restrições as variáveis. A outra forma é usar a rotina fsolve que busca as raízes das funções definidas. Em ambos os casos é necessário informar uma estimativa inicial para os valores que devem ser calculados. A sintaxe das duas instruções é muito similar.

#Resolvendo com fsolv. Os números entre parênteses são o chute inicial para x e y.

x, y =  fsolve(equations, (3, 2))
print('Solução com fsolve= ' +  format(x , ' 6.5f') + format(y , ' 6.5f'))

print()

#Resolvendo com lest_squares. Os números entre parênteses são o chute inicial para x e y.

res = least_squares(equations, (3, 2))

print('Solução com least_squares= ' + format(res.x[0] , ' 6.5f') + format(res.x[1] , ' 6.5f'))

  

Solução com fsolve=  0.47809 0.29277

Solução com least_squares=  0.47809 0.29277

Como pode ser observado os resultados são iguais. A escolha do método de resolução dependerá do problema em questão e se existe a necessidade de impor restrições a variáveis do problema.