Anhang D: Numerisches Lösen von Differentialgleichungen

Die wenigsten Differentialgleichung-Systeme, die man in der Modellierung aufstellt, lassen sich analytisch lösen. Man löst die Gleichungen daher numerisch. Dies bedeutet, das man die Lösung Stück für Stück bestimmt. Um das Grundprinzip einfach darzustellen wird nur das Lösen einer einzelnen Differentialgleichung beschrieben. Das Lösen von Differentialgleichung-Systemen geht aber genauso. Zugrunde gelegt wird jeweils ein Anfangswertproblem:

\[y'=f(t,y(t)), \quad y(t_0)=y_0\]

D.1 Euler Verfahren

Zuerst wird an der Stelle \((t_0,y_0)\) die Steigung bestimmt. Es wird nun angenommen, dass sich die Steigung in einem kleinen Schritt \(\Delta t\) (Schrittweite) nicht ändert. Die Lösung wird also in dem Intervall \([t,t+\Delta t]\) durch eine Gerade approximiert. Man erhält den approximativen Wert \(y_1\) zur Zeit \(t_1\). Nun wird die Steigung an dieser Stelle bestimmt und so weiter:

\[y(t_0+\Delta t) = y(t) + f(t_0,y(t_0))\cdot \Delta t\]

Allgemein wird also zu einer beliebigen Zeit \(t\) der neue Wert zur Zeit \(t+\Delta t\) durch

\[y(t+\Delta t) = y(t) + f(t,y(t)) \cdot \Delta t\]

bestimmt. Ausgedrückt in Zeitschritten bedeutet dies:

\[y_{k+1} = y_k + f(t_k,y_k)\cdot \Delta t\qquad k=0,1,..\]

Dieses Verfahren ist recht ungenau.

Abbildung D.1: Euler-Verfahren zu \(\dot{y}= f(t,y)\), \(y(t_0)=y_0\).

D.1.1 Beispiel \(y'=y\), \(y(0)=1\)

Gegeben ist das Anfangswertproblem \(y'=y\), \(y(0)=1\), dessen explizite Lösung \(y(t)=e^t\) bekannt ist. Bestimmt man die Lösung zur Zeit T=10 mit dem Euler-Verfahren und einer Schrittweite von \(\Delta t=1\), so erhält man als approximative Lösung \(1024\). Der exakte Wert beträgt (auf 2 Nachkommastellen gerundet) 22026.47.

D.1.2 Implementierung in des Euler-Verfahrens in Matlab 

Eine Matlab-Programm zur Lösung des Anfangswertproblems \(\dot{y}= f(t,y)\), \(y(t_0)=y_0\) mit dem Euler-Verfahren ist in Listing D.1 angegeben.

% exponential growth
t(1)=0;     % start time
x(1)=1;         % initial value
dt=0.1;         % step size
tend=10;        % end time
r=1;            % parameter
kend=tend/dt; % number of steps
% loop until tend
for k=1:kend
   t(k+1)=k*dt;
   x(k+1)=x(k)+r*x(k)*dt;
end
plot(t,x)% plotting values

Listing D.1: Matlab-Code zur Lösung einer Differentialgleichung mit dem Euler-Verfahren.

Bemerkung: Die Implementierung geht in MatLab auch eleganter

Ein Verfahren, das bessere Ergebnisse liefert, ist das Runge-Kutta-Verfahren in Anhang D.3.

D.2 Runge-Kutta-Verfahren 2. Ordnung  (Verbessertes Euler-Verfahren, Heun-Verfahren)

Beim Runge-Kutta-Verfahren 2. Ordnung werden die Steigungen in \((t,y)\) und \((t+\Delta t,y+f(t,y)\cdot \Delta t)\) gemittelt, statt nur die Steigung in \((t,y)\) zu berücksichtigen. Dies bedeutet, dass man mit dem Funktionswert an der Stelle \(t+\Delta t\), der sich nach einem Euler-Schritt ergibt, die Steigung an \(t+\Delta t\) noch einmal neu berechnet. Das Mittel der beiden Steigungen wird dann als endgültige Steigung für den Zeitschritt angenommen:

\[y(t+\Delta t) = y(t) + \frac{f(t,y(t)) +f(t+\Delta t,\overbrace{y(t)+f(t,y(t))\cdot \Delta t}^{{\text{Funktionswert nach}} \atop {\text{einem Euler-Schritt}}}) }{2} \cdot \Delta t\]

Ausgedrückt in Zeitschritten

\[y_{k+1} = y_k + \frac{f(t_k,y_k) +f(t_{k+1},y_k+f(t_k,y_k) \cdot \Delta t)}{2}\cdot \Delta t\qquad k=0,1,...\]

Abbildung D.2: Verbessertes Euler-Verfahren zu \(\dot{y}= f(t,y)\), \(y(t_0)=y_0\).

D.3 Runge-Kutta-Verfahren 4. Ordnung 

Das bekannteste Verfahren approximiert die Lösung bis zur 4. Ordnung (Abbildung D.3). Hierbei wird an einer weiteren Stützstelle bei \(t+\frac{\Delta t}{2}\) die Steigung berechnet. Auch dies geschieht wieder mit dem verbesserten Verfahren. Die endgültige Steigung für den Zeitschritt wird aus 4 Einzelsteigungen gewichtet gemittelt, wobei die Steigungen in der Mitte bei \(t+\frac{\Delta t}{2}\) ein höheres Gewicht bekommen. Es wird wie folgt berechnet:

\[y(t+\Delta t) = y(t) + \frac{s_1+2s_2+2s_3+s_4}{6} \cdot \Delta t\]

\[\begin{aligned} s_1 &= f(t,y(t)) & s_2 &= f(t+{ \frac{\Delta t}{2}},y(t)+s_1 \cdot { \frac{\Delta t}{2}}) \\[1ex] s_3 &= f(t+{ \frac{\Delta t}{2}},y(t)+ s_2\cdot { \frac{\Delta t}{2}} ) & s_4 &= f(t+ \Delta t,y(t)+s_3\cdot \Delta t ) \end{aligned}\]

Ausgedrückt in Zeitschritten

\[y_{k+1} = y_k + \frac{s_1+2s_2+2s_3+s_4}{6} \cdot \Delta t\qquad k=0,1,...\]

\[\begin{aligned} s_1 &= f(t_k,y_k) & s_2 &= f(t_k+{ \frac{\Delta t}{2}},y_k+s_1 \cdot { \frac{\Delta t}{2}}) \\[0.5ex] s_3 &= f(t_k+{ \frac{\Delta t}{2}},y_k+ s_2\cdot { \frac{\Delta t}{2}} ) & s_4 &= f(t_k+ \Delta t,y_k+s_3\cdot \Delta t ) \end{aligned}\]

Dieses Verfahren ist in vielen Programmen implementiert. Häufig wird es einfach nur Runge-Kutta-Verfahren genannt. Noch bessere Verfahren sind implizite Verfahren, Verfahren mit einer automatischen Zeitschrittanpassung resp. Prediktor-Korrektor-Verfahren.

Abbildung D.3: Runge-Kutta-Verfahren 4. Ordnung zu \(\dot{y}= f(t,y)\), \(y(t_0)=y_0\).

Das Runge-Kutta-Verfahren 4. Ordnung ist in vielen Programmiersprachen implementiert. Eine Matlab-Programm zur Lösung des Anfangswertproblems \(\dot{y}= f(t,y)\), \(y(t_0)=y_0\) mit dem Runge-Kutta-Verfahren ist in Listing D.2 angegeben.

function myode
t0=0;       % start time
y0=1;       % initial value
tend=10;    % end time
r=1;        % parameter
% call of ode solver
[t,y]=ode45(@f,[t0 tend],y0,[],r);
plot(t,y) %plotting values

%ODE function
function dydt=f(t,y,r)
dydt=r*y;

Listing D.2: Matlab-Code zur Lösung einer Differntialgleichung mit dem Runge-Kutta-Verfahren.

D.4 Fehlerabschätzung 

Wie groß ist der Fehler, den man beim numerischen Lösen macht? Hierzu betrachtet man das Taylorpolynom von \(y\) (siehe A.12):

\[\begin{aligned} y(t+\Delta t) &=& \sum_{k=0}^n \frac{1}{k!}y^{(k)}(t)\cdot \Delta t^k +\frac{1}{(n+1)!}\cdot y^{(n+1)}(\xi) \cdot \Delta t^{n+1} \quad \xi \in [t,t+\Delta t]\\ &=& p_n(t) + R_n(t) \end{aligned}\]

Wird anstelle der exakten Lösung \(y(t)\) das Taylorpolynom n-ten Grades betrachtet, so ist der Abbruchfehler in jedem Schritt proportional zu \(\Delta t^{n+1}\). Im Falle des Euler-Verfahrens wird das Taylorpolynom 1. Grades betrachtet:

\[y(t+\Delta t) = y(t) + y'(t)\cdot \Delta t\]

Der Abbruchfehler ist also in jedem Schritt proportional zu \(\Delta t^2\). Je kleiner der Zeitschritt ist, desto kleiner ist der Fehler. Man kann den Zeitschritt aber nicht beliebig klein machen, da dann nicht nur die Rechenzeit stark anwächst. Bei sehr kleinem Zeitschritt kommt auch zum Tragen, dass die Zahldarstellung einer Programmiersprache begrenzt ist. In der Regel beträgt bei einfacher Genauigkeit(single precision) die Anzahl signifikanter Stellen 8, bei doppelter Genauigkeit (double precision) 16.