Python Sympy, разложение частичных дробей: результат не такой, как ожидалось
Я использую библиотеку Python Sympy для расчета разложения частичных дробей следующей рациональной функции:
Z = (6.43369157032015e-9*s^3 + 1.35203404799555e-5*s^2 + 0.00357538393743079*s + 0.085)/(4.74334912634438e-11*s^4 + 4.09576274286244e-6*s^3 + 0.00334241812250921*s^2 + 0.15406018058983*s + 1.0)
Я вызываю метод apart() так:
Z_f = apart(Z, full=True).doit()
И вот результат:
-1.30434068814287e+23/(s + 85524.0054884464) + 19145168.0/(s + 774.88576677949) + 91.9375/(s + 40.7977016133126) + 2.59375/(s + 7.79746609204661)
Как видно из результата, это остатки:
- -3.60418953263334e+22
- -4789228.25000000
- 11.0468750000000
- 0,787109375000000
Примечание. Я идентифицирую их следующим образом:
Z_f_terms = Z_f.as_ordered_terms()
а затем извлеките их в цикле for.
Проблема в том, что, используя другие инструменты, я получаю разные остатки.
Используя GNU Octave residue(), я получаю следующие остатки:
- 133,6
- 1,0776
- 0,39501
- 0,56426
Через Wolphram Alpha я получаю те же остатки, что и Octave см. здесь результаты (пожалуйста, подождите, пока не увидите разложение неполных дробей).
Почему библиотека Sympy не работает должным образом?
Полюсы всегда правильные. Только остатки не такие, как я ожидаю.
Спасибо за любое предложение.
почему бы не использовать scipy.signal.residue?
@lastchance Я использую «.as_ordered_terms». Пожалуйста, ознакомьтесь с обновленным вопросом.
@Davide_sd Спасибо, я только что попробовал, и это работает. Единственный случай, который не работает, — это Sympy. Может ли кто-нибудь воспроизвести это?
Да, @JackSlater, я могу воспроизвести твою симпатичную проблему. Я также попробовал использовать гораздо более простое выражение: Z = (12*s**2 - 44*s + 36)/(s**3-6*s**2+11*s-6). Посмотрите, что произойдет, если вы измените это число 12 (int) на 12,0 (float). Я думаю, что символическая арифметика Sympy сталкивается с проблемами точности чисел с плавающей запятой, когда у вас есть члены очень разных порядков величины (в вашем случае коэффициенты от 1e-11 до 1,0). Я бы использовал предложение @Davide_sd, а не символьную арифметику. Или вы можете сделать, как я, и найти корни с помощью полиномиального класса, а затем использовать правило перекрытия, чтобы найти остатки.
Покажите, пожалуйста, как вы получили приведенные вами цифры. Детали отсутствуют в «извлеченных в цикле for», и не очевидно, как полученные вами числа связаны с терминами, из которых цикл for их извлек.
@lastchance Спасибо, я попробовал и предложение Davide_sd, и предложение Оскара Бенджамина, и расчет Scipy кажется быстрее.
@JackSlater, кажется, нет. Это! Используйте SymPy только в том случае, если вам действительно нужно символьное решение, в противном случае лучше использовать числовые библиотеки.
Ответы 1
Проблема связана с поплавками. Предположительно, что-то в алгоритме, используемом apart(full=True) или RootSum.doit(), численно нестабильно для приближенных коэффициентов.
Если вы преобразуете коэффициенты в рациональные числа перед вычислением RootSum, то evalf даст ожидаемый результат:
In [27]: apart(nsimplify(Z), full=True).evalf()
Out[27]:
133.599202650992 1.07757928431867 0.395006955518971 0.564264854137341
──────────────────── + ─────────────────── + ──────────────────── + ────────────────────
s + 85524.0054884464 s + 774.88576677949 s + 40.7977016133126 s + 7.79746609204661
In [29]: apart(Z, domain = "QQ", full=True).evalf() # rational field
Out[29]:
133.599202650994 1.07757928431867 0.395006955518971 0.564264854137341
──────────────────── + ─────────────────── + ──────────────────── + ────────────────────
s + 85524.0054884464 s + 774.88576677949 s + 40.7977016133126 s + 7.79746609204661
Обратите внимание, что общий алгоритм и подход, используемые для full=True, на самом деле предназначены для нахождения точного представления разложения, избегая радикалов. Он не предназначен для поплавков. Вероятно, хотя обычный apart без full=True должен вычислить это разложение.
Спасибо за подробное объяснение, это действительно решает проблему. Я также попробовал «отдельно» без «full=True», но это не сработало. У меня нет опыта работы с этой библиотекой. Следует ли сообщить об этой проблеме разработчикам? Кстати, вам следует отредактировать свой ответ следующим образом: domain='QQ'.
Вам не нужно использовать domain='QQ', если вы это сделаете from sympy import QQ. Мой ответ продемонстрирован в сеансе isympy, в котором неявно используется from sympy import *.
Исходная ошибка будет исправлена в SymPy 1.13: github.com/sympy/sympy/pull/26649
Большое спасибо за то, что открыли проблему, исправили ее и очень быстро сообщили об этом. Очень приятно, поздравляю. Супер!
Не могли бы вы объяснить, как вы определяете остатки, которые вы указали, из приведенного вами разложения неполных дробей. Кстати, если ваши корни — это простые полюса, то вы можете определить верхнюю часть дробей по «правилу сокрытия».