SymPy - проблемы при использовании многих параметров в dsolve
Я использую SymPy версии 0.7.3 и столкнулся с некоторыми проблемами при использовании функции 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, но при замене некоторых постоянных параметров он работает быстро.
Вы знаете, в чем причина этого явления? Есть ли способ решить эту проблему автоматически?
2 ответа
Проблема немного яснее, если вы используете "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, была только что выпущена.
Я отметил ту же проблему с 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).