Optymalna metoda ODE dla stałej liczby ocen RHS

W praktyce czas wykonywania rozwiązywania liczbowego IVP

$$\dot{x}(t) = f(t, x(t)) \quad \text{ for } t \in [t_0, t_1]$$ $$x(t_0) = x_0$$ 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 od$f$ do pewnego$N \in \mathbb{N}$ .

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

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$$(t_1 - t_0)/2$$t_1 - t_0$$N$

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$$\{w_i\}$$\{x_i\}$$\sum_{i = 1}^n w_i g(x_i)$$g$$deg(g) \le 2n - 1$

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

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(t_1)$$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 ).$N$$N$$N$

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 = \{ z \in {\mathbb C} : |P(z)|\le 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)$$\lambda h \in S$$\lambda$$f$$h$

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

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ą).$S$$2S^2$$S-1$$S$

$10^{-14}$$f(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.

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 wykonany$O(\mathit{dim}^3)$ lub $O(\mathit{dim}^2)$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.