Способы решения систем уравнений в Python с примерами кода

Пример 1. Решение системы нелинейных уравнений с помощью scipy.optimize.fsolve

Система нелинейных уравнений, например, из двух уравнений с двумя неизвестными, решается численно методом Ньютона. Функция принимает вектор неизвестных и возвращает вектор невязок.

from scipy.optimize import fsolve
import numpy as np

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

initial_guess = [1, 1]
solution = fsolve(equations, initial_guess)
print('Решение:', solution)
print('Проверка:', equations(solution))

Решение: [ 0.62007947 -1.27658617]
Проверка: [4.969432223000867e-10, 4.494532170476073e-10]

Необходимо задать начальное приближение; разные начальные точки могут приводить к разным корням или расходимости. В случае вырожденного якобиана может возникать ошибка.

Пример 2. Решение системы дифференциальных уравнений (через дискретизацию) с использованием разреженных матриц

Пусть требуется решить краевую задачу для уравнения Пуассона в одномерном пространстве. После дискретизации конечными разностями получается трехдиагональная система, которую эффективно решать через scipy.sparse.linalg.spsolve.

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

n = 10
h = 1.0 / (n + 1)
x = np.linspace(0, 1, n+2)[1:-1]  # внутренние точки
# Построение матрицы второй производной
diagonals = [np.ones(n), -2*np.ones(n), np.ones(n)]
A = diags(diagonals, [-1,0,1], shape=(n,n)).toarray() * (1/h**2)
# Правая часть (источник)
f = np.sin(np.pi * x) * np.pi**2
# Граничные условия: u(0)=0, u(1)=0 уже учтены
b = f - 0  # вклад границ отсутствует
u_num = spsolve(csr_matrix(A), b)
print('Численное решение в точках:', u_num)
# Сравнение с аналитическим u_analitic = sin(pi*x)
u_exact = np.sin(np.pi * x)
print('Максимальная ошибка:', np.max(np.abs(u_num - u_exact)))

Численное решение в точках: [0.30901699 0.58778525 0.80901699 0.95105652 1.         0.95105652 0.80901699 0.58778525 0.30901699 0.        ]
Максимальная ошибка: 0.0013... (зависит от n)

Пример 3. Решение системы с комплексными коэффициентами

NumPy и SciPy поддерживают комплексные числа. Если матрица и вектор содержат мнимые единицы, решение будет комплексным.

import numpy as np

A = np.array([[1+2j, 3-1j],
              [0+1j, 2+0j]])
b = np.array([4+0j, 1+1j])
x = np.linalg.solve(A, b)
print('Решение:', x)

Решение: [ 0.48-0.84j -0.26+0.52j]

Пример 4. Решение системы уравнений с параметрами и последующая подстановка значений

Используем SymPy для символьного решения с параметрами, затем подставляем численные значения.

import sympy as sp

a, b, c, x, y = sp.symbols('a b c x y')
eq1 = sp.Eq(a*x + b*y, c)
eq2 = sp.Eq(b*x - a*y, 0)
sol = sp.solve([eq1, eq2], (x, y), dict=True)
print('Символьное решение:', sol)
# Подстановка a=2, b=3, c=1
subs_values = {a:2, b:3, c:1}
num_sol = [{k: v.subs(subs_values).evalf() for k,v in s.items()} for s in sol]
print('Численное решение:', num_sol)

Символьное решение: [{x: a*c/(a**2 + b**2), y: b*c/(a**2 + b**2)}]
Численное решение: [{x: 0.153846153846154, y: 0.230769230769231}]

Пример 5. Использование библиотеки TensorFlow для решения системы (градиентными методами)

Хотя это не классический подход, можно решить систему минимизацией невязки с помощью оптимизаторов TensorFlow. Полезно при необходимости интегрировать решение в нейросетевой пайплайн.

import tensorflow as tf

A = tf.constant([[3.,2.,-1.],[2.,-2.,4.],[-1.,0.5,-1.]])
b = tf.constant([1.,-2.,0.])
x = tf.Variable(tf.zeros(3))

loss_fn = lambda: tf.reduce_sum(tf.square(tf.linalg.matvec(A, x) - b))
optimizer = tf.optimizers.Adam(learning_rate=0.1)
for step in range(1000):
    optimizer.minimize(loss_fn, [x])
print('Приближенное решение:', x.numpy())

Приближенное решение: [ 1.0000001 -2.0000002 -2.0000005]