Tak, z jakiegoś powodu nie ma na to zbyt wielu zasobów. Przez bardzo długi czas standardowym goto było „po prostu stosowanie metod BDF”. Ta mantra została ugruntowana z kilku historycznych powodów: z jednej strony kod Gear był pierwszym powszechnie dostępnym sztywnym solwerem, a z drugiej pakiet MATLAB nie zawierał / nie zawiera niejawnej metody RK. Jednak ta heurystyka nie zawsze jest poprawna i powiedziałbym, że z testowania jest zwykle źle. Pozwól mi wyjaśnić szczegółowo.
Metody BDF mają wysoki koszt stały
Adaptacyjny pomiar czasu i adaptacyjny porządek w metodach BDF ma naprawdę wysoki koszt. Współczynniki muszą zostać ponownie obliczone lub wartości należy interpolować do różnych czasów. Dużo pracy włożono w ulepszenie obecnych kodów BDF (istnieją dwie dobrze znane „formy” implementacji próbujące poradzić sobie z tym problemem inaczej), ale w rzeczywistości jest to po prostu bardzo trudny problem inżynierii oprogramowania. Z tego powodu właściwie wszystkie kody BDF są potomkami oryginalnego kodu Gear: Gear, vode, lsoda, Sundial's CVODE, nawet solwery DAE DASKR i DASSL są potomkami tego samego kodu.
Oznacza to, że jeśli problem jest „zbyt mały”, wysoki stały koszt naprawdę ma znaczenie, a metody Implicit RK będą działać lepiej.
Metody BDF wyższego rzędu są bardzo niestabilne w przypadku złożonych wartości własnych
Metody BDF pozwalają kontrolować maksymalne zamówienie i uczynić go bardziej konserwatywnym z jakiegoś powodu: metody BDF wyższego rzędu nie są w stanie bardzo dobrze poradzić sobie nawet ze „średnimi” złożonymi wartościami własnymi. Dlatego w tych przypadkach, aby być stabilnym, muszą porzucić swoje zamówienie. Jest to powód, dla którego metoda BDF szóstego rzędu, chociaż technicznie stabilna, jest często ignorowana, ponieważ nawet w piątym rzędzie występują już problemy (a szósty rząd jest jeszcze mniej stabilny). Tylko do drugiego rzędu jest stabilny na A, więc zawsze może do niego wrócić, ale wtedy krok jest zdominowany przez błąd.
Tak więc, gdy używasz kodów BDF do nieistotnych problemów, nie dostajesz piątego zamówienia przez cały czas. Oscylacje powodują spadek ich kolejności.
Metody BDF mają wysoki koszt początkowy
i ττ
Metody BDF, będące metodami wieloetapowymi, mają najlepsze skalowanie
faradau
metodami do 1000 zmiennych przestrzennych, więc YMMV. Jedynym sposobem na sprawdzenie jest przetestowanie.
Jakie testy porównawcze są dostępne?
Książka Hairera i DiffEqBenchmark (wyjaśnione poniżej) są prawdopodobnie najlepsze pod względem łatwo dostępnych diagramów precyzji pracy. Rozwiązywanie zwykłych równań różniczkowych różnicowych Hairera II zawiera szereg diagramów precyzji pracy na stronach 154 i 155. Wyniki wybranych przez niego problemów są zgodne z tym, co wskazałem powyżej z powodów, które podałem powyżej: domyślna RK będzie bardziej skuteczna, jeśli problem nie jest "wystarczająco duży". Inną ciekawą rzeczą do odnotowania jest to, że metody Rosenbrocka wyższego rzędu okazują się najbardziej wydajne w wielu jego testach (jak Rodas
) w reżimie, w którym błąd jest wyższy, i niejawne RKradau5
jest najbardziej wydajna przy niższych błędach. Ale jego problemami z testami nie są w większości dyskretyzacje dużych PDE, więc weź pod uwagę powyższe punkty.
Jak testujesz / testujesz?
Lubię to testować w DifferentialEquations.jl w Julii (zastrzeżenie: jestem jednym z programistów). To jest w Julii. Język programowania powinien dostać tutaj notatkę. Pamiętaj, że wraz ze wzrostem kosztu wywołania funkcji metody BDF są coraz lepsze. W R / MATLAB / Python funkcja użytkownika jest jedynym kodem R / MATLAB / Python, którego faktycznie muszą używać zoptymalizowane solwery: jeśli używasz opakowań SciPy lub Sundials, jest to cały kod C / Fortran z wyjątkiem funkcji, którą przekazujesz . Oznacza to, że w dynamicznych językach (które nie są Julią) metody BDF radzą sobie lepiej niż zwykle, ponieważ wywołania funkcji są bardzo niezoptymalizowane (prawdopodobnie jest to powód włączenia Shampineode15s
w pakiecie MATLAB, ale nie ma domyślnej metody RK) .
faODEProblem
@time sol = solve(prob,CVODE_BDF())
@time sol = solve(prob,radau())
gdzie pierwsza wykorzystuje Sundials CVODE
(metoda BDF), a druga używa Hairera radau
(domyślna RK).
W każdym przypadku ODEProblem
możesz użyć narzędzi do analizy porównawczej, aby zobaczyć, jak różne algorytmy skalują się dla różnych tolerancji adaptacyjnych. Niektóre wyniki są publikowane na stronie DiffEqBenchmarks.jl . Na przykład w przypadku problemu ROBER (system 3 sztywnych ODE) widać, że zegar słoneczny jest niestabilny i różni się wystarczająco wysoką tolerancją (podczas gdy inne metody są w porządku), co pokazuje powyższą uwagę na temat problemów ze stabilnością. W przypadku problemu Van Der Pol widać, że jest to bardziej pranie. Mam o wiele więcej niż nie opublikowałem, ale prawdopodobnie nie dostanę się do tego, dopóki nie skończę metod Rosenbrock wyższego rzędu (Rodas jest ich wersją Fortran).
(Uwaga: te sztywne testy porównawcze wymagają aktualizacji. Po pierwsze, tekst wymaga aktualizacji, ponieważ z jakiegoś powodu metody ODE.jl różnią się ...)
Metody ekstrapolacji i RKC
Istnieją również metody ekstrapolacji, seulex
które są stworzone dla sztywnych problemów. Są to „nieskończone porządki adaptacyjne”, ale to oznacza tylko, że są one asymetrycznie dobre, gdy szukasz bardzo niskiego błędu (tj. Szukasz rozwiązania sztywnych problemów na mniejszym 1e-10
lub większym poziomie , w którym to przypadku prawdopodobnie możesz użyć jawnej metody) . Jednak w większości przypadków nie są one tak wydajne i należy ich unikać.
Runge-Kutta Chebyschev to wyraźna metoda, która działa również na sztywne problemy, które należy wziąć pod uwagę. Nie mam go jeszcze w pliku DifferentialEquations.jl, więc nie mam żadnych twardych dowodów na to, jak to się dzieje.
Inne metody do rozważenia: specjalistyczne metody dla sztywnych PDE
Warto zauważyć, że metody Rosenbrocka wyższego rzędu radzą sobie naprawdę dobrze, wiele razy nawet lepiej, w przypadku sztywnych problemów małych i średnich rozmiarów, kiedy można łatwo obliczyć jakobian. Jednak w przypadku niektórych PDE uważam, że problemy z rozpowszechnianiem porad należą do tej kategorii, Rosenbrock może stracić pewne rzędy dokładności. Potrzebują również bardzo dokładnych jakobianów, aby zachować swoją dokładność. W Julii jest to łatwe, ponieważ solwery są wyposażone w symboliczne i automatyczne rozróżnianie, które może być poprawne dla maszyny epsilon. Jednak rzeczy takie jak MATLABode23s
mogą mieć problemy, ponieważ te implementacje używają różnicowania skończonego. W przypadku BDF i niejawnych metod RK błędy w Jakubie prowadzą do wolniejszej zbieżności, natomiast w przypadku Rosenbrocka, ponieważ nie są to równania niejawne i są raczej metodami RK z odwróconymi jakobianami, po prostu tracą porządek dokładności.
Inne metody, na które należy spojrzeć, to wykładnicze metody RK, wykładnicze różnicowanie czasu (ETD), ukryty współczynnik integracji (IIF) i wykładnicze metody Rosenbrocka. Pierwsze trzy wykorzystują fakt, że w wielu dyskretyzacjach PDE
ut= A u + f( u )
ZAA umiZAZA
ZAjotu +g( u )jotfa= Ju + g
Jeszcze inne metody: najnowsze kreacje
Metody, które są w pełni niejawne, oczywiście sprawdzają się w przypadku równań sztywnych. Jeśli PDE nie jest wystarczająco duże, nie można skutecznie „zrównoleglić się w przestrzeni” wystarczającej do wykorzystania HPC. Zamiast tego możesz tworzyć dyskretyzacje równolegle w czasie, które są w pełni niejawne, a zatem dobre dla sztywnych równań, a jednocześnie równoległe, aby w pełni wykorzystać sprzęt. XBraid to solver, który wykorzystuje taką technikę, a głównymi metodami są PFFAST i metody pomocnicze. DifferentialEquations.jl rozwija metodę sieci neuronowej, która działa podobnie.
Znowu jest to optymalne, gdy nie masz wystarczająco dużej dyskretyzacji przestrzennej, aby skutecznie zrównoleglać, ale masz dostępne zasoby do obliczeń równoległych.
Wniosek: rozważ asymptotyczne rozważania z odrobiną soli
Δ t
Metody BDF są asymptotycznie najlepsze, ale w większości przypadków prawdopodobnie nie pracujesz w tym systemie. Ale jeśli dyskretyzacja przestrzenna ma wystarczającą liczbę punktów, metody BDF mogą skutecznie paralelizować w przestrzeni (przez równoległe rozwiązywanie liniowe) i będą miały najmniej wywołań funkcji, a tym samym zrobią to najlepiej. Ale jeśli dyskretyzacja PDE nie jest wystarczająco duża, przyjrzyj się niejawnym metodom RK, Rosenbrock, wykładniczym RK itp. Korzystanie z pakietu oprogramowania, takiego jak DifferentialEquations.jl, który ułatwia przełączanie się między różnymi metodami, może być naprawdę pomocne w zrozumieniu, która metoda najlepiej odpowiada problematycznej domenie, ponieważ w wielu przypadkach nie można jej wcześniej ustalić.
[Jeśli masz jakieś przykładowe problemy, które chcesz dodać do pakietu testowego, nie wahaj się pomóc. Chcę mieć bardzo obszerne zasoby na ten temat. Mam nadzieję, że wszystkie wzorce Hairer w uruchamianych postaciach notebooków będą „wkrótce” i chciałbym, aby inne problemy były reprezentatywne dla dziedzin naukowych. Każda pomoc jest mile widziana!]