Github Repository
Adaptive Schrittweitensteuerung
Contents
Voraussetzungen
Runge-Kutta-Verfahren zur numerischen Integration von gewöhnlichen Differentialgleichungen
Lerninhalte
Grundkenntnisse der adaptiven Schrittweitensteuerung
Adaptive Schrittweitensteuerung¶
ODE-Lösungsroutinen nutzen adaptive Schrittweitensteuerungen, um an sensiblen Stellen der Lösung möglichst nicht zu divergieren und um an weniger sensiblen Stellen den Rechenaufwand gering zu halten. Aufbauend auf den Kapiteln zu Runge-Kutta-Verfahren ist naheliegend, die Schrittweite zu halbieren und zu testen, ob Konvergenz erreicht wurde. Anders gesagt:
“Lohnt es sich noch, die Schrittweite zu verkleinern?”¶
Wir wissen, dass mit einer kleineren Schrittweite \(h\) ein genaueres Ergebnis für \(y_{i+1}\) zu erwarten ist. Beim expliziten Eulerverfahren hat der Gesamtfehler eine Größenordnung \(\mathcal{O}(h)\), also halbiert sich der Fehler bei Halbierung der Schrittweite. Das können wir nutzen, um zu prüfen, ob die bisherige Schrittweite klein genug war.
Aufgabe 1: Adaptive Schrittweitensteuerung durch Halbierung¶
Implementieren Sie eine adaptive Schrittweitensteuerung indem Sie die Schrittweite des expliziten Eulerverfahrens halbieren bis die Lösung hinreichend konvergiert.
function ydot = dgl(y,d)
ydot = -y/sqrt(d^2-y^2);
end
d = 1;
y0 = 0.999*d;
xmax = 10;
i = 1;
h0 = 0.5;
x = 0;
rtol = 1e-5;
hmin = 1e-5;
y(1) = y0;
while x < xmax
...
y_alt = ...
while ...
h = h/2;
y_neu = ...
err = ...
y_alt = y_neu;
end
x = x + h;
xspan(i) = x;
y(i) = y_neu;
end
plot(xspan,y)
“Ist die Ordnung meines Verfahrens ausreichend?”¶
Das Dormand-Prince Verfahren ist in ode45 implementiert. Es fragt anders als wir in Aufgabe 1: “Ist die Ordnung meines Verfahrens ausreichend?” Dafür werden zwei Lösungen berechnet und durch den Vergleich der Fehler geschätzt. Diese nennen sich dann eingebettete Verfahren.
Ein großer Vorteil ist, dass sich die neue Schrittweite schätzen lässt. Das hierfür notwendige Verfahren niedrigerer Ordnung lässt sich für Runge-Kutta-Verfahren immer durch eine Abwandlung der Gewichte finden. Dann lässt sich herleiten, dass für die neue Schrittweite gilt:
Dabei ist \(y_{m+1}\) das Ergebnis des Verfahrens höherer Ordnung, \(\hat{y}_{m+1}\) das des Verfahrens niedrigerer Ordnung, \(\text{tol}\) die Fehlertoleranz und \(p\) die höhere der beiden Ordnungen.
Ein Beispiel für ein eingebettetes Verfahren mit einer adaptiven Schrittweitensteuerung ist im ausklappbaren Fenster angegeben.
function [T,Y] = rk3_simple(fcn,tspan,y0,atol,rtol)
%-------------------------------------------------------------------------
%-- fcn: ode-function
%-- tspan: time interval
%-- y0: initial values
%-- atol,rtol: tolerances
%-------------------------------------------------------------------------
%-- Initialization -------------------------------------------------------
nfevals = 0; nsteps = 0; nrejected = 0; %-- statistics
uround=eps; fac1=0.2; fac2=6.0; %-- could be altered
%-- Method's coefficients
stage = 3; b = [1/6;1/6;2/3]; bd = [1/2;1/2;0];
a = [0;1;1/2]; alpha = [0,0,0;1,0,0;1/4,1/4,0];
% --Preparations for output-parameter [T,Y] ------------------------------
neq = length(y0); t0=tspan(1); tend=tspan(end); y(:,1) = y0; t=t0; T=t0; Y=y';
%-- INITIAL PREPARATIONS for stepsize --------
if (tend <= t0) error('tend<= t0 !!'); end
hmax=abs(tend-t0)/10; hinit=1.0e-6*abs(tend-t0); h=min(hinit,hmax);
%-- The main integration loop ---------------------------------------------
done = false; reject = false; K=zeros(neq,stage);
while ~done
if (t+0.1*h == t) || (abs(h) <= uround) %-- stepsize to small
fprintf('Error exit of rk3_simple at time t=%15g: step size to small h=%15g \n',t,h);
break
end
if (t+h) >= tend
h=tend-t;
else
h=min(h,0.5*(tend-t));
end
%---- integration step --------
if reject==false
K(:,1) = fcn(t,y); nfevals = nfevals+1;
end
for i=2:stage
sum_K = K*alpha(i,:)';
K(:,i) = fcn(t+h*a(i),y+h*sum_K); nfevals=nfevals+1;
end
sum_1 = h*K*b; sum_2 = h*K*bd; ynew = y + sum_1;
%---- error test -----------
SK = atol + rtol.*max(abs(y),abs(ynew));
err = sum( ((sum_1-sum_2)./SK).^2 ); err = max(realmin,sqrt(err/neq));
fac = 0.9 * err^(-1/3); fac=min(fac2,max(fac1,fac)); hnew=h*fac;
if (err <= 1.0) % --- STEP IS ACCEPTED -------
reject = false; nsteps = nsteps+1; told = t; t=t+h;
if (abs(tend-t) < uround) done=true; end %-- succesfull integration --
y = ynew; T=[T;t]; Y=[Y;ynew'];
else % --- STEP IS REJECTED -------
reject = true; nrejected = nrejected+1;
end
h = min(hnew,hmax);
%---------
end
%---------
fprintf('%g successful steps\n', nsteps);
fprintf('%g failed attempts\n', nrejected);
fprintf('%g function evaluations\n', nfevals);
end
Aufgabe 2: Adaptive Schrittweitensteuerung mit eingebetteten Verfahren¶
Gegeben seien für ein eingebettetes Verfahren die Runge-Kutta-Verfahren mit folgenden Butcher-Tableaus (die Herleitung aus den Ordnungsbedingungen bleibt Ihnen hier erspart). Für die Ordnung \(p=3\):
\(\begin{array} {c|ccc} 0\\ 1 & 1\\ \frac{1}{2} &\frac{1}{4} &\frac{1}{4} \\ \hline & \frac{1}{6} &\frac{1}{6} &\frac{4}{6} \end{array}.\)
Für die Ordnung \(\hat{p}=2\):
\(\begin{array} {c|ccc} 0\\ 1 & 1\\ \frac{1}{2} &\frac{1}{4} &\frac{1}{4} \\ \hline & \frac{1}{2} &\frac{1}{2} &0 \end{array}.\)
Butcher-Tabelaus eingebetteter Verfahren
Beide Tableaus lassen sich auch zusammenfassen:
\(\begin{array} {c|ccc} 0 \\ 1 & 1 \\ \frac{1}{2} & \frac{1}{4} & \frac{1}{4} \\ \hline p=3 & \frac{1}{6} & \frac{1}{6} & \frac{4}{6}\\ \hline \hat{p}=2 & \frac{1}{2} & \frac{1}{2} & 0 \end{array}\)
Leiten Sie aus den Butcher-Tableaus die Koeffizienten für \(K_1 ... K_3\), \(y_{i+1}\) und \(\hat{y}_{i+1}\) her.
Ändern Sie Ihren Code aus Integration von ODEs mit expliziten Verfahren zu einem eingebetteten Verfahren ab.
Vergleichen Sie das Ergebnis, insbesondere die Laufzeit, mit Ihrem Code aus Aufgabe 1.
%your code here