Entwurfsaspekte für die Realisierung des numerischen Verfahrens

Bezüglich der Auswahl und der Umsetzung des numerischen Verfahrens stellen sich eine Reihe von Fragen, die einer Antwort bedürfen, bevor man mit der Realisierung eines Programms beginnen kann. Dazu gehört eine genau Festlegung, welche Eingabedaten vorliegen und welche Ausgaben das Programm berechnen soll, aber auch, wie die Anwendungsschnittstelle der Software aufgebaut sein soll, welche Struktur der Programmcode besitzen soll und welche Datenstrukturen erforderlich sind.

Programmtechnische Realisierung einer gewöhnlichen Differentialgleichung

Eine der ersten zu beantwortenden Fragen ist, wie man ein gewöhnliches Differentialgleichungssystem in Programmform realisiert. Nach (1.1) ist ein nichtautonomes Anfangswertproblem einer gewöhnlichen Differentialgleichung gegeben durch

y⁰(t) =f(t,y(t)) , y(t0) =y0 .

2.1 Vom mathematischen Verfahren zum Programm 21 Das Differentialgleichungssystem wird dabei beschrieben durch die Funktionf : R×_Rⁿ → _Rⁿ, die ei-ne reelle Zahl sowie eiei-nenn-dimensionalen reellwertigen Vektor als Argumente besitzt und daraus einen n-dimensionalen reellwertigen Vektor als Ergebnis der Funktionsauswertung berechnet. Es liegt nahe, die mathematische Funktionfdurch eine Funktion im Sinne der verwendeten Programmiersprache zu realisie-ren, die mit entsprechenden Parametern und Rückgabewerten ausgestattet wird. Einige Dinge sind hierbei jedoch zu beachten.

Unterstützung mehrerer Differentialgleichungen. Soll es möglich sein, mehrere verschiedene Differen-tialgleichungen mit rechten Seitenf1,f2, . . . mit Hilfe des implementierten Lösungsverfahrens zu lösen, gibt es verschiedene Möglichkeiten, dies programmtechnisch umzusetzen.

1. Verwendung des gleichen Bezeichnersfür die Namen aller Funktionenf1,f2, . . . innerhalb des Programms.

Da es in einigen Programmiersprachen (u. a. C und FORTRAN) nicht erlaubt ist, daß innerhalb eines Programms mehrere Funktionen mit gleichem Bezeichner vorhanden sind,¹müssen die Funktionen für die verschiedenen Differentialgleichungen in separaten Quelldateien abgelegt werden. Diese Quellda-teien können einzeln in Objektcode übersetzt, jedoch nicht statisch zu einem gemeinsamen Binärpro-gramm gebunden werden. Dieser Ansatz wird in vielen klassischen ODE-Codes verwendet (siehe auch Abschnitt2.1.3).

Möchte man nur ein einziges Binärprogramm bereitstellen, dann ist deshalb zum Wechseln der Diffe-rentialgleichung ein erneutes Binden des Binärprogramms erforderlich. Alternativ könnte man für jede Differentialgleichung ein separates Binärprogramm erzeugen. Eine weitere Möglichkeit wäre, die rech-ten Seirech-ten der Differentialgleichungen als separate dynamische Bibliotheken zu implementieren und die jeweilige Differentialgleichung durch dynamisches Binden der entsprechenden Bibliothek auszuwählen.

Der Vorteil der Verwendung nur eines Bezeichners besteht darin, daß innerhalb des Programmteils, der das Lösungsverfahren realisiert, dieser Bezeichner direkt zum Aufruf der entsprechenden Programm-funktion verwendet werden kann, ohne daß neben den Kosten für den Funktionsaufruf (Sprung zur Adresse der Funktion, evtl. Sichern und Wiederherstellen von Registern, Rücksprung) weitere Kosten entstehen.

2. Verwendung verschiedener Bezeichnerfür die Funktionsnamen, z. B. durch Verwendung eines geeigneten Präfixes wie etwa dem Namen der Differentialgleichung. In diesem Fall ist es problemlos möglich, alle Differentialgleichungen innerhalb eines einzelnen Binärprogramms zur Verfügung zu stellen. Allerdings muß man jetzt dafür Sorge tragen, daß innerhalb des Lösungsalgorithmus die passende Funktion aufge-rufen wird.

Dazu benötigt man eine zusätzliche Variable, die die Information enthält, welches Differentialgleichungs-system gelöst werden soll. Mit Hilfe dieser Information könnte man z. B. jede Funktionsauswertung mit einer Fallunterscheidung verbinden, die allerdings mit unnötig hohen Kosten (mehrfaches Auslesen und Vergleichen der Variablen) verbunden ist. Hinzu kommt, daß, wenn eine zusätzliche Differential-gleichung implementiert werden soll, alle Programmstellen, an denen eine solche Funktionsauswertung aufgerufen wird, modifiziert werden müssen.

Um die mit der Fallunterscheidung verbundenen Probleme zu umgehen, kann man die zusätzliche Va-riable dazu nutzen, die Adresse der Programmfunktion zu speichern, welche die rechte Seite der zu lösenden Differentialgleichung realisiert. Dann ist es möglich, die ausgewählte Funktion über diesen Funktionszeiger indirekt aufzurufen. Im Vergleich zur Verwendung gleicher Bezeichner entstehen hier zwar ebenfalls zusätzliche Kosten, da ein Speicherzugriff und ein indirekter Sprung erforderlich sind, um den Inhalt der Variablen zu ermitteln und die vorgefundene Adresse anzuspringen, diese sind aber wesentlich geringer als bei Verwendung einer Fallunterscheidung.

Bereitstellung von Komponentenfunktionen. Bei der Realisierung der rechten Seite als Programmfunk-tion hat man die Möglichkeit, die ProgrammfunkProgrammfunk-tion entweder als VektorfunkProgrammfunk-tion f : R×Rⁿ → Rⁿ zu gestalten, so daß ein einzelner Funktionsaufruf alle Komponenten des Ergebnisvektors berechnet, oder sie

1In einigen anderen Programmiersprachen ist dies erlaubt, wenn die Funktionen anhand der Parameterliste unterscheidbar sind.

Dies ist hier allerdings nicht der Fall, da die Funktionenf1,f2, . . . die gleichen Argumente und Rückgabewerte besitzen.

als skalare Funktion f :N×_R×_Rⁿ →_Rzu implementieren, die jeweils nur eine einzelne Vektorkompo-nente berechnet.

Wenn die Berechnungsvorschrift für die rechte Seite des Differentialgleichungssystems dies zuläßt, bie-tet die Bereitstellung von skalaren Komponentenfunktionen eine Reihe von Vorteilen gegenüber der Imple-mentierung als Vektorfunktion. So ermöglicht die Aufsplittung in Komponentenfunktionen u. a. eine da-tenparallele Auswertung der rechten Seite, indem auf jedem der verwendeten Prozessoren jeweils eine Teil-menge der Komponentenfunktionen berechnet wird. Beim Entwurf aller Implementierungen – sequentiell oder parallel – erlaubt die Verwendung von Komponentenfunktionen die Durchführung von Programm-transformationen, die andernfalls nicht möglich wären. Die Verwendung von Komponentenfunktionen ist daher die Voraussetzung für einen Großteil der Programmtransformationen zur Optimierung des Speicher-zugriffsverhaltens vieler sequentieller und paralleler Implementierungen, die in dieser Arbeit betrachtet werden.

Programmtechnische Realisierung des Integrationsalgorithmus

Die Realisierung des Integrationsalgorithmus setzt eine Festlegung der gewünschten Ausgaben und der dafür erforderlichen Eingangsinformationen voraus. Darüber hinaus ist zu überlegen, welches konkrete numerische Verfahren geeignet ist, um die an die Qualität der Lösung gestellten Anforderungen zu er-füllen, und welches gleichzeitig hinreichend effizient arbeitet, um die gesuchte Lösung innerhalb einer angemessenen Antwortzeit zu liefern.

Nach (1.1) werden Probleme aus der betrachteten Klasse nichtautonomer Anfangswertprobleme be-schrieben durch

t∈[t0,te]: y⁰(t) =f(t,y(t)) , y(t0) =y0 .

Gegeben ist also ein Intervall[t₀,te], eine Funktionf, deren programmtechnische Realisierung bereits ange-sprochen wurde, sowie der Anfangswerty0. Gesucht wird eine Approximation der Lösungsfunktiony(t), deren Ableitung im Intervall[t₀,te]durch die Funktionfbeschrieben wird. Es kann auch vorkommen, daß der Endpunkt des Intervalls, te, zu Beginn der Integration nicht bekannt ist, da nach einer Stelle gesucht wird, an der ein bestimmtes Ereignis eintritt, wie z. B. das Erfülltsein einer bestimmten Bedingung fürtund y(t).

Die Berechnung der Approximation der Lösungsfunktion erfolgt in der Regel in Form einer Gitterfunk-tion

η_κ =η(tκ)≈y(tκ) , κ=0, . . . ,N , die Näherungswerte für Funktionswerte vony(t)an mehreren Zwischenstellen

t₀<t₁<· · ·<t_N =te (2.1) des Intervalls[t0,te]liefert (vgl. (1.4) und (1.5)). Anzahl und Position der Zwischenstellen werden von der Schrittweitenkontrolle bestimmt. Die Ausgabe der berechneten Lösung kann je nach Anwendungsfall auf sehr unterschiedliche Weise erfolgen. Oft genügt die Ausgabe der Approximationswerte an den Gitterstel-lent0, . . . ,te. In bestimmten Fällen kann sogar nur der einzelne Approximationswert füry(te)am Rand des Intervalls gesucht sein. Einige Anwendungen können jedoch auch die Ausgabe von Approximationswerten an wesentlich dichter liegenden Zwischenstellen erfordern, als sie von der Schrittweitenkontrolle generiert werden, z. B. für eine grafische Darstellung der Lösungsfunktion.

Nahezu jede Implementierung eines Lösungsverfahrens mit Schrittweitenkontrolle gibt dem Anwender die Möglichkeit, durch Vorgabe eines oder mehrerer Parameter, die oft mit TOL bzw. ATOL und RTOL bezeichnet werden, Einfluß darauf zu nehmen, wie groß der geschätzte Betrag des lokalen Fehlers sein darf, damit ein Zeitschritt von der Schrittweitenkontrolle akzeptiert wird (vgl. Abschnitt1.1.1). Dadurch ist es für den Anwender möglich, den für seinen speziellen Anwendungsfall besten Kompromiß zwischen Genauigkeit und Laufzeit auszuwählen.

Neben der Ausgabe der Approximation der Lösungsfunktion können für den Anwender eine Reihe zusätzlicher Informationen über den Integrationsvorgang von Interesse sein, etwa über die Größe des loka-len und globaloka-len Fehlers und über die durchgeführten Zeitschritte. Kann die Integration nicht erfolgreich abgeschlossen werden, weil z. B. aus numerischen Gründen die Integration nicht fortgesetzt werden kann oder weil eine Ausnahmebehandlung erforderlich ist, sollte der Anwender detaillierte Informationen über die Ursache erhalten.

2.1 Vom mathematischen Verfahren zum Programm 23 Wahl der Verfahrenskoeffizienten. Der Kern jeder Implementierung eines eingebetteten Runge-Kutta-Verfahrens besteht in der Umsetzung der mathematischen Berechnungsvorschrift des Runge-Kutta-Verfahrens (vgl. (1.9), (1.16), (1.17)):

vl =_f t_κ+c_lh_κ,η_κ+h_κ

l−1

∑

i=1

a_livi

!

für l=_{1, . . . ,}s , (2.2a)

η_κ+1=η_κ+hκ

∑

s l=1

blvl , (2.2b)

ˆ

η_κ+1=ηˆ_κ+h_κ

∑

s l=1

bˆ_lv_l . (2.2c)

Die Koeffizientena_li,b_l, ˆb_l undc_l definieren das konkrete Verfahren und bestimmen dessen numerische Eigenschaften. Sie sollten deshalb sorgfältig anhand der Eigenschaften der zu lösenden Differentialglei-chungssysteme und der Anforderungen an die zu berechnende Approximation der Lösungsfunktion aus-gewählt werden. Verschiedene eingebettete Verfahren werden z. B. inHairer u. a.(2000) undButcher(2003) bezüglich ihrer numerischen Eigenschaften diskutiert und gegenübergestellt.

Die Entscheidung, welche Verfahrenskoeffizienten unterstützt werden sollen, hat weiterhin einen ent-scheidenden Einfluß auf die Möglichkeiten zur programmtechnischen Realisierung des Verfahrens in Form einer Implementierung. Soll die Implementierung beliebige eingebettete Verfahren unterstützen, müssen die Verfahrenskoeffizienten in dynamischen Datenstrukturen (z. B. ein- bzw. zweidimensionale Felder) ab-gelegt werden, deren Größe und Inhalt erst nach Vorgabe durch den Benutzer zur Laufzeit bekannt ist. Die Berechnungsvorschrift (2.2) muß in diesem Fall mit Hilfe von Schleifen realisiert werden, deren Grenzen von der Stufenzahlsund der Größe des Differentialgleichungssystemsnabhängen und ebenfalls erst zur Laufzeit bekannt sind.

Zusätzliche Möglichkeiten zur Optimierung des Programmcodes kann man durch eine Spezialisierung erreichen, indem nur bestimmte Verfahrensklassen oder sogar nur ein einziges Verfahren zugelassen wer-den. Legt man z. B. allein die Anzahl der Stufen fest, sind die Anzahl der Verfahrenskoeffizienten und die Grenzen der Schleifen, die über die Stufen iterieren, bereits zur Übersetzungszeit bestimmt und können für Optimierungen ausgenutzt werden. So kann man etwa die Verfahrenskoeffizienten in statischen Da-tenstrukturen (z. B. statische Felder oder skalare Variablen) vorhalten und durch Aufrollen von Schleifen zusätzliche Codetransformationen ermöglichen.

Einige eingebettete Runge-Kutta-Verfahren besitzen die sogenannte FSAL-Eigenschaft (englisch: first same as last). Das bedeutet, daß die Koeffizientenb₁, . . . ,bs−1identisch mit den Koeffizientena_s,1, . . . ,as,s−1

sind und bs = 0 gilt. Infolgedessen stimmen auch der Vektor vs aus dem Zeitschritt κ und der Vektor v₁aus dem nachfolgenden Zeitschrittκ+1 überein, so daß eine Auswertung der Funktionfeingespart werden kann. Bei einer Spezialisierung der Implementierung auf Verfahren mit der FSAL-Eigenschaft ist demzufolge ein signifikanter Laufzeitgewinn im Vergleich zur Ausführung derselben Verfahren mittels einer allgemeineren Implementierung zu erwarten.

Legt man sich beim Entwurf der Implementierung auf nur ein spezielles Verfahren fest, ergeben sich die umfangreichsten Möglichkeiten zur Optimierung. Da die Beträge der Verfahrenskoeffizienten bekannt sind, können diese nun direkt als Konstanten in den Quellcode eingefügt werden. Haben einige Koeffizi-enten den Betrag 0 oder 1, können Rechenoperationen eingespart werden, indem Additionen mit 0 und Multiplikationen mit 0 oder 1 vermieden werden. Dies kann dadurch erfolgen, daß man diese Operatio-nen erst gar nicht in den Programmcode aufnimmt oder daß ein optimierender Compiler diese Spezialfälle erkennt und effizienten Maschinencode generiert.

Vermeidung von Rundungsfehlern. Bei der Implementierung jedes numerischen Algorithmus sollte man versuchen, Rundungsfehler so weit als möglich zu minimieren. Die Verfahrensvorschrift für eingebettete Runge-Kutta-Verfahren (2.2) enthält neben den Funktionsauswertungen der rechten Seite des Differential-gleichungssystems ausschließlich Additions- und Multiplikationsoperationen, die zur Entstehung des Run-dungsfehlers beitragen. Besonders kritisch ist die Berechnung vonη_κ+1nach (2.2b) und ˆη_κ+1nach (2.2c), da die Komponenten der Vektoren hκ∑^sl=1b_lv_l bzw.hκ∑^sl=1bˆ_lv_l in der Regel sehr klein gegenüber den

Komponentenη_κ+1bzw. ˆη_κ+1sind. Die schlimmsten Auswirkungen der Akkumulation von Rundungsfeh-lern können jedoch vermieden werden, indem man den bei dieser Addition auftretenden Fehler abschätzt und ihn im nächsten Zeitschritt zum kleineren Summanden hinzuaddiert (Butcher 2003). Diese Technik bezeichnet man als kompensierte Summation (englisch: compensated summation). Higham (1993) gibt einen ausführlichen Überblick über verschiedene Ansätze für eine möglichst genaue Summation von Gleitkom-mazahlen.

Stetige Erweiterungen. Einige Anwendungen benötigen Näherungslösungen für die Werte der Lösungs-funktion außerhalb des von der Schrittweitenkontrolle erzeugten Punktgitters. Beispielsweise könnte für die Erzeugung einer grafischen Darstellung der Lösungsfunktion ein festes Punktgitter vorgegeben sein, dessen Gitterpunkte sehr dicht liegen können. Um die Lösung an allen geforderten Zwischenpunkten zu berechnen, könnte man die Schrittweite entsprechend reduzieren. Da hierdurch gleichzeitig die Anzahl der Zeitschritte und damit verbunden die Anzahl der Funktionsauswertungen ansteigt, ist dieses Vorgehen al-lerdings sehr ineffizient. Weiterhin führt die Erhöhung der Anzahl der Zeitschritte zu einer Zunahme von Rundungsfehlern.

Es wurden deshalb Runge-Kutta-Verfahren konstruiert, die für jeden ZeitschrittκApproximationen für alle Wertey(t_κ+θh_κ)im Bereich 0≤ θ≤ 1 liefern, wobei dafür keine oder nur wenige zusätzliche Funk-tionsauswertungen erforderlich sind. Man bezeichnet derartige Verfahren alsstetige Runge-Kutta-Verfahren.

Zum Teil wird in diesem Zusammenhang auch der Begriffstetige Erweiterungbenutzt, da eine Reihe steti-ger Verfahren aus klassischen Verfahren durch das Hinzufügen zusätzlicher Stufen konstruiert wurden. Die Grundlage für die Berechnung der Zwischenwerte bilden spezielle Interpolationsansätze. Eine ausführli-che Darstellung sowie weiterführende Literaturhinweise geben z. B.Hairer u. a.(2000) sowieEnright u. a.

(1995).

Anfangsschrittweite. Eine Möglichkeit zur Festlegung der Anfangsschrittweite besteht darin, vom An-wender eine Vorgabe der Anfangsschrittweite zu fordern. Dies setzt voraus, daß der AnAn-wender aus seiner Erfahrung oder aus mathematischen Überlegungen heraus eine Vorstellung von der zu wählenden Größe der Anfangsschrittweite besitzt. Um eine ungünstige Wahl der Anfangsschrittweite durch den Anwender zu vermeiden, ist eine bessere Vorgehensweise, die Anfangsschrittweite durch einen geeigneten Algorith-mus festzulegen bzw. dem Anwender einen entsprechenden Vorschlag zu unterbreiten. Ein möglicher Al-gorithmus zur Bestimmung der Anfangsschrittweite wurde vonGladwell u. a.(1987) vorgeschlagen (siehe auchHairer u. a. 2000).

Behandlung des Intervallendes. Zur Erzeugung des Punktgitters (2.1) wird die Schrittweite in jedem Zeitschritt unter Beachtung der Genauigkeitsforderung möglichst groß gewählt, um die Anzahl der Zeit-schritte und damit den Rechenaufwand zu minimieren. Wurden die Gitterpunkte t₀, . . . ,tN−1 auf diese Weise erzeugt, ergibt sich die Schrittweite für den letzten Zeitschritt zur Bestimmung des Ausgabepunktes am Intervallende te = tN als hN−1 = tN−tN−1. Die Schrittweite wird in diesem Zeitschritt also kleiner gewählt, als es aufgrund der Genauigkeitsforderung notwendig wäre, um den Ausgabepunkt am Intervall-ende ermitteln zu können. Dies kann zu einem erhöhten Rechenaufwand führen, insbesondere wenn die Integration auf dem sich ante anschließenden Intervall fortgesetzt wird und dieses Vorgehen systematisch für weitere Intervalle fortgeführt wird.

Um die durch die Verkleinerung der letzten Schrittweite hervorgerufenen Effekte abzumildern, kann man versuchen, das Intervallende „vorauszuahnen“ und die Schrittweite bereits in früheren Zeitschritten zu reduzieren. Ein anderer Ansatz besteht in der Verwendung eines stetigen Runge-Kutta-Verfahrens. In diesem Fall wird die Schrittweite nicht eingeschränkt, und die Integration wird fortgesetzt, bis das Inter-vallendete erstmalig überschritten wird. Die Approximation des Funktionswerts y(te)wird dann durch Interpolation ermittelt (sieheEnright u. a. 1995).

Abbruchkriterien und besondere Ereignisse. In bestimmten Situationen ist es sinnvoll, die Integration bereits vor dem Erreichen des Intervallendes zu beenden bzw. neu zu starten. Eine solche Situation liegt z. B. vor, wenn erkannt wird, daß der Rechenaufwand zur Berechnung der gesuchten Lösung zu groß wäre, um die Berechnung innerhalb einer akzeptablen Zeitspanne zu beenden, oder wenn das Anwachsen von Rundungsfehlern zu befürchten ist. Mögliche Abbruchkriterien können sein:

2.1 Vom mathematischen Verfahren zum Programm 25 1. Eine zu kleine Schrittweite, da diese zu einer großen Anzahl von Zeitschritten und damit zu einem hohen

Berechnungsaufwand und der Akkumulation von Rundungsfehlern führt.

2. Die Berechnung einer vorgegebenen maximalen Anzahl von Zeitschritten.

3. Die Detektion von Steifheit (siehe Hairer u. Wanner 2002), da explizite Verfahren für die Integration steifer Differentialgleichungssysteme ungeeignet sind.

4. Das Auftreten einer Singularität (sieheEnright u. a. 1995).

5. Das Erfülltsein einer vom Nutzer vorgegebenen Bedingung (sieheEnright u. a. 1995;Hairer u. a. 2000;

Shampine u. Thompson 2000).

6. Das Auftreten von Programmausnahmen, verursacht z. B. durch das Fehlschlagen von Betriebssystem-aufrufen (etwa Anforderung von Speicher oder Ein-/Ausgabe-Funktionen).

Andere Situationen können eine spezielle Behandlung, z. B. einen Neustart oder einen Strategiewechsel der Schrittweitenkontrolle, erfordern. Dies ist beispielsweise der Fall, wenn innerhalb des Integrations-intervalls Unstetigkeiten auftreten oder sehr viele aufeinanderfolgende Schritte verworfen werden (siehe auchEnright u. a. 1995).

Statistische Informationen. Zusätzlich zur berechneten Näherungslösung des gestellten Anfangswert-problems können für den Anwender statistische Informationen über die Eigenschaften der Lösung sowie den Integrationsvorgang von Interesse sein. Beispielsweise könnte sich der Anwender für Aussagen über die Genauigkeit der berechneten Näherungslösung interessieren, also etwa den lokalen und globalen Feh-ler. Eine Approximation des lokalen Fehlers wird in jedem Zeitschritt berechnet, da sie für die Schrittwei-tenkontrolle benötigt wird. Zur Bestimmung des globalen Fehlers wurden verschiedene Techniken vorge-schlagen, sieheEnright u. a.(1995) für einen Überblick.

Um die Effizienz des Integrationsvorgangs einschätzen zu können, sind Informationen über die ausge-führten Zeitschritte, d. h. die Gesamtanzahl sowie der Anteil der akzeptierten bzw. verworfenen Schritte und die dazugehörigen Schrittweiten, hilfreich.

Für eine detaillierte Untersuchung des Laufzeit- und Lokalitätsverhaltens ist außerdem die Ausgabe diesbezüglicher Meßdaten erforderlich. Dabei kann es sich beispielsweise um Laufzeiten für verschiedene Programmteile oder Informationen über das Auftreten bestimmter Hardware-Ereignisse, wie z. B. Speicher-zugriffe auf die verschiedenen Ebenen der Speicherhierarchie oder die Ausführung bestimmter Instrukti-onstypen, handeln.

Schnittstelle zum Anwendungsprogramm. Bei Verwendung einer nicht-objektorientierten, imperativen Programmiersprache wie C oder FORTRAN wird man ein Lösungsverfahren in der Regel in Form eines oder mehrerer Unterprogramme realisieren, um den Code des Lösungsverfahrens vom Anwendungspro-gramm zu trennen und so den Code auch für andere AnwendungsproAnwendungspro-gramme nutzbar zu machen.

Für die Kommunikation zwischen dem Unterprogramm des Integrators und dem aufrufenden Anwen-dungsprogramm gibt es verschiedene Möglichkeiten:

I Parameterübergabe:Daten werden beim Aufruf bzw. der Rückkehr des Integrators als Parameter über-geben.

I Globale Variablen:Daten werden vom Anwendungsprogramm vor Aufruf des Integrators in globalen Variablen abgelegt, die später vom Integrator ausgelesen werden. Der Integrator legt Rückgabedaten in globalen Variablen ab, die nach Beendigung dem Anwendungsprogramm zur Verfügung stehen. Im einfachsten Fall kann dem Anwendungsprogramm der direkte Zugriff auf diese globalen Variablen ge-stattet werden. Sicherer ist es aber, dem Anwendungsprogramm das Setzen und Auslesen der Variablen nur über spezielle Funktionen zu gestatten, da so eine Überprüfung der Zugriffe möglich ist.

Soll während der Integration Programmcode des Anwendungsprogramms ausgeführt werden, kann eine der folgende Techniken eingesetzt werden:

I Call-Back-Funktionen:Das Anwendungsprogramm stellt Unterprogramme bereit, die vom Integrator aufgerufen werden. Daten können über Parameter bzw. Rückgabewerte ausgetauscht werden, aber auch der Datenaustausch über globale Variablen ist möglich.

I Rücksprung zum Anwendungsprogramm:Das Unterprogramm des Integrators wird zwischenzeitlich beendet, und das Anwendungsprogramm wird fortgesetzt. Der Status des Integrators wird entweder in globalen Variablen festgehalten oder als Rückgabe an das Anwendungsprogramm übergeben. Zur Fortsetzung der Integration ruft das Anwendungsprogramm den Integrator erneut auf.

I Nebenläufige Ausführung:Integrator und Anwendungsprogramm werden als separate Kontrollflüsse realisiert, die auf geeignete Weise miteinander kommunizieren.

Die Ausführung von Programmcode des Anwendungsprogramms während der Integration wird in der Praxis häufig eingesetzt, um nacheinander anfallende Zwischenergebnisse, wie z. B. die Approxima-tionswerte η_κ, anwendungsspezifisch zu verarbeiten. Existierende Codes lassen sich meist diesbezüglich klassifizieren als:

I Schrittorientierte Codes:Der Integrator berechnet jeweils einen Zeitschritt und kehrt anschließend zum Anwendungsprogramm zurück. Zur Berechnung weiterer Zeitschritte muß der Integrator wiederholt aufgerufen werden.

I Intervallorientierte Codes: Es wird ein Intervall vorgegeben. Der Integrator berechnet so viele Zeit-schritte, wie zum Durchlaufen des Intervalls erforderlich sind. Zur Verarbeitung von Zwischenergebnis-sen durch das Anwendungsprogramm werden Call-Back-Funktionen verwendet.

Shampine u. a.(1976) geben einen Überblick über die wichtigsten damals verbreiteten Codes für die Lösung nichtsteifer Anfangswertprobleme. Diese Autoren haben die Erfahrung gemacht, daß Anwender oft intervallorientierte Codes bevorzugen.

Nebenläufige Codes sind bisher kaum verbreitet, obwohl sie ein signifikantes Potential zur Laufzeit-verbesserung bieten. Beispielsweise führt ein typisches Anwendungsprogramm eine Reihe von Ausgabe-operationen aus, z. B. zur grafischen Darstellung der Lösungsfunktion, zur textuellen Ausgabe auf dem Bildschirm oder zur Speicherung der Approximationswerte in einer Datei für die spätere Verarbeitung.

Die Wartezeiten, die beim Zugriff auf die Ausgabegeräte entstehen, können bei einer nebenläufigen Reali-sierung verdeckt werden. Ein denkbares Vorgehen wäre eine Implementierung auf Basis zweier Threads, wobei der Integrationsthread die schrittweise berechneten Approximationswerte über einen Produzenten-Konsumenten-Mechanismus an den Anwendungsthread übermittelt. Eine derartige Realisierung wäre ins-besondere für moderne Prozessoren, die ein hardwarebasiertes (simultanes) Multithreading oder sogar

Die Wartezeiten, die beim Zugriff auf die Ausgabegeräte entstehen, können bei einer nebenläufigen Reali-sierung verdeckt werden. Ein denkbares Vorgehen wäre eine Implementierung auf Basis zweier Threads, wobei der Integrationsthread die schrittweise berechneten Approximationswerte über einen Produzenten-Konsumenten-Mechanismus an den Anwendungsthread übermittelt. Eine derartige Realisierung wäre ins-besondere für moderne Prozessoren, die ein hardwarebasiertes (simultanes) Multithreading oder sogar

Im Dokument Effiziente Implementierung eingebetteter Runge-Kutta-Verfahren durch Ausnutzung der Speicherzugriffslokalität (Seite 36-42)