{grid}
```{grid-item-card}
:class-header: bg-light
Voraussetzungen
^^^
- Runge-Kutta-Verfahren zur numerischen Integration von gewöhnlichen Differentialgleichungen
```
```{grid-item-card}
:class-header: bg-light
Lerninhalte
^^^
- Grundkenntnisse der adaptiven Schrittweitensteuerung
```
(rk_adaptive_steps)=
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:
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.
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)
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:
$$h_\text{neu} = h \left\lVert\frac{y_{m+1}-\hat{y}_{m+1}}{\text{tol}}\right\rVert^{-\frac{1}{p}}.$$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
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}.$
{admonition}
:class: dropdown
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}$
integration_rk_exp
zu einem eingebetteten Verfahren ab.%your code here