Optymalna metoda ODE dla stałej liczby ocen RHS

14

W praktyce czas wykonywania rozwiązywania liczbowego IVP

x˙(t)=f(t,x(t)) for t[t0,t1]
x(t0)=x0
jest często zdominowane przez czas oceny prawej strony (RHS)f . Załóżmy zatem, że wszystkie inne operacje są natychmiastowe (tj. Bez kosztów obliczeniowych). Jeżeli całkowity czas pracy dla rozwiązania ivp jest ograniczona, to jest równoznaczne z ograniczeniem liczby ocen odf do pewnegoNN .

Interesuje nas tylko końcowa wartość x(t1) .

Szukam teoretycznych i praktycznych wyników, które pomogą mi wybrać najlepszą metodę ODE w takim otoczeniu.

Jeśli na przykład , moglibyśmy rozwiązać IVP za pomocą dwóch wyraźnych kroków Eulera o szerokości ( t 1 - t 0 ) / 2 lub jednego kroku o szerokości t 1 - t 0 przy użyciu metody punktu środkowego. Nie jest od razu jasne, który z nich jest lepszy. W przypadku większego N można oczywiście pomyśleć o metodach wieloetapowych, iterowanych schematach Runge-Kutta itp.N=2(t1t0)/2t1t0N

Szukam wyników podobnych do tych, które istnieją, na przykład dla reguł kwadraturowych: Możemy wybrać wag { w i } i powiązane punkty { x i }, tak że reguła kwadraturowa n i = 1 w i g ( x i ) jest dokładne dla wszystkich wielomianów g tak, że d e g ( g ) 2 n - 1 .n{wi}{xi}i=1nwig(xi)gdeg(g)2n1

Dlatego szukam górnej lub dolnej granicy globalnej dokładności metod ODE, biorąc pod uwagę ograniczoną liczbę dozwolonych ocen RHS . Jest OK, jeśli granice obowiązują tylko dla niektórych klas RHS lub nakładają dodatkowe ograniczenia na rozwiązanie x (podobnie jak wynik dla reguły kwadratury, która obowiązuje tylko dla wielomianów do pewnego stopnia).fx

EDYCJA: Niektóre informacje podstawowe: dotyczy trudnych aplikacji w czasie rzeczywistym, tzn. Wynik musi być dostępny przed znanym terminem. Stąd ograniczenie liczby ocen RHS N jako dominującego czynnika kosztów. Zazwyczaj nasze problemy są sztywne i stosunkowo niewielkie.x(t1)N

EDIT2: Niestety nie mam dokładnych wymagań dotyczących czasu, ale można bezpiecznie założyć, że będzie raczej mały (zdecydowanie <100, prawdopodobnie bliższy 10). Biorąc pod uwagę wymagania w czasie rzeczywistym, musimy znaleźć kompromis między dokładnością modeli (przy lepszych modelach prowadzących do dłuższych czasów wykonania RHS, a tym samym do niższej N ) a dokładnością metody ODE (przy lepszych metodach wymagających wyższych wartości N ).NNN

Florian Brucker
źródło
Zwykła zgodność stałych metod Runge-Kutta z metodami Newtona-Cotesa dotyczy przypadku zastosowania metody RK do IVP ; np. zastosowanie klasycznej metody czwartego rzędu do tego IVP jest równoważne z wykonaniem reguły Simpsona na f ( x ) . y=f(x)f(x)
JM
@JM: Jestem tego świadomy. Chciałem wykorzystać reguły kwadraturowe jako przykład charakteryzujący dokładność metody numerycznej dla pewnego zestawu danych wejściowych, gdy liczba ocen funkcji jest ograniczona. Poza tym interesują mnie „prawdziwe” ODE, czyli takie, które nie ograniczają się do standardowej integracji.
Florian Brucker
1
To staje się coraz bardziej interesujące. Teraz sama liczba nic nie znaczy. Pomocna może być znajomość λ N / T , gdzie T jest długością przedziału całkowania, a λ jest stałą Lipschitza f względem x . Dzięki temu dowiemy się, jak poważny jest problem. Zakładając, że jest sztywny, prawdopodobnym kandydatem jest metoda BDF drugiego rzędu. NλN/TTλfx
David Ketcheson
@DavidKetcheson: Bardziej interesuje mnie ogólne podejście do wyboru odpowiedniej metody dla danego problemu niż optymalna metoda dla konkretnego problemu. Mamy większą liczbę modeli, które różnią się znacznie pod względem sztywności i wymagań dotyczących rozrządu.
Florian Brucker
Mówisz, że ocena jest bardzo droga. Czy potrafisz w ogóle obliczyć jakobian? Co powiesz na jakieś przybliżenie, które może skorygować podstawową sztywność? Nie jesteś w dobrej formie, jeśli twój problem jest bardzo sztywny i nie możesz go naprawić. f
Jed Brown

Odpowiedzi:

7

Myślę, że kluczowym odniesieniem do odpowiedzi na twoje pytanie jest ten artykuł autorstwa Ozeasza i Szampina . Teraz dam trochę tła.

Zasadniczo rozmiar kroku, którego można użyć podczas numerycznej integracji IVP, może być ograniczony przez stabilność lub dokładność. Jeśli chcesz wybrać najlepszy solver pod względem stabilności, musisz wziąć pod uwagę region absolutnej stabilności . W przypadku metody jednoetapowej jest to

S={zC:|P(z)|1}.

Tutaj jest funkcją stabilności metody; patrz np. tekst Hairer i in. glin. Niezbędnym warunkiem stabilności jest to, że λ h S, gdzie λ jest wyższe niż wartości własne jakobianu f i h, jest wielkością kroku. Nie zawsze jest to wystarczający warunek dla problemów nieliniowych, ale zwykle jest dobrą zasadą praktyczną i jest stosowany w praktyce.P(z)λhSλfh

Aby zapoznać się z obszernym podejściem do problemu znajdowania (jednoznacznych) metod, które pozwalają na uzyskanie dużych stabilnych rozmiarów kroków, zobacz ten mój artykuł na temat wielomianów stabilności i ten na temat optymalizacji metod Runge-Kutta dla symulacji płynów ściśliwych .

Stabilność jest istotnym problemem, jeśli okaże się, że największy stabilny rozmiar kroku już zapewnia wystarczającą dokładność. Z drugiej strony rozmiar kroku może zamiast tego być ograniczony twoimi wymaganiami dokładności. Zazwyczaj wykonuje się lokalną kontrolę błędów. Rozwiązanie oblicza się przy użyciu dwóch metod, a ich różnicę stosuje się jako oszacowanie błędu w mniej dokładnym. Rozmiar stopnia jest dobierany adaptacyjnie, aby osiągnąć jak najściślej zalecaną tolerancję.

Dwie miary teoretyczne są kluczowe do przewidywania wydajności dokładności. Pierwszy to rząd dokładności metody, który opisuje szybkość, z jaką błąd zbliża się do zera, gdy rozmiar kroku jest zmniejszany. Drugi to wskaźnik efektywności dokładności (patrz artykuł Ozeasza i Szampina połączony w pierwszym zdaniu powyżej), który uwzględnia stałe występujące w terminach błędów i umożliwia porównania między metodami tego samego rzędu.

Dokładność i efektywność stabilności szerokiej gamy metod można obliczyć w prosty i zautomatyzowany sposób przy użyciu NodePy (zastrzeżenie: NodePy został opracowany przeze mnie).

David Ketcheson
źródło
Dziękuję Ci. Artykuł Hosei i Shampine jest rzeczywiście bardzo interesujący. Czy znasz podobne wyniki dla sztywnych problemów? Wiem, że zwykle stosuje się dla nich metody niejawne, ale nie mają one a priori związku z liczbą ocen RHS, więc w moim przypadku są mało przydatne.
Florian Brucker
Nic takiego nie wiem z powodu sztywnych problemów, ale podejrzewam, że coś istnieje. Jak mówisz, pytanie jest bardziej subtelne, gdy używasz metod niejawnych. Jednym z podejść może być zastosowanie metod Rosenbrocka, które dobrze radzą sobie ze sztywnymi problemami, ale mają ustaloną liczbę ocen RHS.
David Ketcheson
6

Nie ma zbyt wielu wyników w tym kierunku, ponieważ jest to trudniejsze niż zwykłe ustalanie dokładności, ponieważ względy stabilności często mogą wymagać wybrania kroków czasowych, które są mniejsze niż byłoby to potrzebne dla pożądanej dokładności. Dlatego wyniki są podzielone między przypadki sztywne i niesztywne. W pierwszym przypadku wymagania dotyczące etapów czasowych i oceny RHS zasadniczo nie podlegają dokładności, aw drugim przypadku tak.

Skupię się na jawnych metodach, ponieważ domniemany przypadek jest o wiele mniej oczywisty, ile ocen RHS będziesz musiał użyć ... to całkowicie zależy od tego, jak zdecydujesz się rozwiązać wynikowy system.

W przypadku układów niesztywnych:

Istnieją wyraźne ograniczenia etapowe dla jednoznacznych metod Runge-Kutta, dla których powiedzieć, ile etapów (oceny RHS) jest potrzebnych do osiągnięcia określonej kolejności dokładności. Po czwartym rzędzie liczba etapów przekracza porządek dokładności, a różnica nadal rośnie. Duża książka ODE Butchera: http://books.google.com/books/about/Numerical_Methods_for_Ordinary_Different.html?id=opd2NkBmMxsC

wykonuje dobrą robotę, tłumacząc niektóre z tych dowodów „nieistnienia”.

Twój przykład reguły kwadratury prowadzi albo do metody wieloetapowej, takiej jak Adams-bashforth, lub do metod zwanych obecnie metodami korekcji odroczonej spektralnej. W przypadku Adamsa-bashfortha potrzebujesz tylko jednej oceny RHS na krok, ale ponieważ regiony stabilności są ogólnie tak małe dla tych metod, na ogół kończysz tyle samo pracy pod względem oceny RHS jak metoda Runge-Kutta z tym samym zamówienie.

Oto artykuł na temat odroczonej korekcji spektralnej:

https://www.google.com/search?q=spectral+deferred+correction&aq=f&oq=spectral+deferred+correction&aqs=chrome.0.57j0l2j62.3336j0&sourceid=chrome&ie=UTF-8

Nie jestem pewien, w jaki sposób te metody integracji sprawdzają się w porównaniu ze standardowymi metodami jawnymi, często wymagają one dużo więcej pamięci, aby zapisać stany rozwiązania w węzłach kwadraturowych, więc nigdy ich nie użyłem.

W przypadku sztywnych systemów:

istnieją „zoptymalizowane” narzędzia krokowe, ale dokładne wyniki teoretyczne dotyczące tego, jak dobre mogą być one, są niestety ograniczone do kilku prostych przypadków (a nawet te okazały się nie trywialne). Trzy standardowe wyniki mówią, że dla metod Runge-Kutta ze stopniami : Najbardziej ujemną osią rzeczywistą, jaką może ona zawierać w swoim obszarze stabilności, jest przedział długości 2 S 2 , najbardziej wyobrażoną osią, jaką może zawierać, jest przedział długości S - 1 , a największy okrąg styczny do wyobrażonej osi, który może zawierać, ma promień S (wszystkie one również wzajemnie się wykluczają).S2S2S1S

Reid.Atcheson
źródło
2
Może się zdarzyć, że zastosowanie metody kroku zmiennego (lub nawet zmiennej kolejności) może być bardziej wydajne niż uporczywe trzymanie się metody kroku stałego. Można na przykład rozważyć zastosowanie metody ekstrapolacyjnej, takiej jak Bulirsch-Stoer: wykonaj kilka ocen na niektórych etapach, a następnie zbuduj (pozornie) dokładniejsze szacunki na podstawie wyników tych kroków.
JM
Prawdziwe. W rzeczywistości wiele optymalnych metod jest w pewnym sensie równoważnych zmiennej wersji krokowej innego steppera czasowego. Na przykład Runge-Kutta-Chebshev może być postrzegany jako zastosowany Euler do przodu, przy czym zmienne stopnie czasowe są punktami Czebyszewa.
Reid.Atcheson
@JM: Dokładnie. Ale czy istnieje sposób oceny dokładności tych podejść w stosunku do liczby ocen RHS, oprócz eksperymentów numerycznych (które byłyby bardzo zaangażowane, biorąc pod uwagę dużą liczbę możliwych podejść)?
Florian Brucker
@Florian, nie ogólnie. Sądzę, że słyszałeś o równaniach Lorenza?
JM
1
@JM: Tak :) Dlatego wspomniałem o przykładzie kwadratury, w którym dokładność jest mierzona względem podzbioru (wielomianów) pierwotnej przestrzeni problemowej. Byłbym zadowolony z wyników, które działają tylko w przypadku określonego podzbioru problemów.
Florian Brucker
3

1014f(x)

Istnieją oczywiście wyjątki (bardzo duże systemy, bardzo sztywne systemy), ale powszechne jest przekonanie społeczności, że kwestia projektowania solverów ODE dla „standardowych” systemów jest rozwiązana. W związku z tym uważam, że pytanie, które stawiasz, nie jest bardzo interesujące - dotyczy elementu konstrukcji solvera ODE, który nie jest już ważny. Może to również tłumaczyć brak literatury na ten temat.

Wolfgang Bangerth
źródło
+1. Ilekroć ktoś pyta o wydajne solvera ODE, po prostu zakładam, że jest zainteresowany ogromnymi systemami ODE pochodzącymi z półdyskretyzacji PDE lub dużych problemów n-ciała.
David Ketcheson
Czy możesz wyjaśnić, w jaki sposób odnosi się to do mojego pytania? Nie widzę związku, ponieważ interesuje mnie przypadek, w którym ocena f(x)nie jest darmowa, ale jest tak droga, że ​​liczba ocen jest ograniczona.
Florian Brucker
@DavidKetcheson: W tym przypadku tak nie jest. Chodzi raczej o to, że mamy bardzo surowe wymagania dotyczące czasu (w czasie rzeczywistym) na słabym sprzęcie (urządzenia wbudowane). Same systemy ODE są stosunkowo małe.
Florian Brucker
NNNN.mniejszy błąd.
Wolfgang Bangerth
@FlorianBrucker: Myślę, że nie jest całkiem jasne, czy te same metody, które znamy, działają dobrze dla asymptotycznie dużych N. są również najlepsze dla małych N., mówić N.<1000. Innymi słowy, musisz powiedzieć nam coś o swoim budżecie, biorąc pod uwagę, ile wycen funkcji możesz sobie pozwolić.
Wolfgang Bangerth
1

Moje rozumowanie jest następujące. Ponieważ twoje problemy są sztywne i nie są duże, najprawdopodobniej powinieneś użyć jakiejś niejawnej metody. Jeśli zastosujesz metodę niejawną, rozwiążesz przynajmniej jeden układ liniowy na każdym kroku. Oznacza to, że jeśli jakobian nie ma odpowiedniej struktury, każdy krok zostanie wykonanyO(rejam3)) lub O(rejam2))klapy na każdym kroku (w zależności od tego, jak często LU rozkłada się jakobian), oprócz oceny RHS.

Pierwszą kwestią jest upewnienie się, czy RHS jest naprawdę droższy niż leżąca u podstaw algebra liniowa.

Drugi punkt: z literatury wiadomo, że solwery oparte na „drogich” metodach (tj. Jawnych metodach RK) czasami działają szybciej niż „tańsze” (jawne metody wieloetapowe).

Podsumowując, uważam, że należy nie tylko brać pod uwagę oceny RHS.

faleichik
źródło
Masz rację, że czas realizacji i dokładność to coś więcej niż liczba ocen RHS. Jednak w przypadku innych czynników (np. Algebry liniowej) jest dużo materiału na temat szybkości wykonania w porównaniu z rozmiarem problemu. Dlatego skupiam się naN.w tym pytaniu.
Florian Brucker