Czy mogę użyć jawnego schematu krokowego w celu ustalenia liczbowego, czy ODE jest sztywny?

10

Mam ODE:

$u'=-1000u+sin(t)$
$u(0)=-\frac{1}{1000001}$

Wiem, że ten konkretny ODE jest sztywny, analitycznie. Wiem również, że jeśli użyjemy jawnej (krokowej) metody krokowej (Euler, Runge-Kutta, Adams itp.), Metoda powinna zwrócić bardzo duże błędy, jeśli krok czasowy jest zbyt duży. Mam więc dwa pytania:

Czy w ten sposób określa się sztywne wartości ODE na ogół, gdy wyrażenie analityczne dla terminu błędu jest niedostępne lub dające się wyprowadzić?
Ogólnie rzecz biorąc, kiedy ODE jest sztywny, w jaki sposób mogę określić „wystarczająco mały” przedział czasowy?

ode stiffness time-integration

— Paweł
źródło

Istnieją standardowe metody wykrywania sztywności za pomocą metod jawnych. Zamieszczam ten komentarz tutaj, ponieważ może być trudno znaleźć moją bardziej szczegółową odpowiedź daleko poniżej.

— David Ketcheson

6

Aby odpowiedzieć na twoje pytania:

O ile mi wiadomo, w praktyce, jeśli jawne metody wymagają wyjątkowo małych kroków czasowych w stosunku do twojej skali czasu zainteresowania (zobacz odpowiedzi na to pytanie, co oznacza sztywność ODE ), aby uzyskać dokładne wyniki, to dla wszystkie intencje i cele, twój problem jest sztywny. Aby określić wymagania dotyczące wielkości kroku, skorzystaj z jednej z wielu bibliotek napisanych przez ekspertów (pakiet MATLAB jest jednym z przykładów, także SUNDIALS, VODE, DASPK, DASSL, LSODE itp.), Które mają adaptacyjną heurystykę krokową. Podręcznik SUNDIALS wyjaśnia zasady decyzyjne, których używają do określenia wielkości kroków czasowych, które pakiet zajmuje, aby dać przykład reguł, które są stosowane w praktyce.
Ponownie użyłbym biblioteki z adaptacyjnym przyspieszaniem czasu w praktyce, ponieważ jest to bardziej wydajne. Jeśli jednak sam kodowałeś metodę, używając stałych rozmiarów kroków, jeśli zauważyłeś duże oscylacje lub „wysadzenie” rozwiązania, podejrzewasz, że twój krok czasowy był zbyt duży, i zredukuj go. Powtarzaj tę czynność, dopóki nie uzyskasz rozsądnego rozwiązania numerycznego. Podręczniki takie jak Ascher i Petzold oraz Hairer i Wanner mają dobre przykłady tego zjawiska.

— Geoff Oxberry
źródło

9

Lepszym sposobem na to jest to, że w przypadku sztywnego problemu każde stabilne, jawne obliczenie prowadzi do błędu, który jest znacznie mniejszy niż wymagana tolerancja błędu .

Istnieje wiele dobrych metod automatycznego wykrywania sztywności za pomocą jawnych schematów, zwłaszcza osadzonych par Runge-Kutta. Zobacz na przykład:

W drugim przykładzie faleichik, gdy wielkość kroku jest zmniejszona, można zauważyć nagły dramatyczny spadek błędu do poziomów znacznie poniżej typowej pożądanej tolerancji, gdy przekroczony zostanie próg stabilnego pomiaru czasu. Dobry estymator błędów rzeczywiście ujawniłby sztywność problemu. W pierwszym problemie błąd uzyskany przy stabilnym rozmiarze kroku mieściłby się w zakresie typowej pożądanej tolerancji, co wskazuje na brak sztywności.

W rezultacie zauważ, że każdy problem staje się niesztywny, jeśli wymagana jest wystarczająco ścisła tolerancja błędu.

— David Ketcheson
źródło

2

To były papiery, z którymi miałem zamiar się połączyć, zanim zobaczyłem twoją odpowiedź. +1 oczywiście. :) Pozwolę sobie również dodać to , to i wreszcie to . Jest to zdecydowanie dobrze zbadany problem ...

— JM

9

1. Czy możemy liczbowo wykryć sztywność, stosując tylko wyraźne metody?

Załóżmy, że masz problem z wartością początkową dla niektórych ODE na . Bierzesz znacznie duży stepize i jawną metodę Eulera, wykonujesz swoje obliczenia przy stałym rozmiarze kroku i otrzymujesz następujące punkty: $[0,10]$ $\tau=1$ $\tau$

Szacujesz błąd i wydaje się on duży. Okej, weź i uzyskaj $\tau=0.1$

Szacunkowy błąd jest teraz akceptowalny. Wielkość kroku jest niewielka w stosunku do . $\tau=0.1$ $[0,10]^\star$

Czy problem jest więc sztywny? Odpowiedź brzmi NIE ! W tym przypadku wymagany jest niewielki rozmiar, aby poprawnie odtworzyć oscylacje roztworu .

Problem, który rozwiązaliśmy, to
$y^{'} (t) = - 2 \cos π t, y (0) = 1.$ $y'(t)=-2\cos \pi t,\quad y(0)=1.$

Teraz bierzesz kolejne ODE w tym samym przedziale, a jawny Euler daje prawie takie samo rozwiązanie numeryczne: $\tau=1$

Bierzesz a teraz rozwiązaniem numerycznym jest $\tau=0.1$

Szacunkowy błąd jest teraz niewielki. Wielkość kroku jest niewielka w stosunku do . $\tau=0.1$ $[0,10]^\star$

Czy ten problem jest sztywny? TAK ! Zrobiliśmy bardzo małe kroki, aby odtworzyć rozwiązanie, które zmienia się bardzo powoli. To irracjonalne! Wielkość kroku czasowego jest tutaj ograniczona przez właściwości stabilności wyraźnego Eulera .

Ten problem to

$y^{'} (t) = - 2 y (t) + \sin t / 2, y (0) = 1.$ $y'(t)=-2 y(t)+\sin t/2,\quad y(0)=1.$

$^\star$ Zauważ, że liczba rozmiarów kroków może być znacznie większa, jeśli weźmiemy dłuższy interwał integracji.

Wniosek: informacja o krokach czasowych i odpowiadających im błędach nie jest wystarczająca do wykrycia sztywności. Powinieneś także spojrzeć na uzyskane rozwiązanie. Jeśli zmienia się powoli, a wielkość kroku jest bardzo mała, najprawdopodobniej problem będzie sztywny. Jeśli rozwiązanie oscyluje szybko i ufasz swojej technice szacowania błędów, problem nie jest sztywny.

2. Jak określić maksymalny rozmiar kroku, który pozwala zintegrować sztywny problem z jawną metodą?

Jeśli używasz jakiegoś jawnego solvera z czarną skrzynką z automatyczną kontrolą kroków, nie musisz nic robić: oprogramowanie dostosuje wymaganą wielkość kroków.

Załóżmy jednak, że chcesz ręcznie uzyskać rozwiązanie ze stałym krokiem lub po prostu chcesz oszacować, ile godzin powinieneś poczekać, aż Twoja jawna metoda rozwiąże problem. Następnie powinieneś znać spektrum swojej matrycy Jacobi. Załóżmy, że jest prawdziwy i leży w (w twoim przykładzie ). $[\Lambda,0]$ $\Lambda=-1000$

Następnie powinieneś obliczyć rzeczywisty przedział stabilności ( domena stabilności w złożonym przypadku) swojej jawnej metody. Nie jest to zbyt trudne, wystarczy zajrzeć do dowolnego podręcznika na ten temat. W przypadku jawnego Eulera odstęp ten wynosi . Teraz, jeśli chcesz, aby twoje rozwiązanie nie wybuchło, powinieneś tak, że leży w przedziale stabilności, tj. W naszym przypadku $[-2,0]$ $\tau$ $\Lambda\tau$

τ \leq \frac{2}{| Λ |} .

$\tau\leq \frac{2}{|\Lambda|}.$

Jeśli chcesz uzyskać większą spójność, powinieneś wziąć ponieważ dlaTwoje rozwiązanie najprawdopodobniej spowoduje nienaturalne (ale zanikające) oscylacje.

τ \leq \frac{1}{| Λ |},

$\tau\leq \frac{1}{|\Lambda|},$

1 / | Λ | < τ \leq 2 / | Λ |

$1/|\Lambda|<\tau \leq 2/|\Lambda|$

Oczywiście taka analiza ma zastosowanie głównie w przypadku problemów liniowych o znanym spektrum. W przypadku bardziej praktycznych problemów powinniśmy polegać na numerycznych metodach wykrywania sztywności (patrz referencje i komentarze w innych odpowiedziach).

— faleichik
źródło

Jak wspomniano w niektórych artykułach, do których David się przyłączył, metoda mocy do znajdowania dominujących wartości własnych (odpowiednio zmodyfikowanych) jest zwykłym wyborem dla wykrywaczy sztywności na bazie jakobianów.

— JM