2025, Sep 26 23:17

Матрица вторых производных в SymPy: как отфильтровать нули и сохранить символику

Разбираем, почему Derivative(x, y).is_constant() в SymPy даёт True, и как построить «второй якобиан» с обнулением нулей, сохранив символические производные.

Собрать матрицу частных производных второго порядка в SymPy кажется простой задачей — пока внезапно весь результат не оказывается «константным». Если вы пытались отсечь тривиальные нули с помощью is_constant() и увидели сплошные True, вы столкнулись с тем, как SymPy трактует производные независимых символов. Разберёмся, что происходит, и как получить практичный «якобиан второго порядка» с корректным обнулением явных нулей.

Постановка задачи

Цель — сформировать матрицу частных производных первого порядка для набора символов, а затем снова продифференцировать её по каждому символу, получив структуру второго порядка. Ниже — минимальный пример с ASCII-именами и той же логикой.

from sympy import symbols, Matrix, zeros, Derivative
sym_names = ["Delta", "Gamma", "v"]
def make_vars():
    return symbols(" ".join(sym_names))
def first_order_partials():
    grid = []
    vars_seq = make_vars()
    for i in range(len(sym_names)):
        left = vars_seq[i]
        row = []
        for j in range(len(sym_names)):
            right = vars_seq[j]
            if left == right:
                row.append(1)
            else:
                row.append(Derivative(left, right))
        grid.append(row)
    return vars_seq, Matrix(grid)
def second_order_partials():
    vars_seq, jmat = first_order_partials()
    out = zeros(len(vars_seq)**2, len(vars_seq))
    nrows, ncols = jmat.shape
    block_idx = 0
    for k in range(len(vars_seq)):
        denom = vars_seq[k]
        for r in range(nrows):
            for c in range(ncols):
                num = jmat[r, c]
                d2 = Derivative(num, denom)
                out[block_idx * len(vars_seq) + r, c] = d2
        block_idx += 1
    return out

После этого при попытке отфильтровать элементы по is_constant() всё выглядит как константа.

Почему всё кажется константным

Ключевой момент: частная производная одного независимого символа по другому равна нулю. Именно это и отражает SymPy:

from sympy import symbols, Derivative, diff
x, y = symbols('x y')
d = Derivative(x, y)
# Показывает невычисленную производную
print(d)
# SymPy знает, что это константа (равна 0)
print(d.is_constant())
# Вычисление производной даёт 0
print(d.doit())
print(diff(x, y))

Поэтому Derivative(x, y).is_constant() возвращает True: производная равна 0. Та же логика для Derivative(1, x), она тоже сводится к 0. Если в матрице первого порядка стоят константы или производные независимых символов, любые вторые производные по ним останутся константами — отсюда сплошные True.

Как получить неконстантные производные

Если вы ожидаете ненулевые производные, выражение должно действительно зависеть от переменной. В SymPy для этого используют функцию, а не одиночный символ:

from sympy import symbols, Function, Derivative
x = symbols('x')
f = Function('f')
u = Derivative(f(x), x)
print(u)          # остаётся невычисленной
print(u.doit())   # по-прежнему символическая производная f(x)
print(u.is_constant())  # False

Это единственный способ задать выражение, которое действительно зависит от x, а не отдельный независимый символ. Нельзя заставить набор символов быть функциями друг друга.

Корректное обнуление тривиальных нулей

Если вы хотите оставлять производные невычисленными, но убирать гарантированные нули, не проверяйте сам объект Derivative через is_constant(). Вместо этого выясняйте, является ли константой выражение, которое дифференцируют, и действуйте исходя из этого. Такой подход удаляет, например, Derivative(1, x), но оставляет Derivative(x, y) в невычисленном виде, если нужно сохранять символику. Когда требуется свернуть до нуля, просто ставьте 0 там, где d2.expr является константой.

from sympy import symbols, Matrix, zeros, Derivative
sym_names = ["Delta", "Gamma", "v"]
def make_vars():
    return symbols(" ".join(sym_names))
def first_order_partials():
    grid = []
    vars_seq = make_vars()
    for i in range(len(sym_names)):
        left = vars_seq[i]
        row = []
        for j in range(len(sym_names)):
            right = vars_seq[j]
            if left == right:
                row.append(1)
            else:
                row.append(Derivative(left, right))
        grid.append(row)
    return vars_seq, Matrix(grid)
def second_order_partials_collapsed():
    vars_seq, jmat = first_order_partials()
    out = zeros(len(vars_seq)**2, len(vars_seq))
    nrows, ncols = jmat.shape
    block_idx = 0
    for k in range(len(vars_seq)):
        denom = vars_seq[k]
        for r in range(nrows):
            for c in range(ncols):
                d2 = Derivative(jmat[r, c], denom)
                # Убираем тривиальные нули (выражение — константа)
                out_idx = block_idx * len(vars_seq) + r
                out[out_idx, c] = 0 if d2.expr.is_constant() else d2
        block_idx += 1
    return out

С таким подходом всё вроде Derivative(1, Delta) сразу превращается в 0. Элементы, у которых дифференцируемое выражение не является константой, сохраняются как символические производные.

Почему это важно

Понимание того, что именно SymPy считает «константой», избавляет от ложных выводов. Вызов is_constant() для Derivative говорит о том, зависит ли значение самой производной от переменных, а не о том, переменно ли внутреннее выражение. Поэтому Derivative(x, y).is_constant() возвращает True: производная равна нулю. Если нужно решить, будет ли производная тривиальным нулём ещё до вычисления, проверяйте выражение, которое дифференцируют, через di.expr.is_constant().

Выводы

При построении матриц производных более высокого порядка в SymPy важно различать независимые символы и функции от символов. Ожидаете изменяемость — задавайте её явно через Function. Если нужно сохранить символическую структуру и при этом убрать гарантированные нули, проверяйте .expr у производной, а не сам объект Derivative, и подставляйте 0 в соответствующие позиции. Этот небольшой сдвиг согласует код с семантикой SymPy и даёт более чистую и информативную матрицу второго порядка.

Статья основана на вопросе на StackOverflow от James Strieter и ответе Oscar Benjamin.