Решение большой (150 переменных) системы линейных обыкновенных дифференциальных уравнений; столкнулись с проблемами округления и/или жесткости с плавающей запятой

Обновлено: Исходный пост слишком расплывчатый. Я ищу алгоритм для решения большой системы, решаемой, линейной 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 или даже эта дополнительная точность не поможет?

Я рад предоставить дополнительную информацию по запросу. Я открыт для смены языка, если это необходимо.

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

Community 06.01.2023 23:01

@itprorh66 моя правка лучше? Честно говоря, я не уверен, что есть способ уточнить этот вопрос.

Kibeth 07.01.2023 00:27

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

Nick ODell 07.01.2023 00:56
Вот менее экстремальный пример, когда numpy.linalg.eig() производил комплексные числа для другого пользователя, а виновником, похоже, является округление маленьких чисел. Мои самые большие и самые маленькие числа указаны в исходном посте. Я уверен, что эти функции не являются периодическими; они моделируют перенос радиоактивных масс, а распады и движения экспоненциальны с реальными показателями.
Kibeth 07.01.2023 02:01

Вы пытались использовать разреженные методы из scipy.sparse.linalg? В приведенном примере уже при построении ковариационной матрицы вносятся ошибки округления, делающие матрицу несимметричной и имеющей более высокий или полный ранг. Здесь с общим алгоритмом начальное преобразование к форме Хессенберга разрушит разреженную структуру и внесет случайные ошибки с плавающей запятой. // Кроме того, почему бы не использовать встроенные функции expm?

Lutz Lehmann 07.01.2023 08:32

Спасибо, @LutzLehmann! Это именно те ответы, которые я ищу. Я не был знаком с экспоненциальной матрицей как с решением для систем линейных ОДУ, но даже беглый взгляд на Википедию говорит мне, что ее стоит использовать здесь. Я также не был знаком с scipy.sparse.linalg() и очень надеялся, что что-то подобное существует, но мой поиск в гугле ничего не дал. Я попробую оба и вернусь с результатами. Еще раз спасибо!

Kibeth 07.01.2023 09:17

Добро пожаловать! ⚠️️🛒 Вопросы о рекомендациях по программному обеспечению здесь не по теме по причинам, описанным в № 3 раздела /help/on-topic.

user 12.01.2023 02:36
Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
0
8
79
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

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

В моем случае я использовал scipy.sparse.linalg.expm(), так как моя система разрежена. Это быстро, точно и просто. Моя единственная жалоба - потеря оценки на бесконечности, но это достаточно легко обойти.

Другие вопросы по теме