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

9

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.

wprowadź opis zdjęcia tutaj

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

Fomite
źródło
Czy potrafisz zapisać równania? Rozwiązywanie dużego, nieograniczonego układu równań nieliniowych często odbywa się numerycznie, stosując metodę Newtona lub jeden z jej wariantów. Wybór zależy od tego, ile masz informacji na temat oryginalnego układu równań - co najważniejsze, czy jakobian z układu równań jest dostępny, obliczalny lub łatwo przybliżony?
Aron Ahmadia
ach! Widzę, że twoje równania są wyszczególnione na stronie Mathematica. Czy masz coś przeciwko sprowadzeniu ich tutaj? (To nie jest post-cross-post, szczególnie jeśli zamierzamy zaproponować Ci rozwiązania numeryczne wykraczające poza zakres tego, co potrafi Mathematica).
Aron Ahmadia,
Później przyniosę równania z Mathematiki - po 5 godzinach jazdy muszę zejść z drogi.
Fomite,
1
Nie jest dUsdt=dHdt. Wydaje się, że tak jest z powyższych równań. Czy coś brakuje?
ja72
1
@GeoffOxberry: więc gdy pochodne są zrównane do zera, oba równania # 1 i # 2 są identyczne i jedno można pominąć.
ja72

Odpowiedzi:

13

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 (nth) root itp. Słynny przykład (w jednej zmiennej) to x5x+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);
Erik P.
źródło
7

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.

Arnold Neumaier
źródło
6

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

US 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 US pod względem wszystkich innych zmiennych, a równanie 2 można odrzucić.

US=HNθ(γ+ζ)CAKA+Cp+KD
gdzie KA=γρAσA+κρDσDτ+ρAσAζ i KD=γρpσp+κρDσD+ρpσpζ

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

D=κ(CAτ+Cp)γ+ζ.

Wynikowe wyrażenie sugeruje, że powinniśmy spróbować rozwiązać wszystkie pozostałe zmienne pod względem CA i CP; 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ć UA lub UP. Dodanie równań 4 i 6 razem daje również równanie liniowe we wszystkich zmiennych i można je rozwiązaćUA lub UP (które nie zostały rozwiązane podczas dodawania równań 3 i 5 razem).

W tym momencie powinniśmy mieć wyrażenia dla UA i UP pod względem H, CA, i CP (ponieważ możesz wyeliminować Duż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 CA i CP. Następnie, dokonując wszystkich niezbędnych podstawień, pozostałe równanie powinno być tylko pod względemCA i CP. Korzenie tego równania określą stany ustalone układu; symboliczne może być odnalezienie tych korzeni.

Powodzenia!

ja72
źródło
USnie pojawia się w żadnym równaniu innym niż równania 1 i 2. Zakładam, że wyeliminowałeś równania 1, 2 i 7. (Są to najłatwiejsze do rozwiązania.) Dodanie równań 3 i 5 daje równanie liniowe we wszystkich zmiennych, a zatem łatwe do rozwiązania. Podobnie dodanie równań 4 i 6 daje równanie liniowe we wszystkich zmiennych, a zatem łatwe do rozwiązania. To zajmuje 4 zmienne o wartości 7 (D, UA, UP, i US), aby wszystko było pod względem H, CA, i CP.
Geoff Oxberry
W tym momencie pozostały równania 3 i 4. (Równania 5 i 6 mają więcej terminów, więc wyrzućmy je.) Możesz użyć jednego z nich, aby rozwiązaćH pod względem CA i CP, i w tym momencie istnieje jedno równanie w odniesieniu do dwóch zmiennych: CA i CP, w którym to przypadku łatwiej byłoby znaleźć rozwiązania symboliczne ... chyba że całkowicie źle odczytam równania.
Geoff Oxberry
Tak jest. @GeoffOxberry, myślę, że powinieneś po prostu dodać swoje komentarze bezpośrednio do odpowiedzi ja72.
David Ketcheson
@DavidKetcheson: Gotowe; Nie martwię się o wikipikowanie, ponieważ przedstawiciel nie jest ważny. Nie wróciłem i nie wprowadziłem jeszcze symbolicznych manipulacji.
Geoff Oxberry
3

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ć F(x)=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

minxSF(x),

gdzie Sjest 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.
Geoff Oxberry
źródło
1

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

aterrel
źródło
Tak! Metody kontynuacji homotopii są wykładniczo trudne (trzeba to rozważyć2npoczątkowe warunki „rozruchu”), ale w przypadku tak małego problemu będzie on wykonalny obliczeniowo i możesz zagwarantować globalną optymalizację problemu minimalizacji.
Aron Ahmadia
Dr Ahmadia, najwyraźniej nie nadążałeś za literaturą na temat metod homotopii. Przeczytaj publikacje Jana i popraw ten numer.
aterrel