Adaptive 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:

\[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.

Show code cell source Hide code cell source

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