BDF vs niejawne przyspieszenie czasu Runge Kutta

Czy są jakieś powody, dla których należy wybierać domyślny Runge Kutta (IMRK) wyższego rzędu zamiast skokowego czasu BDF? BDF wydaje mi się znacznie łatwiejszy, ponieważ etap IMRK potrzebuje $q$$q$ rozwiązań liniowych na krok czasowy. Stabilność BDF i IMRK wydaje się być kwestią sporną. Nie mogę znaleźć żadnych zasobów porównujących / kontrastujących niejawne kroki krokowe.

Jeśli to pomoże, ostatecznym celem jest dla mnie wybór steppera czasowego domyślnego wysokiego rzędu dla PDE z dyfuzyjnym doradztwem.

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\tau$$\tau$

Metody BDF, będące metodami wieloetapowymi, mają najlepsze skalowanie

$f$radau 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) .

$f$ODEProblem

@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 ODEProblemmoż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, seulexktó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-10lub 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 MATLABode23smogą 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

$A$$Au$$e^A$$A$

$A$$Ju + g(u)$$J$$f = 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

$\Delta 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!]