Symboliczne rozwiązanie układu 7 równań nieliniowych

Mam układ zwykłych równań różniczkowych - 7 równań i ~ 30 parametrów rządzących ich zachowaniem w ramach matematycznego modelu przenoszenia choroby. Chciałbym jak znaleźć stabilne stany tych równań Zmiana dx/dt = rest of the equationaby 0 = equationdla każdego z równań sprawia, że prosta algebra problem. Można to zrobić ręcznie, ale jestem absurdalnie zły w tego rodzaju obliczeniach.

Próbowałem użyć Mathematiki, która może obsługiwać mniejsze wersje tego problemu ( patrz tutaj ), ale Mathematica zatrzymuje się na ten problem. Czy istnieje bardziej wydajny / skuteczny sposób podejścia do tego? Bardziej wydajny symboliczny system matematyczny? Inne sugestie?

Kilka aktualizacji (21 marca):

• Celem jest rzeczywiście ich rozwiązanie symboliczne - odpowiedzi liczbowe są fajne, ale w tej chwili celem końcowym jest wersja symboliczna.
• Jest co najmniej jedna równowaga. Właściwie to nie usiadłem i nie udowodniłem tego, ale z założenia powinien on mieć co najmniej jeden trywialny, w którym żadne nie jest zainfekowane na początku. Poza tym może nie być nic , ale to sprawiłoby, że jestem równie zadowolony jak wszystko inne.
• Poniżej znajduje się rzeczywisty zestaw równań, o których mowa.

Podsumowując, szukam wyrażeń symbolicznych dla rozwiązań układu 7 równań kwadratowych w 7 zmiennych.

Wygląda na to, że równania, z którymi masz do czynienia, są wielomianowe po usunięciu mianowników. To dobrze (funkcje transcendentalne są często nieco trudniejsze do algebraicznie). Nie jest to jednak gwarancją, że twoje równania mają rozwiązanie w formie zamkniętej. Jest to istotny punkt, którego wielu ludzi tak naprawdę „nie rozumie”, nawet jeśli znają to w teorii, więc trzeba to powtórzyć: istnieją dość proste układy równań wielomianowych, dla których nie ma możliwości podania rozwiązań w kategoriach ($n$th) root itp. Słynny przykład (w jednej zmiennej) to $x^5-x+1=0$. Zobacz także tę stronę wikipedii .

To powiedziawszy, oczywiście istnieją również układy równań, które można rozwiązać, i warto sprawdzić, czy twój układ jest jednym z nich. I nawet jeśli twojego układu nie da się rozwiązać, nadal może być możliwe znalezienie formy układu równań, która jest w pewnym sensie prostsza. Na przykład znajdź jedno równanie obejmujące tylko pierwszą zmienną (nawet jeśli nie można go rozwiązać algebraicznie), a następnie drugie równanie obejmujące tylko pierwszą i drugą zmienną itp. Istnieje kilka konkurencyjnych teorii na temat znajdowania takich „normalnych form” układów wielomianowych; najbardziej znana jest teoria podstawowa Groebnera, a konkurencyjna to teoria regularnych łańcuchów.

W systemie algebry komputerowej Maple (pełne ujawnienie: pracuję dla nich) oba są zaimplementowane. Wydaje mi się, że to solvepolecenie zwykle wywołuje metodę Groebnera, co szybko zatrzymuje się na moim laptopie. Próbowałem uruchomić regularne obliczanie łańcuchów i trwa to dłużej niż mam cierpliwość, ale wydaje się, że nie jest tak kiepskie pod względem pamięci. Jeśli jesteś zainteresowany, strona pomocy dla polecenia, którego użyłem, jest tutaj , a oto kod, którego użyłem:

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);

Profesjonalnym sposobem jest napisanie równań w języku modelowania, takim jak AMPL lub GAMS, i rozwiązanie ich za pomocą solwera, takiego jak IPOPT.

AMPL to system komercyjny, ale darmowa wersja studencka AMPL może stwarzać problemy z maksymalnie 300 równaniami i zmiennymi.

Jeśli chcesz tylko rozwiązać jeden lub kilka problemów, możesz go rozwiązać online za pomocą serwera NEOS do optymalizacji - po prostu prześlij opis AMPL i poczekaj na odpowiedź.

Jeśli musisz wielokrotnie rozwiązywać takie systemy w ramach większego badania (np. Zmieniając parametry), powinieneś pobrać IPOPT (oprogramowanie na bardzo liberalnej licencji).

Edycja: Zwróć uwagę, że zrozumiałe rozwiązania symboliczne są zwykle ograniczone do dość małych problemów - zwykle rozmiar podstawy Groebnera rośnie gwałtownie wraz z liczbą zmiennych lub stopniem wielomianów, a czas przetwarzania jest jeszcze większy. Zatem czas oczekiwania wynoszący godzinę lub dłużej z Mathematica jest znakiem (choć nie dowodem), że twoje symboliczne rozwiązanie byłoby całkowicie niezrozumiałe. Co więcej, ocena tak długiego wyrażenia może być niestabilna numerycznie, więc aby uzyskać znaczące wyniki, potrzebna byłaby wysoka precyzja oceny.

Napisanie całego rozwiązania jest niemożliwe z uzasadnionego powodu. Ale oto kilka równań, aby nieco zmniejszyć system:

$U_{S}$ nie pojawia się w żadnym równaniu innym niż równanie 1 i 2. Ponadto równania te są zbiorem zależnym (równanie 1 to -1 razy równanie 2), więc równanie 1 można rozwiązać dla $U_{S}$ pod względem wszystkich innych zmiennych, a równanie 2 można odrzucić.

Równanie 7 jest liniowe we wszystkich zmiennych i można je zmienić w celu rozwiązania $D$:

Wynikowe wyrażenie sugeruje, że powinniśmy spróbować rozwiązać wszystkie pozostałe zmienne pod względem $C_{A}$ i $C_{P}$; ponieważ mamy tylko 6 niezależnych równań, najlepsze, co możemy zrobić, to sprowadzić układ do jednego równania w dwóch zmiennych.

Na szczęście dodanie równań 3 i 5 razem daje równanie liniowe we wszystkich zmiennych i można je rozwiązać $U_{A}$ lub $U_{P}$. Dodanie równań 4 i 6 razem daje również równanie liniowe we wszystkich zmiennych i można je rozwiązać$U_{A}$ lub $U_{P}$ (które nie zostały rozwiązane podczas dodawania równań 3 i 5 razem).

W tym momencie powinniśmy mieć wyrażenia dla $U_{A}$ i $U_{P}$ pod względem $H$, $C_{A}$, i $C_{P}$ (ponieważ możesz wyeliminować $D$używając powyższego wyrażenia). Użyliśmy równań 1, 2, 5, 6 i 7; zachowamy równania 3 i 4, ponieważ są prostsze.

Do rozwiązania możemy użyć równania 3 lub 4 $H$ pod względem $C_{A}$ i $C_{P}$. Następnie, dokonując wszystkich niezbędnych podstawień, pozostałe równanie powinno być tylko pod względem$C_{A}$ i $C_{P}$. Korzenie tego równania określą stany ustalone układu; symboliczne może być odnalezienie tych korzeni.

Powodzenia!

To zależy od struktury twoich równań.

Jeśli szukasz wszystkich stałych stanów z zestawu równań i możesz zmienić ich kolejność, jak ErikP mówi na wielomiany, możesz użyć metod z prawdziwej geometrii algebraicznej, aby obliczyć wszystkie rozwiązania numeryczne z wysoką precyzją. Bertini jest jednym z takich pakietów, które znam, ale są też inne. Kilka lat temu poszedłem na konferencję w Notre Dame, gdzie Bertini był używany do znajdowania stałych stanów ODE na podstawie kinetyki chemicznej; Bertini został opracowany w Notre Dame.

Inną możliwością jest zastosowanie metod zaproponowanych w „Teście wykluczenia niehoothowego do znajdowania wszystkich rozwiązań równań nieliniowych” autorstwa MD Stubera, V. Kumara i PI Bartona, BIT Mathematics 50 (4), 885-917, DOI: DOI: 10.1007 / s10543-010-0280-6 ; metody te nie wymagają, aby układ równań był wielomianami. Paul Barton jest moim doradcą, a Matt Stuber jest moim kolegą; jeśli chcesz, mogę poprosić go o oprogramowanie i wysłać je do ciebie. W pracy wykorzystano metody z globalnej optymalizacji i arytmetyki przedziałowej (cytuje książkę Arnolda Neumaiera), a także metodę Newtona. Zaletą tej metody jest to, że powinna ona zlokalizować wszystkie rozwiązania; wadą jest to, że jest skomplikowane.

W przypadku, gdy nie jest to jasne, ArnoldNeumaier sugeruje, aby zamiast rozwiązać $\mathbf{F}(\mathbf{x}) = \mathbf{0}$ bezpośrednio używając czegoś takiego jak metoda Newtona (która na ogół zadziała, jeśli dobrze zgadniesz, wystarczająco blisko rozwiązania), rozwiązujesz

gdzie $S$jest jakimś możliwym zestawem zdefiniowanym przez ograniczenia twojego problemu zamiast próbować go rozwiązać. Na bardzo prymitywnym poziomie korzystanie z płynnego nieliniowego solvera programistycznego przypomina korzystanie z metody Newtona, z dodatkowym algorytmicznym wyrafinowaniem zapewniającym niezawodność i wydajność. IPOPT jest naprawdę dobrym oprogramowaniem do tego celu; istnieje litania innych solverów (wystarczy spojrzeć na listę dostępnych solverów dla GAMS , AMPL lub NEOS . Jeśli wybierzesz taką metodę, pamiętaj o kilku zastrzeżeniach:

• Lokalizuje tylko najwyżej jedno rozwiązanie na raz. Aby znaleźć dodatkowe rozwiązania, musisz dodać ograniczenia wykluczające wszystkie poprzednie znalezione rozwiązania.
• Jeśli twój problem optymalizacji nie jest wypukły, aby użyć IPOPT lub podobnych solverów, będziesz musiał albo dobrze zgadnąć początkowo, blisko rozwiązania swoich równań (ta sama podstawowa zasada, co metoda Newtona), lub niekonwersyjny solver optymalizacji, taki jak BARON , Couenne , Bonmin itp. Powinieneś wypróbować każdy solver, na który masz ochotę, ponieważ wydajność każdego nie wypukłego nieliniowego solvera programistycznego zależy od problemu.

Sugerowałbym spojrzenie na metodę homotopii. Chociaż nie jest to symboliczne, stworzy wszystkie rozwiązania twojego problemu. Aby łatwo sprawdzić bibliotekę:

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