Символьное решение системы из 7 нелинейных уравнений

9

У меня есть система обыкновенных дифференциальных уравнений - 7 уравнений и ~ 30 параметров, управляющих их поведением, как часть математической модели передачи болезни. Я бы хотел , чтобы найти устойчивые состояния для этих уравнений Changing dx/dt = rest of the equationдля 0 = equationкаждого из уравнений делает простую проблему алгебры. Это можно сделать вручную, но я смехотворно плох в таких вычислениях.

Я пытался использовать Mathematica, которая может обрабатывать меньшие версии этой проблемы ( см. Здесь ), но Mathematica останавливает эту проблему. Есть ли более эффективный / эффективный способ приблизиться к этому? Более эффективная символическая математическая система? Другие предложения?

Несколько обновлений (21 марта):

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

введите описание изображения здесь

Таким образом, я ищу символические выражения для решений системы 7 квадратных уравнений в 7 переменных.

фомиты
источник
Вы можете записать уравнения? Решение большой неограниченной системы нелинейных уравнений часто выполняется численно с использованием метода Ньютона или одного из его вариантов. Выбор здесь будет зависеть от того, сколько информации у вас есть об исходной системе уравнений - наиболее важно, является ли якобиан системы уравнений доступным, вычислимым или легко аппроксимируемым?
Арон Ахмадиа
ааа! Я вижу, что ваши уравнения подробно изложены на сайте Mathematica. Ты не против привезти их сюда? (Это не перекрестная публикация, особенно если мы собираемся предложить вам численные решения, выходящие за рамки возможностей Mathematica).
Арон Ахмадиа
Я принесу уравнения из Mathematica позже сегодня - после 5-часовой поездки мне нужно убраться с дороги.
Fomite
1
Разве . Похоже, что так из приведенных выше уравнений. Я что-то пропустил? dUsdTзнак равно-dЧАСdT
ja72
1
@GeoffOxberry: поэтому, когда производные равны нулю, оба уравнения # 1 и # 2 идентичны, и одно можно опустить.
ja72

Ответы:

13

Похоже, что уравнения, с которыми вы имеете дело, являются полиномиальными после очистки знаменателей. Это хорошо (трансцендентные функции часто немного сложнее решать алгебраически). Однако это не гарантия того, что ваши уравнения имеют решение в замкнутой форме. Это важный момент, который многие люди на самом деле не «понимают», даже если они знают это теоретически, поэтому стоит повторить: существуют довольно простые системы полиномиальных уравнений, для которых нет способа дать решения в терминах (Nth) корни и т. д. Известный пример (по одной переменной) Икс5-Икс+1знак равно0, Смотрите также эту страницу википедии .

Сказав это, конечно, есть также системы уравнений, которые можно решить, и стоит проверить, является ли ваша система одной из них. И даже если ваша система не может быть решена, все же можно найти форму для вашей системы уравнений, которая в некотором смысле проще. Например, найдите одно уравнение, включающее только первую переменную (даже если оно не может быть решено алгебраически), затем второе уравнение, включающее только первую и вторую переменную и т. Д. Существует несколько конкурирующих теорий о том, как найти такие «нормальные формы» полиномиальных систем; наиболее известной является теория базиса Гребнера, а конкурирующей является теория регулярных цепей.

В системе компьютерной алгебры Maple (полное раскрытие: я работаю на них) оба они реализованы. solveКоманда , как правило , вызывает метод базисов Groebner, я считаю, что и быстро перемалывает к остановке на моем ноутбуке. Я попытался запустить обычное вычисление цепочек, и это занимает больше времени, чем у меня есть терпение, но, похоже, не так плохо, как память. Если вам интересно, страница справки для команды, которую я использовал, находится здесь , и вот код, который я использовал:

restart;
sys, vars := {theta*H - rho_p*sigma_p*
       Cp*(Us/N) - rho_d*sigma_d*D*(Us/N)*rho_a*sigma_a*
       Ca*(Us/N) = 0, 
         rho_p*sigma_p*Cp*(Us/N) + rho_d*sigma_d*
       D*(Us/N)*rho_a*sigma_a*Ca*(Us/N) + theta*H = 0, 
         (1/omega)*Ua - alpha*Up - rho_p*psi_p*
       Up*(H/N) - Mu_p*sigma_p*Up*(Cp/N) - 
             Mu_a*sigma_a*Up*(Ca/N) - Theta_p*
       Up + Nu_up*(Theta_*M + Zeta_*D) = 0, 
         alpha*Up - (1/omega)*Ua - rho_a*psi_a*
       Ua*(H/N) - Mu_p*sigma_p*Ua*(Cp/N) - 
             Mu_a*sigma_a*Ua*(Ca/N) - Theta_a*
       Ua + Nu_ua*(Theta_*M + Zeta_*D) = 0, 
         (1/omega)*Ca + Gamma_*Phi_*D + rho_p*psi_p*
       Up*(H/N) + Mu_p*sigma_p*Up*(Cp/N) + 
             Mu_a*sigma_a*Up*(Ca/N) - alpha*Cp - Kappa_*
       Cp - Theta_p*Cp + Nu_cp*(Theta_*M + Zeta_*D) = 0, 
         alpha*Cp + Gamma_*(1 - Phi_)*D + rho_a*psi_a*
       Ua*(H/N) + Mu_p*sigma_p*Ua*(Cp/N) + 
             Mu_a*sigma_a*Ua*(Ca/N) - (1/omega)*
       Ca - Kappa_*Tau_*Ca - Theta_a*Ca + 
             Nu_ca*(Theta_*M + Zeta_*D) = 
     0, Kappa_*Cp + Kappa_*Tau_*Ca - Gamma_*Phi_*
       D - Gamma_*(1 - Phi_)*D - 
             Zeta_*D + Nu_d*(Theta_*M + Zeta_*D) = 0, 
    Us + H + Up + Ua + Cp + Ca + D = 0, 
         Up + Ua + Cp + Ca + D = 0}, {Us, H, Up, Ua, Cp, Ca, D, N, 
    M}:

sys := subs(D = DD, sys):
vars := subs(D = DD, vars):
params := indets(sys, name) minus vars:
ineqs := [theta > 0 , rho_p > 0 , sigma_p > 
       0 , rho_d > 0 , sigma_d > 0 , 
            rho_a > 0 , sigma_a > 0 , 
      omega > 0 , alpha > 0 , psi_p > 0 , Mu_p > 0 , 
            Mu_a > 0 , Theta_p > 0 , Nu_up > 0 , Theta_ > 
       0 , Zeta_ > 0 , psi_a > 0 , 
            Theta_a > 0 , Nu_ua > 0 , Gamma_ > 0 , Phi_ > 
       0 , Kappa_ > 0 , Nu_cp > 0 , 
            Tau_ > 0 , Nu_ca > 0]:
with(RegularChains):
R := PolynomialRing([vars[], params[]]):
sys2 := map(numer, map(lhs - rhs, normal([sys[]]))):
sol := LazyRealTriangularize(sys2,[],map(rhs, ineqs),[],R);
Эрик П.
источник
7

Профессиональный способ - написать свои уравнения на языке моделирования, таком как AMPL или GAMS, и решить их с помощью решателя, такого как IPOPT.

AMPL - это коммерческая система, но бесплатная версия AMPL для студентов может создавать проблемы с 300 уравнениями и переменными.

Если вы просто хотите решить одну или несколько проблем, вы можете решить их в режиме онлайн бесплатно с помощью сервера NEOS для оптимизации - просто отправьте описание AMPL и подождите, пока ответ не будет возвращен вам.

Если вам необходимо повторно решить такие системы в рамках более масштабного исследования (например, для изменения параметров), вам следует загрузить IPOPT (программное обеспечение по очень либеральной лицензии).

Редактировать: Обратите внимание, что приемлемые символьные решения обычно ограничиваются довольно небольшими задачами - обычно размер базиса Гребнера растет с ростом числа переменных или степени полиномов, а время обработки еще больше. Таким образом, ожидание часа или более с Mathematica является признаком (хотя и не доказательством) того, что ваше символическое решение было бы совершенно непостижимым. Более того, оценка такого длинного выражения, вероятно, будет численно нестабильной, поэтому для получения значимых результатов вам потребуется высокая точность оценки.

Арнольд Ноймайер
источник
6

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

US не фигурирует ни в одном уравнении, кроме уравнений 1 и 2. Кроме того, эти уравнения являются зависимым множеством (уравнение 1 равно -1, умноженное на уравнение 2), поэтому уравнение 1 можно решить для US с точки зрения всех других переменных, и уравнение 2 можно отбросить.

USзнак равноЧАСNθ(γ+ζ)СAКA+Сп+КD
где КAзнак равноγρAσA+κρDσDτ+ρAσAζ а также КDзнак равноγρпσп+κρDσD+ρпσпζ

Уравнение 7 является линейным по всем переменным и может быть переставлено для решения D:

Dзнак равноκ(СAτ+Сп)γ+ζ,

Полученное выражение предполагает, что мы должны попытаться найти оставшиеся переменные в терминах СA а также Сп; поскольку у нас всего 6 независимых уравнений, лучшее, что мы можем сделать, - это свести систему к одному уравнению с двумя переменными.

К счастью, сложение уравнений 3 и 5 дает уравнение, которое является линейным по всем переменным и может быть решено для UA или Uп, Сложение уравнений 4 и 6 вместе также дает уравнение, которое является линейным по всем переменным и может быть решено дляUA или Uп (что не было решено при сложении уравнений 3 и 5 вместе).

На данный момент мы должны иметь выражения для UA а также Uп с точки зрения ЧАС, СA, а также Сп (потому что вы можете устранить Dиспользуя выражение выше). Мы использовали уравнения 1, 2, 5, 6 и 7; мы сохраним уравнения 3 и 4, потому что они проще.

Мы можем использовать уравнение 3 или 4 для решения ЧАС с точки зрения СA а также Сп, Тогда, делая все необходимые замены, оставшееся уравнение должно быть только в терминахСA а также Сп, Корни этого уравнения будут определять стационарные состояния системы; может быть или не быть возможно найти эти корни символически.

Удачи!

ja72
источник
USне фигурирует ни в одном уравнении, кроме уравнений 1 и 2. Я предполагаю, что вы исключили уравнения 1, 2 и 7. (Решить их проще всего). Добавление уравнений 3 и 5 дает уравнение, линейное по всем переменным, и, таким образом, легко решить. Точно так же добавление уравнений 4 и 6 дает вам уравнение, которое является линейным по всем переменным и, следовательно, легко решаемым. Так что заботится о 4 переменных из 7 (D, UA, Uп, а также US), так что все с точки зрения ЧАС, СA, а также Сп,
Джефф Оксберри
Таким образом, на данный момент у нас есть уравнения 3 и 4 слева. (Уравнения 5 и 6 имеют больше терминов, поэтому давайте выбросим их.) Вы можете использовать одно из них дляЧАС с точки зрения СA а также Спи в этот момент существует одно уравнение в терминах двух переменных: СA а также Сп, в этом случае может быть легче найти символические решения ... если я полностью не прочитал уравнения.
Джефф Оксберри
Право на. @ GeoffOxberry, я думаю, что вы должны просто добавить свои комментарии непосредственно к ответу ja72.
Дэвид Кетчесон
@DavidKetcheson: Готово; Я не беспокоюсь по поводу викификации, потому что репутация не важна. Я еще не вернулся и не заполнил символические манипуляции.
Джефф Оксберри
3

Это зависит от структуры ваших уравнений.

Если вы ищете все стационарные состояния вашего набора уравнений и можете перегруппировать их, как говорит ЭрикП, в полиномы, вы можете использовать методы из реальной алгебраической геометрии для вычисления всех численных решений с высокой точностью. Бертини - один из таких пакетов, о котором я знаю, но есть и другие. Несколько лет назад я отправился на конференцию в Нотр-Дам, где Бертини использовался для нахождения устойчивых состояний ОДУ из химической кинетики; Бертини был разработан в Нотр-Дам.

Другая возможность заключается в использовании методов, предложенных в «Негладком тесте исключения для нахождения всех решений нелинейных уравнений» М. Д. Стубером, В. Кумаром и П. И. Бартоном, BIT Численная математика 50 (4), 885-917, DOI: DOI: 10.1007 / s10543-010-0280-6 ; эти методы не требуют, чтобы система уравнений была полиномами. Пол Бартон - мой советник, а Мэтт Стубер - мой коллега; если хотите, я могу попросить у него программное обеспечение и отправить его вам. В работе используются методы из глобальной оптимизации и интервальной арифметики (в ней цитируется книга АрнольдНумаера), а также метод Ньютона. Преимущество этого метода заключается в том, что он должен найти все решения; недостаток в том, что это сложно.

В случае, если это не ясно, ArnoldNeumaier предлагает вместо решения F(Икс)знак равно0 непосредственно используя что-то вроде метода Ньютона (который обычно будет работать, если вы дадите ему хорошее начальное предположение, достаточно близкое к решению), вы решаете

минИксS| |F(Икс)| |,

где Sэто некоторый выполнимый набор, определенный ограничениями вашей проблемы вместо того, чтобы пытаться ее решить. На очень грубом уровне, использование гладкого решателя нелинейного программирования во многом похоже на использование метода Ньютона, с дополнительным алгоритмическим изощрением для надежности и производительности. IPOPT - действительно хорошее программное обеспечение для этой цели; Есть целый ряд других решателей (просто посмотрите на список доступных решателей для GAMS , AMPL или NEOS . Если вы выберете такой метод, помните о нескольких предостережениях:

  • За один раз будет найдено не более одного решения. Чтобы найти дополнительные решения, вам нужно добавить ограничения, которые исключают все предыдущие решения, которые вы нашли.
  • Если ваша задача оптимизации невыпуклая, для использования IPOPT или подобных решателей вам потребуется либо хорошее начальное предположение, близкое к решению ваших уравнений (тот же базовый принцип, что и метод Ньютона), либо невыпуклый решатель оптимизации, такой как BARON , Couenne , Bonmin и т. д. Вы должны попробовать каждый решатель, который у вас есть, так как производительность каждого невыпуклого решателя нелинейного программирования зависит от задачи.
Джефф Оксберри
источник
1

Я бы предложил посмотреть на гомотопический метод. Хотя это и не символично, оно даст все решения вашей проблемы. Для легкой библиотеки, чтобы проверить:

http://homepages.math.uic.edu/~jan/PHCpack/phcpack.html

aterrel
источник
Да! Методы гомотопического продолжения экспоненциально сложны (вам нужно будет рассмотреть2Nначальные «начальные» условия), но для такой маленькой задачи она будет вычислительно управляемой, и вы можете гарантировать глобальную оптимальность задачи минимизации.
Арон Ахмадиа
Доктор Ахмадиа, вы явно не поспеваете за литературой о гомотопических методах. Пожалуйста, прочитайте публикации Яна и исправьте это число.
Aterrel