Vorlesung Applikationssoftware in der Physik 2, SS 1998
Auf dieser Page möchten wir dem Benutzer zeigen, wie man mit Hilfe von Bibliotheksroutinen numerische Probleme behandeln kann. Im Vordergrund steht dabei die prinzipielle Vorgehensweise, nicht das mathematische Problem an sich. Dieses ist vielmehr Gegenstand der Vorlesungen über numerische Methoden von Prof. Sormann. Wir haben daher bewußt einfache Problemstellungen gewählt, von denen wir die Lösungen kennen. Verwendet werden die folgenden Bibliotheken:
Als Programmiersprache wählen wir C.Im folgenden soll kurz das Vorgehen behandelt werden, wie man bei der Verwendung der verschiedenen Bibliotheken von einem gegeben Problem zu einem lauffähigen C-Programm gelangt. Die angegebenen Befehle gehen davon aus, daß man auf dem Physikserver der TU-Graz fubphpc.tu-graz.ac.at eingeloggt ist. Auf anderen Rechnern kann es sein, daß für das Compilieren des Programmes einige Verzeichnispfade erst richtig gesetzt werden müssen (der Compiler muß die entsprechenden header-files finden).
Gegeben sei ein numerisches Problem, z.B. eine transzendente Gleichung, deren Lösungen gesucht sind. Die grobe Vorgehensweise ist immer dieselbe:
Auf dem Server des Institutes für theoretische Physik gibt es eine komplette Online-Hilfe für die NAG-Bibliothek. Man loggt sich dort mit dem Usernamen 'naghelp' ein und das Hilfeprogramm wird automatisch gestartet:
> xhost + > telnet faepal02Die Verwendung der Dyna-Text gesteuerten Hilfe ist selbsterklärend, sie ähnelt in der Bedienung einem Web-Browser und sollte daher den meisten Benutzern keine Probleme bereiten. Man klickt mit der Maus einfach die Rubriken an, die man lesen will. Es ist zu empfehlen, sich vor dem ersten Arbeiten mit der NAG das einführende Kapitel ein wenig durchzulesen.
Zunächst muß man sich entscheiden, ob man die C oder die FORTRAN Bibliothek durchforsten möchte. Dann gelangt man in ein zweigeteiltes Fenster, das links einen gegliederten Inhaltsbaum und rechts den Text des gewählten Kapitels enthält.
Die einzelnen Kapitel fassen mathematisch zusammengehörige Algorithmen zusammen. Nach einer kurzen Einleitung, welche Probleme im gewählten Kapitel behandelt werden, folgt eine Darstellung des mathematischen Hintergrundes. Dort erhält man wichtige Informationen darüber, in welche spezielle Form man sein Problem gegebenenfalls umformulieren muß.
Es folgt eine Beschreibung der verwendeten Parameter und Fehlermeldungen, die vor allem dann wichtig ist, wenn man selbst eine Hauptroutine zum Aufruf der Funktion schreiben will. Zum Abschluß ist ein Beispielprogramm für so einen Aufruf in der jeweiligen Programmiersprache abgedruckt. Dieses Listing ist aber auch separat als Source-File vorhanden und kann heruntergeladen werden.
Es gibt zwei Möglichkeiten ein lauffähiges C-Programm zu generieren. Die erste besteht darin, ein Hauptprogramm selbst zu schreiben.
Wir empfehlen die zweite Möglichkeit. Sie besteht darin, das mitgelieferte Beispielprogramm so abzuändern, daß es zur eigenen Problemstellung paßt. Über ein gewisses Grundmaß an Programmierkenntnissen in C wird man aber auch hier verfügen müssen.
Folgender Befehl lädt das entsprechende Beispielprogramm und gegebenenfalls ein Datenfile mit Inputdaten herunter, compiliert den Code und führt das Programm mit den Beispieldaten aus:
> nagexample_c nag_kurz´nag_kurz´ ist der kurze NAG-Prozedurname. Der Code und die Daten befinden sich dann in den Files ´nag_kurze.c´ und ´nag_kurze.d´ (´e´ anhängen nicht vergessen). Das ausführbare Beispielprogramm ist das File ´a.out´. Der Code kann nun in einem Editor bearbeitet werden.
Die erste Aufgabe des Benutzers ist, zunächst alle Funktionen, Parameter, Startwerte etc. des Beispieles durch die des eigenen Problemes zu ersetzen.
Das File wird am besten unter einem anderen Namen mit der Endung ´.c´ abgespeichert, am Terminal compiliert und aufgerufen:
> gcc -g file.c -o outputfile -lm -lnagc > outputfileBereits jetzt verfügt man über ein lauffähiges C-Programm, das aber nicht unbedingt der Weisheit letzter Schluß sein muß, da erst hier die Schwierigkeiten auftreten, die mit der Problemstellung und dem verwendeten Algorithmus direkt zu tun haben (siehe die unten angeführten Beispiele).
Der zweite Schritt wird daher darin bestehen, das C-Programm weiter an seine Problemstellung anzupassen, z.B. eine andere als die vorgegebene Genaugikeit zu fordern, oder gewisse Eingaben während des Programmablaufes zu ermöglichen.
Es soll betont werden, daß numerische Algorithmen natürlich keine Garantie für korrekte oder umfassende Ergebnisse sind, und diese daher immer auf ihre Richtigkeit kontrolliert werden sollten.
Der Grundgedanke, der hinter der Numerical Recipies - Programmbibliothek steht, ist ein etwas anderer als jener der NAG: Statt eine Online-Hilfe anzubieten, erwartet man hier, daß der Anwender die Dokumentation ´Numerical Recipies in C´ zur Hand hat. Die jeweiligen Bibliotheksroutinen können dann über unseren Server erhalten werden (Numerical Recipies ). Die Vorgehensweise ist folgende:
Es ist wesentlich darauf zu achten, auch die in der Bibliotheksroutine aufgerufenen Unterprogramme und Bibliotheksfiles bei der Compilierung einzubinden, wie zum Beispiel:
> gcc -g cfile.c linkfiles.c -lm -o cfile.outDas lauffähige Programm heißt nun ´cfile.out´. In unserem Fall ist ´cfile.c´ ein nach unseren Bedürfnissien abgewandeltes Beispielprogramm, ´linkfiles.c´ ist eine Auflistung der vom Hauptprogramm aufgerufenen Unterprogramme, ´-lm´ weist auf <mathlib.h> hin und ´-o´ erklärt ´cfile.out´ zum Outputfile. Trotzdem muß, wie bei jedem numerischen Algorithmus, die Richtigkeit der Ergebnisse je nach Problem hinterfragt bzw. überprüft werden.
Die allgemeine Vorgehensweise ist genau die gleiche, wie die zur Erhaltung einer Routine aus der Bibliothek der Numerical Recipies. Allerdings sind die Beispielprogramme im Verzeichnis 'tst/' des jeweiligen Kapitels. Die Source-Codes finden sich wieder auf unserem Server: Numerical Algorithms
Hier ist Kap.1.2.1. zu Rate zu ziehen.
In diesem Abschnitt stellen wir drei konkrete sehr einfache Probleme vor. Zur Kontrolle wurde das Problem jeweils auch mit Mathematica gelöst. Mit der NAG wurden alle drei Probleme behandelt, mit den Numerical Recipies die letzten zwei und mit den Numerical Algorithms exemplarisch nur das erste.
Es sollen die Nullstellen eines reellen Polynoms 5-ter Ordnung berechnet werden:
| z5 + 2 z4 + 3 z3 + 4 z2 + 5 z + 6 = 0 | (1) |
In[1]:= NSolve[z^5 + 2 z^4 + 3 z^3 + 4 z^2 + 5 z +6 == 0, z]Wir erhalten nach minimaler Rechenzeit folgende fünf Nullstellen:
z = -1.4918 z = -8.0579e-01 - 1.2229 I, z = -8.0579 + 1.2229 I z = 5.5169e-01 - 1.2533 I, z = 5.5169e-01 + 1.2533 I
Wir haben die NAG-Routine ´c02agc()´ verwendet. Entsprechend des obigen Abschnittes wurde ein Beispielprogramm generiert. Dieses kann getestet werden:
> a.out < c02agce.dLediglich die Benutzerschnittstelle wurde von uns um ein wenig Text ergänzt. Die Ordnung und die Koeffizienten werden vom Benutzer eingegeben und nicht aus einem File gelesen. Das Listing findet sich in ´poly.c´ . Folgendes Ergebnis wird für unser konkretes Polynom ausgegeben:
Degree of polynomial = 5 Roots of polynomial z = -1.4918e+00 z = 5.5169e-01 +/- 1.2533e+00 z = -8.0579e-01 +/- 1.2229e+00Die Übereinstimmung mit den Ergebnissen von Mathematica ist optimal. Die Rechenzeit ist minimal.
In unserem Fall wurde die Routine ´fmueller.c´ verwendet. Diese, sowie sämtliche anderen benötigten Dateien finden sich unter Numerical Algorithms . Das Beispielprogramm dazu (´tmueller.c´) wurde in das eigene Verzeichnis kopiert und dann mit den entsprechenden Unterprogrammen und Links zusammen generiert. Es mußten folgende Headerfiles in unser Verzeichnis kopiert werden:´basis.h´, ´vmblock.h´, ´u_proto.h´. Weitere Dateien mußten ins Verzeichnis geschrieben werden: ´fmueller.c´, ´basis.c´, ´vmblock.c´, ´u_proto.c´. Auch hier wurde die Benutzerschnittstelle entsprechend ergänzt bzw. abgeändert. Das abgeänderte Listing findet sich unter ´tmueller.c´ . Folgende Lösungen ergaben sich für unser Problem:
Without scaling z=-8.0578646938903120e-01 +/- I*1.2229047133744098e+00 z= 5.5168546345898151e-01 +/- I*1.25334886022772060e+00 z=-1.4917979881399006e+00 With scaling z=-8.0578646938903131e-01 +/- I*1.2229047133744098e+00 z= 5.5168546345898162e-01 +/- I*1.25334886022772063e+00 z=-1.4917979881399006e+00
Die Option ´scaling´ gestattet es dem aufgerufenen Programm, die zu untersuchende Funktion zu skalieren. Dies kann sich für den Algorithmus unter Umständen vorteilhaft auswirken. Die von Mathematica gelieferten Ergebnisse stimmen verblüffend, um nicht zu sagen geradezu aufs Erstaunlichste, mit den unseren überein.
Es soll folgende gewöhnliche Differentialgleichung dritter Ordnung numerisch gelöst werden:
| y''' - 2 y'' - y' + 2 y = 2 t2 - 6 t + 4 | (2) |
y1' = y2
y2' = y3
y3' = 2 y3 + y2 - 2 y1 + 2 t2 -6 t + 4
Mathematica liefert auf die Eingabe
In[1]:= DSolve[ {y'''[t]- 2 y''[t] - y'[t] + 2 y[t]
== 2 t^2 - 6 t + 4, y[0]==5, y'[0]==-5,
y''[0]==1},y[t],t]
die exakte Lösung:
f(0) = 5 f(0.8) = 0.211166 f(0.2) = 4.00704 f(1.0) = -1.93502 f(0.25)= 3.77541 f(1.25)= -6.05664 f(0.4) = 2.96692 f(1.5) = -12.9076 f(0.5) = 2.3935 f(1.75)= -24.4508 f(0.6) = 1.75963 f(2.0) = -43.9384 f(0.75)= 0.642544Ein numerischer Algorithmus kann natürlich nur Funktionswerte an bestimmten Stellen von t liefern.
Für die Realisierung wurde die NAG-Routine ´d02pcc()´ verwendet. Sie verwendet Runge-Kutta Methoden zur Lösung. Das Listing des generierten und abgeänderten Beispielprogrammes findet sich in ´runge-kutta.c´ . Es mußte die Anzahl der Gleichungen auf drei erhöht werden (das Beispielprogramm war für eine Differentialgleichung zweiter Ordnung geschrieben). Weiters müssen die eigenen Funktionen, d.h. die rechten Seiten der Differentialgleichungen erster Ordnung im Unterprogramm ´f´ realisiert werden. Auch der Bereich, in dem die Lösung berechnet wird, muß neu festgelegt werden.
Das Programm macht zwei Durchläufe für jeweils eine andere Genauigkeit. Es liefert nach kurzer Rechenzeit folgende Ergebnisse:
d02pcc Example Program Results
Calculation with tol = 1.0e-03
t y1 y2
-0.000 5.000 -5.000
0.250 3.755 -5.071
0.500 2.394 -6.001
0.750 0.643 -8.291
1.000 -1.935 -12.795
1.250 -6.056 -20.946
1.500 -12.906 -35.132
1.750 -24.446 -59.314
2.000 -43.930 -100.061
Cost of the integration in evaluations of f is 45
Calculation with tol = 1.0e-04
t y1 y2
-0.000 5.000 -5.000
0.250 3.755 -5.071
0.500 2.394 -6.001
0.750 0.643 -8.291
1.000 -1.935 -12.795
1.250 -6.056 -20.947
1.500 -12.907 -35.135
1.750 -24.450 -59.321
2.000 -43.936 -100.073
Cost of the integration in evaluations of f is 55
Es zeigt sich eine sehr gute Übereinstimmung mit Mathematica.
Wir haben hiezu die Routine `rk4.c` (Runge-Kutta-Methode) verwendet. Das Beispielprogramm dazu ('xrk4.c') wurde in das eigene Verzeichnis kopiert und dann mit den entsprechenden Unterprogrammen und Links zusammengeneriert. Es mußten folgende Headerfiles in unser Verzeichnis kopiert werden: 'nrutil.h', 'nr.h'. Das Listing findet sich in ´xrk4.c´ . Wir haben folgende Ergebnisse erhalten:
x=0.2 f(x)=4.007147 x=0.4 f(x)=2.984089 x=0.6 f(x)=1.854119 x=0.8 f(x)=0.513944 x=1.0 f(x)=-1.177578Wie man sehen kann, gibt es hier eine Abweichung von den Ergebnissen aus der analytischen allgemeinen Lösung (siehe oben, Mathematica!).
Als letztes Problem haben wir ein nichtlineares Gleichungssystem behandelt:
x2 + y2 - 2 = 0
e(x-1) + y3 - 2 = 0
Dieses System bereitet trotz seiner Einfachheit Mathematica schon Probleme. Mit ´NSolve´ kommt man nicht ans Ziel. Mit ´FindRoot´ findet man Lösungen, allerdings nur wenn man die Anfangsbedingungen schon in die Nähe der Lösung setzt. Wir haben das Problem letzten Endes mit Mathematica graphisch gelöst, indem wir es ins dreidimensionale erweiterten und die Flächen miteinander zum Schnitt brachten:z = x2 + y2 - 2
z = e(x-1) + y3 -2
z = 0
Diese Vorgehensweise liefert uns nämlich sofort eine explizite Gleichung, die wir grafisch darstellen können. Wir sehen, daß es zwei Lösungen gibt und wissen wo sie liegen. Nun haben wir mit ´FindRoot´ und den richtigen Anfangswerten die Lösungen gefunden:x1 = 1 ; y1 = 1 x2 = -0.713747; y2 = 1.22089
Es wurde die NAG-Routine ´c05pbc()´ verwendet. Sie löst ein System von nichtlinearen Gleichungen, wobei die ersten Ableitungen der Gleichungen in Form der Jacobi-Determinante vom Benutzer selbst berechnet werden müssen.
Das Listing befindet sich in ´exp.c´ . Es hat sich herausgestellt, daß die gefundenen Lösungen stark von den gewählten Anfangswerten abhängen, daher gibt der Nutzer dem Programm ein Netz von Startwerten vor. Außerdem muß jedes Ergebnis auf seine Richtigkeit überprüft werden, da die Routine selbst den Iterationsfortschritt nicht gut genug bewerten kann. Für einen Startwertbereich von -10 <= x <= 10, -10 <= y <= 10 liefert mir das Programm schon in den ersten sechs Ergebnisse die beiden richtigen Nullstellen. Man sieht auch schön, daß die zweite ``Lösung'' eine falsche Lösung ist, da sie die Gleichung nicht erfüllt:
x[i] f[i]
-------------------------------------
1.00000 -0.00000
1.00000 -0.00000
1.48506 0.20541
-0.00006 -0.37572
-0.71375 -0.00000
1.22089 -0.00000
-0.71375 0.00000
1.22089 0.00000
1.00000 0.00000
1.00000 -0.00000
1.00000 -0.00000
1.00000 0.00000
Die Lösungen stimmen wiederum bestens mit Mathematica überein.
Für diese Aufgabe verwendeten wir die Routine `mnewt.c`. Das entsprechende Beispielprogramm (`xmnewt.c`) wurde kopiert und an unser Problem angepaßt. Als Headerfiles wurden `nr.h`, `nrutil.h` benötigt. Des weiteren wurden `nr.c`, `nrutil.c`, `newt.c`, 'ludcmp.c' und 'lubksb.c' in unser Verzeichnis geschrieben. Das Listing befindet sich in ´xmnewt.c´ . Wir erhielten folgende Ergebnisse:
1.Lsg.: x= 1, y=1 2.Lsg.: x= -0.713747, y=1.220887Wieder zeigt sich eine gute Übereinstimmung mit Mathematica.
This document was generated using the LaTeX2HTML translator Version 97.1 (release) (July 13th, 1997)
Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
The command line arguments were:
latex2html -html_version 3.2 -split 0 -no_navigation texfiles/appl2c.tex.
The translation was initiated by Eherer Christian Gerhard on 5/27/1998