Adaptiv stegstorlek

Inom matematik och numerisk analys används en adaptiv stegstorlek i vissa metoder för numerisk lösning av vanliga differentialekvationer (inklusive specialfallet med numerisk integration ) för att kontrollera metodens fel och för att säkerställa stabilitetsegenskaper som A- stabilitet . Att använda en adaptiv stegstorlek är särskilt viktigt när det finns en stor variation i storleken på derivatan. Till exempel, när man modellerar en satellits rörelse runt jorden som en vanlig Kepler-bana , kan en fast tidsstegningsmetod som Euler-metoden vara tillräcklig. Men saker och ting är svårare om man vill modellera rörelsen hos en rymdfarkost med hänsyn till både jorden och månen som i trekroppsproblemet . Där dyker det upp scenarier där man kan ta stora tidssteg när rymdfarkosten är långt från jorden och månen, men om farkosten kommer nära att kollidera med en av planetkropparna så behövs små tidssteg. Rombergs metod och Runge–Kutta–Fehlberg är exempel på numeriska integrationsmetoder som använder en adaptiv stegstorlek.

Exempel

För enkelhetens skull använder följande exempel den enklaste integrationsmetoden, Euler-metoden ; i praktiken föredras metoder av högre ordning som Runge–Kutta -metoder på grund av deras överlägsna konvergens- och stabilitetsegenskaper.

Tänk på problemet med initialvärde

$y'(t)=f(t,y(t)),\qquad y(a)=y_ {a}$

där y och f kan beteckna vektorer (i vilket fall denna ekvation representerar ett system av kopplade ODEs i flera variabler).

Vi får funktionen f ( t , y ) och initialvillkoren ( a , y _a ), och vi är intresserade av att hitta lösningen vid t = b . Låt y ( b ) beteckna den exakta lösningen vid b , och låt y _b beteckna den lösning som vi beräknar. Vi skriver ${\displaystyle y_{b}+\varepsilon =y(b)} ,$ där $\varepsilon$ är felet i den numeriska lösningen.

För en sekvens ( t _n ) av värden på t , med t _n = a + nh , ger Euler - metoden approximationer till motsvarande värden på y ( t _n ) som

$y_{n+1}^{(0)}=y_{n}+hf(t_{n},y_{n })$

Det lokala trunkeringsfelet för denna approximation definieras av

$\tau _{n+1}^{(0)}=y(t_{n+1})-y_ {n+1}^{(0)}$

och med Taylors teorem kan det visas att (förutsatt att f är tillräckligt jämnt) det lokala trunkeringsfelet är proportionellt mot kvadraten på stegstorleken:

$\tau _{n+1}^{(0)}=ch^{2}$

där c är någon proportionalitetskonstant.

Vi har markerat den här lösningen och dess fel med en $(0)$ .

Värdet på c är inte känt för oss. Låt oss nu tillämpa Eulers metod igen med en annan stegstorlek för att generera en andra approximation till y ( t _{n +1} ). Vi får en andra lösning, som vi märker med en $(1)$ . Anta att den nya stegstorleken är hälften av den ursprungliga stegstorleken och tillämpa två steg av Eulers metod. Denna andra lösning är förmodligen mer exakt. Eftersom vi måste tillämpa Eulers metod två gånger är det lokala felet (i värsta fall) två gånger det ursprungliga felet.

$y_{n+{\frac {1}{2}}}=y_{n}+{\frac {h}{2 }}f(t_{n},y_{n})$

$y_{ n+1}^{(1)}=y_{n+{\frac {1}{2}}}+{\frac {h}{2}}f(t_{n+{\frac {1}{2} }},y_{n+{\frac {1}{2}}})$

$\tau _{n+1}^{(1)}=c\left({\frac {h}{2}}\right)^{ 2}+c\left({\frac {h}{2}}\right)^{2}=2c\left({\frac {h}{2}}\right)^{2}={\frac {1}{2}}ch^{2}={\frac {1}{2}}\tau _{n+1}^{(0)}$

$y_{n+1}^{(1)}+\tau _{n+1}^{(1)}=y(t+h)$

Här antar vi att felfaktorn $c$ är konstant över intervallet $[t,t+h]$ . I verkligheten är dess förändringshastighet proportionell mot $y^{(3)}(t)$ . Att subtrahera lösningar ger feluppskattningen:

$y_{n+1}^{(1)}-y_{n+1}^{(0)} =\tau _{n+1}^{(1)}$

Denna lokala feluppskattning är tredje ordningens korrekt.

Den lokala feluppskattningen kan användas för att bestämma hur stegstorlek $h$ ska modifieras för att uppnå önskad noggrannhet. Till exempel, om en lokal tolerans av ${\text{tol}}$ är tillåten, kan vi låta h utvecklas som:

$h\rightarrow 0,9\times h\times \min \left(\ max \left(\left({\frac {\text{tol}}{2\left|\tau _{n+1}^{(1)}\right|}}\right)^{1/2} ,0.3\right),2\right)$

0.9 $0.9$ är en säkerhetsfaktor för att säkerställa framgång vid nästa försök Minsta och högsta är för att förhindra extrema förändringar från föregående stegstorlek. Detta bör i princip ge ett fel på cirka $0,9\ gånger {\text{tol}}$ i nästa försök. Om ${\displaystyle |\tau _{n+1}^{(1)}|<{\text{tol}}} ,$ vi anser att steget är lyckat och feluppskattningen används för att förbättra lösningen:

$y_{n+1}^{(2)}=y_{n+1}^{(1) }+\tau _{n+1}^{(1)}$

Denna lösning är faktiskt tredje ordningens korrekt i det lokala omfånget (andra ordningen i det globala omfattningen), men eftersom det inte finns någon feluppskattning för den, hjälper det inte till att minska antalet steg. Denna teknik kallas Richardson extrapolation .

Börjar med en initial stegstorlek på $h=ba$ , denna teori underlättar vår kontrollerbara integration av ODE från punkt $a$ till $b$ , med ett optimalt antal steg ges en lokal feltolerans. En nackdel är att stegstorleken kan bli oöverkomligt liten, speciellt när man använder Euler-metoden av låg ordning .

Liknande metoder kan utvecklas för metoder av högre ordning, såsom 4:e ordningens Runge–Kutta-metoden. Dessutom kan en global feltolerans uppnås genom att skala det lokala felet till global omfattning.

Inbäddade feluppskattningar

Adaptiva stegstorleksmetoder som använder en så kallad "inbäddad" feluppskattning inkluderar metoderna Bogacki–Shampine , Runge–Kutta–Fehlberg , Cash–Karp och Dormand–Prince . Dessa metoder anses vara mer beräkningseffektiva, men har lägre noggrannhet i sina feluppskattningar.

För att illustrera idéerna med den inbäddade metoden, överväg följande schema som uppdaterar $y_{n}$ :

$y_{n+1}=y_{n}+h_{n}\psi (t_{n},y_{ n},h_{n})$

$t_{n+1}=t_{n}+h_{n}$

Nästa steg $h_{n}$ förutsägs från föregående information $h_{n}=g(t_{n}, y_{n},h_{n-1})$ .

För den inbäddade RK-metoden inkluderar beräkningen av $\psi$ en lägre ordnings RK-metod ${\tilde {\psi }}$ . Felet kan då enkelt skrivas som

${\textrm {err}}_{n}(h)={\tilde {y}}_{n+1}-y_{n+1}=h({\tilde {\psi }}(t_{n},y_ {n},h_{n})-\psi (t_{n},y_{n},h_{n}))$

${\textrm {err}}_{n}$ är det onormaliserade felet. För att normalisera det jämför vi det med en användardefinierad tolerans, som består av den absoluta toleransen och den relativa toleransen:

${\textrm {tol}}_{n}={\textrm {Atol}}+{\textrm {Rtol} }\cdot \max(|y_{n}|,|y_{n-1}|)$

$E_{n}={\textrm {norm}}( {\textrm {err}}_{n}/{\textrm {tol}}_{n})$

Sedan jämför vi det normaliserade felet $E_{n}$ mot 1 för att få det förutsagda $h_{n}$ :

$h_{n}=h_{n-1}(1/E_{n})^{1/(q+1 )}$

Parametern q är den ordning som motsvarar RK-metoden ${\displaystyle {\tilde {\psi }}} ,$ som har lägre ordning. Ovanstående prediktionsformel är rimlig i den meningen att den förstorar steget om det uppskattade lokala felet är mindre än toleransen och det krymper steget annars.

Beskrivningen ovan är en förenklad procedur som används i stegstorlekskontrollen för explicita RK-lösare. En mer detaljerad behandling finns i Hairers lärobok. ODE-lösaren i många programmeringsspråk använder denna procedur som standardstrategi för adaptiv stegstorlekskontroll, vilket lägger till andra tekniska parametrar för att göra systemet mer stabilt.

Se även

• Adaptiv kvadratur
• Adaptiv numerisk differentiering

Vidare läsning

•   William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, Numerical Recipes in C , Second Edition, CAMBRIDGE UNIVERSITY PRESS, 1992. ISBN 0-521-43108-5
•   Kendall E. Atkinson, Numerical Analysis , andra upplagan, John Wiley & Sons, 1989. ISBN 0-471-62489-6