Numerische Lösung von Delaydifferentialgleichungen

In diesem Abschnitt wird das Verfahren beschrieben, mit dem in dieser Arbeit die Delaydifferentialgleichungen numerisch gelöst wurden. Es basiert auf einem Verfahren, das von BAKER et al. 1994 in Manchester in [4] beschrieben wurde.

Eine gewöhnliche Differentialgleichung hat die Form

$$\dot{y}(t) = F(t, y(t), \vec{\mu}), \quad \text{mit einem Parametervektor } \vec{\mu}.$$ (49)

Dagegen hat die entsprechende Delaydifferentialgleichung die Form

$$\dot{y}(t) = F(t, y(t), y(t-\tau) ,\vec{\mu}),$$ (50)

mit einem (in dieser Arbeit) konstanten Delay $\tau > 0$. Wenn die Ableitung $\dot{y}(t)$ zudem noch von der Ableitung zu einem früheren Zeitpunkt abhängt,

$$\dot{y}(t) = F(t, y(t), y(t-\tau), \dot{y}(t-\tau) ,\vec{\mu}),$$ (51)

spricht man von einer neutralen Delaydifferentialgleichung ([32],[12],[4]). Da im Keller-Miksis-Modell im Kopplungsterm auch die Ableitung des Blasenradius der ankoppelnden Blase auftaucht, liegt hier also eine solche neutrale Delaydifferentialgleichung vor. Der Zustand eines solchen Systems wird nur durch die abhängigen Variablen nicht ausreichend beschrieben. Vielmehr benötigt man noch die Funktionswerte eines Intervalls, das bis zum Zeitpunkt $t-\tau$ zurückreicht. Somit hat man also eine unendlichdimensionale Differentialgleichung vorliegen, für die eine aufwendige Behandlung notwendig ist. Im Folgenden wird das Lösungsverfahren für eine einfache Delaydifferentialgleichung gezeigt, wie sie BAKER et al. beschrieben haben. Dieses Verfahren lässt sich leicht auf neutrale Delaydifferentialgleichungen anwenden.

Die Lösung der Gleichung erfolgt auf kleinen Intervallen [ $\sigma_k,\sigma_{k+1}$], auf denen die Gleichung in eine gewöhnliche Differentialgleichung transformiert wird. Dies ist der wesentliche Schritt zur Lösung einer Delaydifferentialgleichung nach der Methode, die Baker die ,,Method Of Steps`` nennt.

Sei $\{\sigma_i\}$ die Menge aller Punkte, die durch $t_0 +i \tau$ beschrieben werden. Nun setzt man auf jedem so entstandenem Teilintervall [ $\sigma_k,\sigma_{k+1}$]:

$$\dot{y}_1(t) = F(t, y_1(t), \Psi(t-\tau),\vec{\mu}) (\sigma_0 \leq t \leq \sigma_1)$$

$$\dot{y}_2(t) = F(t, y_2(t), y_1(t-\tau),\vec{\mu}) (\sigma_1 \leq t \leq \sigma_2)$$

$$\vdots \notag$$ (52)

$$\dot{y}_m(t) = F(t, y_m(t), y_{m-1}(t-\tau),\vec{\mu}) (\sigma_{m-1} \leq t \leq \sigma_m),$$ (53)

wobei $\Psi(t)$ eine Funktion ist, die $y(t)$ für $t_0-\tau \leq t \leq t_0$ festlegt.
Damit ist $y(t) = y_l(t) \text{ f\uml {u}r } t \in [\sigma_{l-1};\sigma_l]$ und $y(t) = \Psi(t) \text{ f\uml {u}r } t_0-\tau \leq t \leq t_0$. Wenn nun $y_{m-1}(t)$ bereits bekannt ist, ist

$$\dot{y}_m(t) = F(t, y(t),y_{m-1}(t-\tau), \vec{\mu}) \quad (\sigma_{m-1} \leq \sigma_m)$$ (54)

eine gewöhnliche Differentialgleichung, die mit einer beliebigen um eine Interpolationsroutine zur Beschaffung der verzögerten Funktionswerte erweiterten Routine (s. u. ) numerisch gelöst werden kann.

Wie BAKER et al. beschrieben haben, funktioniert dieses Verfahren auch für ein variables Delay solange es nur nicht verschwindet. In einem solchen Fall wird lediglich die Beschreibung der $\sigma_i$ etwas schwieriger.

Man muß nun dafür sorgen, dass die Lösung dicht ist. Benötigt man einen verzögerten Funktionswert, so darf man nicht davon ausgehen, diesen bereits berechnet zu haben, schließlich erfolgt eine numerische Lösung immer nur in diskreten Schritten endlicher Länge. Somit wird eine Interpolationsroutine bei der Beschaffung der verzögerten Werte eingesetzt, um diese Anforderung zu erfüllen. Dabei sollte die Ordnung der Interpolationsroutine höher sein als die der Integrationsroutine ([27]), damit deren Ordnung sichergestellt wird. Diese Interpolationsroutine bewältigt zudem die unendlich vielen Dimensionen, die durch die Funktionswerte des zurückliegenden Intervalls in die Delaydifferentialgleichung kommen.

Ein weiter Punkt, der besonderer Beachtung bedarf, ist die Fortpflanzung von Diskontinuitäten. Tritt ( durch entsprechende Wahl der Initialisierungsfunktion $\Psi(t)$) an der Stelle $t_0=0$ eine solche Unstetigkeit auf, so setzt sie sich über die Ableitung fort und führt zum Zeitpunkt $t_1=\tau$ zu einem Knick in der Lösung. Zum Zeitpunkt $t_2=2 \tau$ findet sich die Unstetigkeit in der 2. Ableitung wieder, usw. Die Lösung wird jedoch während der Fortpflanzung immer glatter, so dass sich der Einfluss dieser Unstetigkeit immer mehr verliert. Dieses Verhalten gilt selbstverständlich nur für ,,normale`` Delaydifferentialgleichungen. Bei neutralen Gleichungen bleibt die nicht differenzierbare Stelle stets erhalten und ,,wandert`` nicht in höhere Ableitungen, da sie ja immer wieder zu einem Sprung in der Ableitung führt. In den Rechnungen der hier vorliegenden Arbeit wurde die Initialisierungsfunktion meist einem bereits (meist ohne Delay) eingeschwungenen System entnommen oder aber konstant auf die Startwerte der betroffenen Variablen gesetzt, um dieses Problem zu umgehen einzuschränken. Lesenswert sind in dieser Richtung noch die Ergebnisse aus [24] von LOSSON et al. , wo gezeigt wird, dass eine um ein Delay erweiterte Form der kontinuierlichen logistischen Abbildung bei bestimmten Parametern zwei koexistierende Attraktoren besitzt, deren Erreichen empfindlich von der Initialisierungsfunktion abhängt.

Es liegt auf der Hand, dass die Schrittweite des Integrators kleiner als das Delay sein muß, damit die Delaydifferentialgleichung nicht zu einer impliziten Gleichung wird, für die eine andere Form der Behandlung notwendig wäre, oder aber für die Interpolationsroutine nicht mehr ausreichend Stützstellen zur Verfügung stehen, um die nötige Ordnung zu erreichen.