Решение уравнений в Python (SciPy)

Решение уравнений в Python выполняется средствами библиотеки SciPy: модуль scipy.optimize находит корни нелинейных уравнений и систем и подбирает параметры моделей, а scipy.integrate вычисляет определённые интегралы и решает дифференциальные уравнения. Для линейных систем используется NumPy. Это стандартный набор инструментов инженерных расчётов, заменяющий ручные итерации и электронные таблицы.

Ниже разобраны основные задачи с рабочими примерами кода: поиск корня нелинейного уравнения через fsolve, решение систем уравнений, численное интегрирование и подбор параметров по экспериментальным данным. Все функции приведены по действующей документации SciPy.

Инструменты SciPy и NumPy

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

Численное решение даёт приближённый ответ с контролируемой точностью там, где аналитическое решение недоступно или неудобно.

Корни нелинейных уравнений: fsolve

Функция fsolve находит корни уравнения вида f(x) = 0 по начальному приближению. Уравнение записывают как функцию, возвращающую значение левой части, а fsolve подбирает x, при котором оно обращается в ноль.

import numpy as np
from scipy.optimize import fsolve

def f(x):
    return x - np.cos(x)      # уравнение x = cos(x)

root = fsolve(f, x0=0.5)      # x0 — начальное приближение
print(root)                  # [0.73908513]

Если известен интервал, на концах которого функция меняет знак, надёжнее применять brentq — метод Брента гарантированно сходится на таком интервале. Для гладкой функции с известной производной подойдёт newton (метод Ньютона или секущих).

from scipy.optimize import brentq

# f(0) = -1 < 0, f(1) = 1 - cos(1) > 0 — знак меняется
root = brentq(f, 0, 1)
print(root)                  # 0.7390851332151607

Инженерный пример: коэффициент трения

Уравнение Колбрука для коэффициента гидравлического трения задано неявно и не решается аналитически — классический случай для fsolve.

import numpy as np
from scipy.optimize import fsolve

Re = 1e5          # число Рейнольдса
eps_D = 0.0002    # относительная шероховатость

def colebrook(lam):
    return 1/np.sqrt(lam) + 2*np.log10(eps_D/3.7 + 2.51/(Re*np.sqrt(lam)))

lam = fsolve(colebrook, x0=0.02)
print(lam)        # коэффициент трения, около 0.019

Результат fsolve зависит от начального приближения: при нескольких корнях функция вернёт тот, что ближе к x0. Начальное приближение стоит выбирать из физического смысла задачи.

Системы уравнений

Нелинейные системы

Для системы нелинейных уравнений функция принимает вектор неизвестных и возвращает список значений левых частей. Начальное приближение задаётся вектором.

import numpy as np
from scipy.optimize import fsolve

def system(v):
    x, y = v
    return [x**2 + y**2 - 25,   # x^2 + y^2 = 25
            x**2 - y]           # x^2 = y

sol = fsolve(system, x0=[1, 1])
print(sol)                     # около [2.13, 4.52]

Более гибкая альтернатива — root: она позволяет выбрать метод решения (по умолчанию гибридный алгоритм Пауэлла hybr, доступны также lm, broyden1 и другие) и возвращает объект с полем x и признаком сходимости success.

Линейные системы

Систему линейных уравнений вида A·x = b решают напрямую через NumPy, без итераций и начального приближения.

import numpy as np

A = np.array([[3, 2, -1],
              [2, -2, 4],
              [-1, 0.5, -1]])
b = np.array([1, -2, 0])

x = np.linalg.solve(A, b)
print(x)                       # [ 1. -2. -2.]

Функция numpy.linalg.solve предназначена именно для линейных систем с квадратной невырожденной матрицей. Применять для них fsolve не нужно — прямое решение точнее и быстрее.

Численное интегрирование

Определённый интеграл: quad

Функция quad вычисляет определённый интеграл функции на заданном отрезке методом адаптивной квадратуры и возвращает кортеж из двух значений: результат и оценку погрешности. Пределы могут быть бесконечными (np.inf).

import numpy as np
from scipy.integrate import quad

def f(x):
    return np.exp(-x**2)

value, err = quad(f, 0, np.inf)
print(value)                   # 0.8862269254527579
print(err)                     # оценка погрешности

Дифференциальные уравнения: solve_ivp

Обыкновенные дифференциальные уравнения интегрирует solve_ivp. Правую часть записывают как функцию f(t, y), задают интервал интегрирования и начальные условия. По умолчанию используется явный метод Рунге — Кутты. Результат содержит массивы времени t и решения y.

import numpy as np
from scipy.integrate import solve_ivp

k, T_env = 0.1, 20.0           # закон охлаждения Ньютона

def cooling(t, T):
    return -k * (T - T_env)

sol = solve_ivp(cooling, t_span=[0, 60], y0=[90],
                t_eval=np.linspace(0, 60, 7))
print(sol.t)                   # моменты времени
print(sol.y[0])                # температура в эти моменты

Уравнение более высокого порядка предварительно сводят к системе уравнений первого порядка, вводя производные как дополнительные компоненты вектора состояния.

Подбор параметров

Функция curve_fit подбирает параметры заданной модели по экспериментальным данным методом наименьших квадратов. Модель записывают как функцию, где первый аргумент — независимая переменная, а остальные — искомые параметры.

import numpy as np
from scipy.optimize import curve_fit

def model(x, a, b, c):
    return a * np.exp(-b * x) + c

x_data = np.array([0, 1, 2, 3, 4, 5])
y_data = np.array([9.9, 6.1, 3.9, 2.6, 2.1, 1.8])

popt, pcov = curve_fit(model, x_data, y_data, p0=[8, 0.5, 1])
print(popt)                    # найденные параметры a, b, c

Функция возвращает два массива: popt — оптимальные параметры и pcov — ковариационную матрицу. Диагональ pcov даёт дисперсии оценок параметров; квадратный корень из них — стандартные отклонения. Аргумент p0 задаёт начальные приближения параметров: для устойчивой сходимости их полезно указывать исходя из ожидаемого вида зависимости.

Практические рекомендации

1. Задавайте разумное начальное приближение. Для fsolve, newton и curve_fit оно определяет, к какому корню или набору параметров сойдётся алгоритм.
2. Используйте brentq при известном интервале. Если на концах отрезка функция меняет знак, метод Брента надёжнее методов, стартующих из одной точки.
3. Проверяйте сходимость. У root и solve_ivp есть поля success и status; подставьте найденное решение обратно в уравнение и убедитесь, что невязка близка к нулю.
4. Векторизуйте функции. Операции NumPy (np.exp, np.sqrt) работают с массивами и ускоряют вычисления по сравнению с поэлементными циклами.
5. Управляйте точностью. Параметры допусков (xtol, rtol, atol) регулируют точность и число итераций под конкретную задачу.

Частые вопросы