SymPy — проблемы при использовании многих параметров в dsolve

Я использую SymPy версии 0.7.3 и столкнулся с некоторыми проблемами при использовании функции dsolve. Кажется, что у dsolve возникают трудности, когда входное уравнение имеет слишком много параметров.

Я попытался решить следующее уравнение:

from sympy import *
p = Function('p')
t, u1, u2, u3, u4, u5 = symbols('t u1 u2 u3 u4 u5')
eq = Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2)))
eq
     Out: Derivative(p(t), t) == (-u3 - u4)*p(t) + exp(-t*(u1 + u2))
%time dsolve(eq)

И получил:

CPU times: user 213.11 s, sys: 0.00 s, total: 213.11 s
Wall time: 213.12 s
       p(t) == (C1 + Piecewise((t*u1/(u1 + u2) + t*u2/(u1 + u2), u3 == u1 + u2 - u4), (-exp(t*u3)*exp(t*u4)/(u1*exp(t*u1)*exp(t*u2) + u2*exp(t*u1)*exp(t*u2) - u3*exp(t*u1)*exp(t*u2) - u4*exp(t*u1)*ex
p(t*u2)), True)))*exp(-t*(u3 + u4))

(Заняло 213,12 секунды!)

Затем я заменил u1+u2 на u5:

eq = Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))).subs(u1+u2, u5)
eq
     Out:Derivative(p(t), t) == (-u3 - u4)*p(t) + exp(-t*u5)
%time dsolve(eq)

и получил:

CPU times: user 1.62 s, sys: 0.00 s, total: 1.62 s
Wall time: 1.62 s
        p(t) == (C1 + Piecewise((t, u3 == -u4 + u5), (exp(t*u3)*exp(t*u4)/(u3*exp(t*u5) + u4*exp(t*u5) - u5*exp(t*u5)), True)))*exp(-t*(u3 + u4))

(Всего 1,62 секунды!)

Я пробовал использовать разные подсказки, но это не помогло...

Я также заметил, что в гораздо более сложных функциях dsolve падает, но при замене некоторых постоянных параметров работает быстро.

Знаете ли вы, в чем причина этого явления? Есть ли способ решить эту проблему автоматически?


person user5497    schedule 26.02.2014    source источник


Ответы (2)


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

В общем решении упрощающему может быть очень трудно восстановить, что определенные символы встречаются только в заданной комбинации и, таким образом, могут рассматриваться как один символ. Так что упрощенцу без необходимости приходится иметь дело с большим пространством поиска. Кроме того, он должен убедиться, что правильно обрабатывает все граничные случаи (u1 может быть ‹ 0, > 0, = 0, сложный, ...?)

Я нашел очень выгодным сделать как можно больше упрощений вручную, прежде чем отдать задачу символическому решателю. Группировка переменных вручную (как в вашем примере) является одним из полезных методов. Другой метод заключается в разработке нормализации или установке одного параметра в 1. Например, при работе с полиномом a x^2 + b x + c большую часть времени проблема x^2 + B x + C эквивалентна нам. (Поскольку мы на самом деле уверены, что a != 0, но забыли сказать решателю, верно?...) Но для решателя имеет большое значение, если количество символьных переменных уменьшается на 1.

В какой-то момент эти символьные решатели, безусловно, станут достаточно умными, чтобы сгруппировать переменные вместе, прежде чем пытаться решить проблему. Однако в настоящее время кажется, что вмешательство человека по-прежнему требуется.

С другой стороны, трудно представить, что символьный решатель автоматически распознает более сложные преобразования, упрощающие задачу, такие как переход от декартовых координат к полярным или изменение переменных, таких как l=x+iy, r=x-iy, которые вовсе не нужны. очевидно, но известно, что оно полезно для решения определенных задач.

Обновить

Кажется, лучшее, что мы можем сделать для этой проблемы, это установить a = u3 + u4 и a + b = u1 + u2. Помимо уменьшения количества параметров с 4 до 2, теперь появляется особый случай, когда b == 0, поэтому мы можем легко указать решателю игнорировать его:

from sympy import *
p = Function('p')
t, a = symbols('t a')
b = symbols('b', nonzero=True)
eq = Eq(Derivative(p(t),t), -a*p(t) + exp(-(a + b)*t))
dsolve(eq)
# -> p(t) == (C1 - exp(-b*t)/b)*exp(-a*t)
# (0.75 s)

Таким образом, помогая решателю, избегая особого случая, время решения снова сократилось вдвое (время в моей системе было похоже на ваше). В случае, если частный случай b == 0 действительно имеет значение, его можно легко восстановить из разложения в ряд Тейлора exp(-b*t) ~ 1 - b*t.

В общем, указание того, что переменные являются действительными, ненулевыми, строго больше нуля и т. д., также очень полезно, чтобы избежать зависания решателя в особых случаях. Иногда действительно лучше решать задачу отдельно для x < 0, x > 0 и x == 0 (избегая пресловутых sqrt(x^2) выражений, которые решатель не может упростить дальше без знания знака x).

person Stefan    schedule 26.02.2014

Проблема становится немного яснее, если вы используете подсказку "1st_linear_Integral", которая возвращает то, что было бы сделано как невычисленные интегралы (я использую 1st_linear, потому что это первый метод, возвращаемый classify_ode(eq), то есть тот, который используется dsolve по умолчанию):

In [61]: dsolve(Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))), hint='1st_linear_Integral')
Out[61]:
       ⎛     ⌠                                 ⎞
       ⎜     ⎮                ⌠                ⎟   ⌠
       ⎜     ⎮                ⎮ (u₃ + u₄) dt   ⎟  -⎮ (u₃ + u₄) dt
       ⎜     ⎮  -t⋅u₁  -t⋅u₂  ⌡                ⎟   ⌡
p(t) = ⎜C₁ + ⎮ ℯ     ⋅ℯ     ⋅ℯ               dt⎟⋅ℯ
       ⎝     ⌡                                 ⎠

In [62]: dsolve(Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))).subs(u1+u2, u5), hint='1st_linear_Integral')
Out[62]:
       ⎛     ⌠                          ⎞
       ⎜     ⎮         ⌠                ⎟   ⌠
       ⎜     ⎮         ⎮ (u₃ + u₄) dt   ⎟  -⎮ (u₃ + u₄) dt
       ⎜     ⎮  -t⋅u₅  ⌡                ⎟   ⌡
p(t) = ⎜C₁ + ⎮ ℯ     ⋅ℯ               dt⎟⋅ℯ
       ⎝     ⌡                          ⎠

У алгоритма интеграции есть проблема с exp(a*x)*exp(b*x), тогда как у него нет проблем с exp((a + b)*x). В основном существует более быстрый алгоритм с более ограниченной областью действия, который вызывается первым, а затем более медленный алгоритм с большей областью действия, который вызывается в случае сбоя быстрого алгоритма. Быстрый может обрабатывать exp((a + b)*x), но пока не exp(a*x)*exp(b*x).

Простой обходной путь в этом конкретном случае - объединить экспоненциалы, используя powsimp:

In [67]: a = dsolve(Eq(Derivative(p(t),t), -(u3 + u4)*p(t) + exp(-t*(u1 + u2))), hint='1st_linear_Integral')

In [68]: powsimp(a, deep=True)
Out[68]:
       ⎛     ⌠                                  ⎞
       ⎜     ⎮                 ⌠                ⎟   ⌠
       ⎜     ⎮  -t⋅u₁ - t⋅u₂ + ⎮ (u₃ + u₄) dt   ⎟  -⎮ (u₃ + u₄) dt
       ⎜     ⎮                 ⌡                ⎟   ⌡
p(t) = ⎜C₁ + ⎮ ℯ                              dt⎟⋅ℯ
       ⎝     ⌡                                  ⎠

In [69]: %time powsimp(a, deep=True).doit()
CPU times: user 261 ms, sys: 2.11 ms, total: 263 ms
Wall time: 262 ms
Out[69]:
       ⎛     ⎛⎧              t                for u₁ + u₂ - u₃ - u₄ = 0⎞⎞
       ⎜     ⎜⎪                                                        ⎟⎟
       ⎜     ⎜⎪  -t⋅u₁ - t⋅u₂ + t⋅(u₃ + u₄)                            ⎟⎟  -t⋅(u₃ + u₄)
p(t) = ⎜C₁ + ⎜⎨-ℯ                                                      ⎟⎟⋅ℯ
       ⎜     ⎜⎪─────────────────────────────          otherwise        ⎟⎟
       ⎜     ⎜⎪      u₁ + u₂ - u₃ - u₄                                 ⎟⎟
       ⎝     ⎝⎩                                                        ⎠⎠

В общем, предложения Стефана могут быть верными, а могут и не быть. Теоретически, CAS не должно волновать, насколько сложны ваши символьные константы, потому что это просто константы. На самом деле есть проблемы, потому что он объединяет константы, а затем ему нужно посмотреть, не отменяются ли вещи, и так далее. Кроме того, небольшие различия, подобные приведенным выше, могут привести к большим различиям в реальных путях алгоритмов. Математически два выражения могут быть одинаковыми, но алгоритмы зависят от того, как они выглядят структурно. Как правило, более простые выражения работают лучше.

Если вам часто приходится заменять подвыражения символами, вы можете использовать cse, чтобы помочь.

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

Кстати, SymPy 0.7.3 немного устарел. Только что вышла самая новая версия 0.7.5.

person asmeurer    schedule 01.03.2014