Обновлено: Исходный пост слишком расплывчатый. Я ищу алгоритм для решения большой системы, решаемой, линейной IVP, которая может обрабатывать очень маленькие значения с плавающей запятой. Решение для собственных векторов и собственных значений невозможно с помощью numpy.linalg.eig()
, поскольку возвращаемые значения являются сложными и не должны быть, он также не поддерживает numpy.float128
, а матрица не симметрична, поэтому numpy.linalg.eigh()
не будет работать. Sympy мог сделать это, учитывая бесконечное количество времени, но после 5 часов работы я сдался. scipy.integrate.solve_ivp()
работает с неявными методами (пробовали Радау и BDF), но вывод совершенно неверный. Существуют ли какие-либо библиотеки, методы, алгоритмы или решения для работы с таким количеством очень маленьких чисел?
Не стесняйтесь игнорировать все остальное.
У меня есть разреженная (~ 500 ненулевых элементов из 22500) матрица размером 150x150, представляющая систему линейных дифференциальных уравнений первого порядка. Я пытаюсь найти собственные значения и собственные векторы этой матрицы, чтобы построить функцию, которая служит аналитическим решением системы, так что я могу просто дать ей время, и она даст мне значения для каждой переменной. Я использовал этот метод в прошлом для аналогичных матриц 40x40, и он намного (в десятки, а в некоторых случаях в сотни раз) быстрее, чем scipy.integrate.solve_ivp()
, а также значительно упрощает постмодельный анализ, поскольку я могу найти максимальные значения и максимальные скорости изменения, используя scipy.optimize.fmin()
или оцените мою функцию в inf
, чтобы увидеть, как все уляжется, если оставить ее достаточно долго.
Однако на этот раз numpy.linalg.eig()
, похоже, не нравится моя матрица, и он дает мне комплексные значения, которые, как я знаю, неверны, потому что я моделирую физическую систему, которая не может иметь сложные скорости роста или распада (или синусоидальные решения). ), гораздо менее сложные значения для его переменных. Я считаю, что это проблема жесткости или округления с плавающей запятой, когда базовый алгоритм LAPACK не может обрабатывать либо очень маленькие значения (наименьшее значение ~ 3e-14, а большинство ненулевых значений имеют аналогичный масштаб), либо несоответствие между некоторыми значениями (самое большое составляет ~ 4000, но значения больше 1 появляются только несколько раз).
Я видел предложения для похожих проблем пользователей использовать sympy для решения собственных значений, но когда он не решил мою матрицу через 5 часов, я понял, что это не жизнеспособное решение для моей большой системы. Я также видел предложения использовать numpy.real_if_close()
для удаления мнимых частей комплексных значений, но я не уверен, что это тоже хорошее решение; несколько собственных значений из numpy.linalg.eig()
равны 0, что для меня является признаком ошибки, но, кроме того, почти все действительные части имеют тот же масштаб, что и мнимые части (чрезвычайно малы), что также заставляет меня сомневаться в их достоверности. Моя матрица реальная, но, к сожалению, несимметричная, поэтому numpy.linalg.eigh()
тоже нежизнеспособна.
Я нахожусь в точке, где я могу просто запускать scipy.integrate.solve_ivp()
в течение сколь угодно долгого времени (несколько тысяч часов), что, вероятно, потребует много времени для вычислений, а затем использовать scipy.optimize.curve_fit()
для аппроксимации аналитических решений, которые я хочу, поскольку у меня есть хороший представление об их формах. Это не идеально, так как делает мою программу намного медленнее, и я даже не уверен, что она будет работать с проблемами жесткости и округления, с которыми я столкнулся с numpy.linalg.eig()
; Я подозреваю, что Радау или BDF смогут справиться с жесткостью, но не с закруглением.
У кого-нибудь есть идеи? Любые другие алгоритмы поиска собственных значений, которые могли бы справиться с этим? Может ли numpy.linalg.eig()
работать с numpy.float128
вместо numpy.float64
или даже эта дополнительная точность не поможет?
Я рад предоставить дополнительную информацию по запросу. Я открыт для смены языка, если это необходимо.
Предоставьте достаточно кода, чтобы другие могли лучше понять или воспроизвести проблему.
@itprorh66 моя правка лучше? Честно говоря, я не уверен, что есть способ уточнить этот вопрос.
Учитывая природу чисел с плавающей запятой, я не думаю, что было бы проблемой решить матрицу для собственных значений с очень маленькими значениями. Однако потенциально это может быть проблемой, если у вас есть очень большие и очень маленькие значения вместе. Насколько вы уверены, что сложные решения неуместны? Многие физические системы демонстрируют синусоидальное поведение при некоторых условиях, хотя в нормальных условиях это синусоидальное поведение заглушается другими терминами.
numpy.linalg.eig()
производил комплексные числа для другого пользователя, а виновником, похоже, является округление маленьких чисел. Мои самые большие и самые маленькие числа указаны в исходном посте. Я уверен, что эти функции не являются периодическими; они моделируют перенос радиоактивных масс, а распады и движения экспоненциальны с реальными показателями.
Вы пытались использовать разреженные методы из scipy.sparse.linalg? В приведенном примере уже при построении ковариационной матрицы вносятся ошибки округления, делающие матрицу несимметричной и имеющей более высокий или полный ранг. Здесь с общим алгоритмом начальное преобразование к форме Хессенберга разрушит разреженную структуру и внесет случайные ошибки с плавающей запятой. // Кроме того, почему бы не использовать встроенные функции expm
?
Спасибо, @LutzLehmann! Это именно те ответы, которые я ищу. Я не был знаком с экспоненциальной матрицей как с решением для систем линейных ОДУ, но даже беглый взгляд на Википедию говорит мне, что ее стоит использовать здесь. Я также не был знаком с scipy.sparse.linalg()
и очень надеялся, что что-то подобное существует, но мой поиск в гугле ничего не дал. Я попробую оба и вернусь с результатами. Еще раз спасибо!
Добро пожаловать! ⚠️️🛒 Вопросы о рекомендациях по программному обеспечению здесь не по теме по причинам, описанным в № 3 раздела /help/on-topic.
Как упоминалось в цепочке комментариев выше, лучшим решением для этого является использование Матричной экспоненты, которая намного проще (и, по-видимому, менее подвержена ошибкам), чем диагонализация вашей системы с собственными векторами и собственными значениями.
В моем случае я использовал scipy.sparse.linalg.expm()
, так как моя система разрежена. Это быстро, точно и просто. Моя единственная жалоба - потеря оценки на бесконечности, но это достаточно легко обойти.
Добро пожаловать в Stack Overflow. ! Вопросы, которые требуют общего руководства относительно подхода к проблеме, как правило, не подходят для этого сайта. Разделите свою задачу на более мелкие и начните работать над ними. И задавайте конкретные вопросы, связанные с конкретными проблемами, с которыми вы сталкиваетесь. Вот как: как спросить