Es ist im Rahmen dieser Lehrveranstaltung völlig unmöglich, das sehr umfangreiche Gebiet der numerischen Behandlung von Randwertproblemen (RWP) umfassend zu behandeln. Ich habe daher aus der großen Mannigfaltigkeit der Probleme jene Gruppe herausgegriffen, die in der Praxis von Physik und Technik die wichtigste Rolle spielt, nämlich
Dieses Problem besteht in mathematischer Hinsicht aus der gewöhnlichen
linearen Differentialgleichung (Dgl.) zweiter Ordnung
Man nennt ein Problem (9.1-9.3) homogen, wenn
gilt:
Der wichtigste Schritt beim Differenzenverfahren besteht darin, daß
das kontinuierliche Intervall, innerhalb dessen die Lösungsfunktion
gesucht wird, durch Aufteilung in N (gleich-breite)
Subintervalle diskretisiert wird. Man erhält auf diese Weise die
äquidistanten Stützpunkte
An diesen Stützpunkten werden nun alle in der Dgl. und in den Randbedingungen auftretenden Differentialquotienten durch geeignete Differenzenquotienten ersetzt.
Es gibt nun in der Literatur zahlreiche Formeln für solche
Differenzenquotienten (s. z. B. [20], S. 267 und 267).
In vielen Fällen sind jedoch die einfachsten
derartigen Formeln die numerisch praktikabelsten. Aus diesem Grund
werden in diesem Kapitel aussschließlich die sogenannten
'symmetrischen Dreipunktsformeln'9.1
Die grafische Interpretation von (9.5) und (9.6)
ist besonders einfach: Die Steigung der Kurve im Punkt wird
approximiert durch die Steigung der Verbindungsgeraden, die durch den rechten
und den linken Nachbarpunkt hindurchgeht (s. Abb. 9.1).
Entsprechend wird die zweite Ableitung der Funktion im Punkt
durch die zweite Ableitung jener Parabel angenähert, die durch den Punkt
selbst sowie durch seinen rechten und linken Nachbarn
hindurchgeht.
Ersetzt man nun die erste und zweite Ableitung in (9.1) durch
die Ausdrücke (9.5) und (9.6), so erhält man für
die Differenzengleichungen
Diese Randbedingung beschreibt das Verhalten der Lösungsfunktion
an der Stelle d.h. für :
Die Behandlung der Randbedingung (9.3) geschieht äquivalent
zu den obigen Überlegungen und führt unter der Bedingung
zur folgenden neuen letzten (d. h. ()-ten) Differenzengleichung
Auf diese Weise hat man das gegebene RWP in ein inhomogenes, lineares
Gleichungssystem von Gleichungen für die Unbekannten
.
Die erste Gleichung dieses Systems ist gegeben durch (9.9) bzw.
(9.10), die zweite bis -te Gleichung hat die Form (9.7),
und die ()-te Gleichung lautet (9.11) bzw. (9.12).
Wie bereits erwähnt, stellen die Lösungen dieses Systems die
approximativen Werte der Lösungsfunktion des RWP an den
Stützpunkten
dar.
Dieses System muß also (i. a. mit numerischen Mitteln) gelöst werden.
Eine Analyse der Differenzengleichungen zeigt Ihnen sofort, daß die
Koeffizientenmatrix des linearen Gleichungssystems eine tridiagonale
Matrix ist. Diese einfache Form ist der speziellen Wahl der finiten
Elemente (9.5) und (9.6) zu danken!
Wie Sie sich erinnern, gibt es für die Lösung von linearen Gleichungssystemen mit tridiagonaler Matrix sehr effizente numerische Methoden. Eine solche finden Sie im Abschnitt 2.6.1 dieses Skriptums.
Wiederholt man die Berechnung der mit der doppelten Zahl von
Subintervallen (d.h.:
), so entspricht wegen der
äquidistanten Stützpunkte jeder zweite Punkt der '2-Rechnung' einem
Punkt der '-Rechnung'. Für diese Punkte kann man nun eine sehr
einfache und äquidistante Reduktion des beim Differenzenverfahren
auftretenden Verfahrensfehlers vornehmen; diese Methode beruht auf dem
in der Literatur beschriebenen Faktum, daß der Verfahrensfehler
für die Form
Dieses Programm löst das RWP (9.1)-(9.3) mittels des Differenzenverfahrens. Die Hauptpunkte des Programmes sind:
INPUT-Parameter:
OUTPUT-Parameter:
INTERNE FELDER:
Das Programm DIFF1 benötigt die Prozeduren FCT1 und TRID.
=16cm
Struktogramm 29DIFF1(ANF,AEND,ALP0,ALP1,BET0,BET1,GAMMA,
DELTA,N0,Y)T=1(1)2N:=T*N0
H:=(AEND-ANF)/N
HQ:=H*HI=2(1)N X:=ANF+(I-1)*H FCT1(X,R,S,Q,F) A(I):=R/HQ-S/2.0/H
B(I):=Q-2.0*R/HQ
C(I):=R/HQ+S/2.0/H
INHOM(I):=FALP1=0.0B(1):=1.0
C(1):=0.0
INHOM(1):=GAMMA/ALP0 FCT1(ANF,R,S,Q,F)B(1):=Q-2.0*R/HQ - ALP0/ALP1*(S-2.0*R/H)
C(1):=2.0*R/HQ
INHOM(1):=F - GAMMA/ALP1*(S-2.0*R/H)BET1=0.0A(N+1):=0.0
B(N+1):=1.0
INHOM(N+1):=DELTA/BET0 FCT1(AEND,R,S,Q,F)A(N+1):=2.0*R/HQ
B(N+1):=Q-2.0*R/HQ - BET0/BET1*(S+2.0*R/H)
INHOM(N+1):=F - DELTA/BET1*(S+2.0*R/H) NP1:=N+1TRID(A,B,C,INHOM,NP1,Y) I=1(1)N0+1T=1 INDEX:=IINDEX:=2*I-1 YY(I,T):=Y(INDEX) H:=(AEND-ANF)/N0 I=1(1)N0+1 X:=ANF+(I-1)*H
Y(I):=(4.0*YY(I,2)-YY(I,1))/3.0(return)
Das nun folgende Beispiel soll Ihnen die Leistungsfähigkeit von DIFF1
demonstrieren. Es ist der Referenz [19], S. 252f entnommen und
lautet:
Die exakte Lösung dieses RWP lautet:
Der Benutzer hat also eine Prozedur FCT1 zu schreiben, welche z. B. die Form hat:
. PROCEDURE FCT1(x: real; VAR r,s,q,f: real); BEGIN r:=1.0; s:=-2.0*x; q:=-2.0; f:=-4.0*x END;
Die weiteren Parameter lauten:
. anf=0.0 aend=1.0 alp0=1.0 alp1=-1.0 bet0=1.0 bet1=0.0 gamma=0.0 delta=1.0 + e
Die Auswertung mittels DIFF1 wurde für N0=10 Subintervalle durchgeführt,
d.h. das Programm DIFF1 führte die Rechnung einmal mit 10 Subintervallen
und ein zweitesmal mit 20 Subintervallen durch; dann wurde die
Korrekturformel (9.14) angewendet.
Die Ergebnisse sehen Sie in den folgenden Tabellen sowie in der
Abb. 9.2.
. Testergebnisse fuer DIFF1: absolute Fehler zwischen Naeherungswerten und exakten Werten: x abs. Fehler n=10 n=20 korr. mittels (9.14) 0.0 1.9162406497E-03 4.7755794913E-04 -2.0029510779E-06 0.1 2.0768600371E-03 5.1777355657E-04 -1.9219369278E-06 0.2 2.1803285636E-03 5.4365274082E-04 -1.9058661564E-06 0.3 2.2261855847E-03 5.5508881815E-04 -1.9434355636E-06 0.4 2.2088322276E-03 5.5069313748E-04 -2.0198913262E-06 0.5 2.1176171067E-03 5.2782016428E-04 -2.1121486498E-06 0.6 1.9371537346E-03 4.8265275109E-04 -2.1809100872E-06 0.7 1.6485163578E-03 4.1051038716E-04 -2.1582709451E-06 0.8 1.2326037795E-03 3.0670399065E-04 -1.9292747311E-06 0.9 6.7823138306E-04 1.6857800802E-04 -1.3064491213E-06 1.0 1.0913936421E-11 1.0913936421E-11 1.0913936421E-11
Im folgenden soll nun die numerische Auswertung des homogenen RWP
(9.1)-(9.3) diskutiert werden, und zwar in der Form eines
Eigenwertproblems:
Für ein solches Problem gibt es immer die triviale Lösung sowie für bestimmte Werte von nicht-triviale Lösungen. Man spricht von den Eigenwerten und Eigenfunktionen (Eigenlösungen) des homogenen RWP.
Die mathematische Behandlung von (9.15)-(9.17) ist
identisch mit der Vorgangsweise in Abschnitt 9.2. Wieder wird das RWP
durch Diskretisierung des Intervalles [,] in ein lineares, diesmal
natürlich homogenes Gleichungssystem umgewandelt. Unter der
Voraussetzung
und
erhält man
In jedem Fall hat man nun die Eigenwerte einer tridiagonalen Matrix zu ermitteln. Dies kann z. B. mittels des Hyman-Verfahrens in Zusammenarbeit mit einem Nullstellen-Suchprogramm wie INTSCH geschehen. Genaueres zu diesem Thema finden Sie im Abschnitt 7.5.5 dieses Skriptums.
Ohne weitere theoretische Ableitungen soll hier eine Formel für die
Fehlerkorrektur der so erhaltenen Näherungen der Eigenwerte des RWP
angegeben werden. Bezeichnet man mit bzw. mit
einen Eigenwert der tridiagonalen Matrix in (9.21), wobei das
Intervall in bzw. Subintervalle geteilt wurde, so
kann der Verfahrensfehler bzgl. dieses Eigenwertes mittels der Formel
Es soll hier aber ein ganz wesentlicher Nachteil dieser Korrekturformel nicht verschwiegen werden: man hat keinerlei Information über die Genauigkeit der korrigierten Eigenwerte!
Dieses Programm berechnet einige reelle Eigenwerte des homogenen RWP
(9.15)-(9.17) mittels des Differenzenverfahrens sowie des
Verfahrens von
Hyman (s. Abschnitt 7.5.5).
Die Hauptpunkte des Programmes sind:
INPUT-Parameter:
Die folgenden 4 Input-Parameter beziehen sich auf die von DIFF2 aufgerufene Prozedur INTSCH:
OUTPUT-Parameter:
INTERNE FELDER:
Das Programm DIFF2 benötigt die Prozeduren FCT1, INTSCH
und HYMAN.
Anmerkung:
Um die Verfahrensfehler von INTSCH möglichst klein zu halten, sollte
die Fehlerschranke GENEIG sehr klein gewählt werden,
z.B.: GEN = 0.0000001.
Um dieser Anforderung gerecht zu werden, sollte DIFF2 - wenn möglich -
mit doppelter Genauigkeit gerechnet werden.
=16cm
Struktogramm 30DIFF2(ANF,AEND,ALP0,ALP1,BET0,BET1,ANFEIG,
ENDEIG,HEIG,GENEIG,N0,ANZEIG,EIGW,FEHL1)T=1(1)2N:=T*N0
H:=(AEND-ANF)/N
HQ:=H*HI=2(1)N X:=ANF+(I-1)*H FCT1(X,R,S,Q) A(I):=R/HQ-S/2.0/H
B(I):=Q-2.0*R/HQ
C(I):=R/HQ+S/2.0/HALP1=0.0IANF:=2IANF:=1 FCT1(ANF,R,S,Q)B(1):=Q-2.0*R/HQ - ALP0/ALP1*(S-2.0*R/H)
C(1):=2.0*R/HQBET1=0.0IEND:=NIEND:=N+1 FCT1(AEND,R,S,Q)A(N+1):=2.0*R/HQ
B(N+1):=Q-2.0*R/HQ - BET0/BET1*(S+2.0*R/H)N:=IEND-IANF+1ALP1=0.0I=1(1)NA(I):=A(I+1)
B(I):=B(I+1)
C(I):=C(I+1)......INTSCH(ANFEIG,ENDEIG,HEIG,GENEIG,ANZMAX,EIGW,ANZ(T))ANZ(T)=0print: DIFF2 'keine Nullst. in INTSCH'I=1(1)ANZ(T)EWERTE(I,T):=EIGW(I) ANZEIG:=ANZ(2)ANZ(1)=ANZ(2)I=1(1)ANZEIG EIGW(I):=(4.0*EWERTE(I,2) - EWERTE(I,1))/3.0 FEHL1:=FALSEprint: DIFF2 'Keine Korrektur'I=1(1)ANZEIGEIGW(I):=EWERTE(I,2)FEHL1:=TRUE(return)
Als Test für die Leistungsfähigkeit des eben besprochenen Algorithmus
soll nun das lineare, homogene EWP zweiter Ordnung
Die Eigenwerte dieses Problems gehorchen der
Formel
Im Bereich
gibt es also die 5 Eigenwerte
Die INPUT-Parameter und die Prozedur FCT1 lauten in diesem Fall:
. anf:=1.57079633 aend:=3.14159265 alp0:=1.0 alp1:=0.0 bet0:=0.0 bet1:=1.0 n0:=50 anfeig:=0.0 endeig:=100.0 heig:=1.0 geneig:=0.0000001 PROCEDURE FCT1(x: double; VAR r,s,q: double); BEGIN r:=-1.0; IF(abs(x-PI) < 1.e-15)THEN s:=0.0 ELSE s:=-cos(x)/sin(x); q:=0.0; END;
Anmerkung zu FCT1:
Bei der Auswertung der Koeffizientenfunktion ist darauf zu achten,
daß sie für eine Singularität hat. In der obigen
Routine wird für diesen Fall einfach Null gesetzt. Dies ist aber nur deshalb
möglich, weil der in der Dgl. (9.23) vorkommende Term
Löst man das vorliegende Problem mit Hilfe der im vorigen Abschnitt präsentierten Prozedur DIFF2, so erhält man das folgende Ergebnis:
. Zahl der Subintervalle <= 50 GENEIG = 0.0000001 N=50 N=100 korr. nach exakt (9.22) 1 1.99901 1.99975 2.00000 2 2 11.97850 11.99463 12.00001 12 3 29.88705 29.97181 30.00006 30 4 55.64148 55.91053 56.00021 56 5 89.12707 89.78212 90.00047 90
Anschließend soll hier noch die in der Praxis häufig verwendete
Methode besprochen werden, das homogene RWP (9.15-9.17)
durch Überführung in ein Anfangswertproblem zu lösen.
Zu diesem Zweck bringt man (9.15) in die Form
Man benötigt für die Dgl. zweiter Ordnung natürlich
zwei Anfangsbedingungen.
Diese folgen sofort aus der ersten Randbedingung (für ), welche die
Form
Ist hingegen , so bleibt für nur mehr die
Möglichkeit
Die Gleichungen (9.25)-(9.29) definieren ein Anfangswertproblem, das z. B. mittels der im Kapitel 8 besprochenen Verfahren gelöst werden kann.
Die Eigenwerte sind noch unbekannt. Man kann nun wie folgt vorgehen: Man löst das Cauchy-Problem (9.25)-(9.29) im Intervall für ein beliebiges 'Probe'-.
Nun überprüft man, ob die erhaltene Lösungsfunktion
an der Stelle des rechten Randes, also für , die
zweite Randbedingung (9.17) erfüllt, d.h. ob gilt:
Abb. 9.3 zeigt das Prinzip dieser Methode: und 'treffen daneben', 'trifft' die Randbedingung bei und stellt daher einen Eigenwert dar. Die Formulierung 'treffen' bzw. 'daneben treffen' ist bewußt gewählt, um Ihnen den Namen dieser Methode (shooting method) klar zu machen.
Die shooting method hat jedoch einen gravierenden Nachteil: um das Eigenwertproblem mit genügender Genauigkeit zu lösen, müssen sehr viele 'Probe'- verwendet werden d. h. es müssen sehr viele Cauchy-Probleme gelöst werden, was natürlich häufig große Rechenzeiten bedeutet.
Aus diesem Grund wird die shooting method in der Praxis vor allem in Kombination mit einer sehr schnellen und dennoch genauen Lösungsmethode für das Cauchy-Problem kombiniert, die auf Numerov zurückgeht.
Diese Methode ist nur anwendbar, wenn es gelingt, die Dgl. des homogenen RWP
auf die Form
Die Numerov-Methode ist wieder eine Differenzenmethode, d.h. man teilt
das gegebene Intervall in Subintervalle mit den Stützpunkten
Nun noch einige wichtige Anmerkungen zum shooting-Numerov-Verfahren:
Als ersten Schritt des Numerov-Prozesses muß die Glg. (9.36) für
ausgewertet werden. Dies führt zur Formel
Es wäre nun ausgesprochen ungeschickt, bei jedem Aufruf des Numerov-Programms die Potentialwerte aufs Neue zu berechnen, denn wenn die Potentialfunktion kompliziert ist, würde das eine Menge Rechenzeit kosten. Viel besser ist es, das Feld der Werte , , ..., vor dem ersten Numerov-Aufruf ein für alle Mal zu berechnen und dieses Feld über die Parameterliste an Numerov zu übergeben. Eine solche Vorgangsweise ist in der folgenden Beschreibung des Numerov-Programms berücksichtigt.
Die Prozedur NUMEROV führt eine numerische Integration der Dgl. (9.31) durch, wobei diese Gleichung den Eigenwert-Parameter enthält.
INPUT-Parameter:
OUTPUT-Parameter:
=14cm
Struktogramm 31NUMEROV(ANF,AEND,Y1,Y2,V1Y1,EIGW,POTF,N,YF) H:=(AEND-ANF)/N
HH12:=H*H/12.0 YF(1):=Y1
YF(2):=Y2
YLL:=Y1
YL:=Y2
KL:=EIGW-POTF(2)
K:=EIGW-POTF(3) Y:=(2.0*(1.0-5.0*HH12*KL)*YL - YLL - HH12*(EIGW*YLL-V1Y1)) / (1.0+HH12*K) YF(3):=YI=4(1)N+1YLL:=YL
YL:=Y
KLL:=KL
KL:=K
K:=EIGW-POTF(I) Y:=(2.0*(1.0-5.0*HH12*KL)*YL - YLL - HH12*KLL*YLL) / (1.0+HH12*K) YF(I):=Y(return)
Im folgenden soll nun diese Methode auf ein Problem der Metallphysik angewendet werden. Gesucht ist die Grundzustandsenergie und die entsprechende Wellenfunktion eines Valenzelektrons in metallischem Natrium in der Wigner-Seitz-Näherung.
Die Wigner-Seitz-Näherung besteht darin, daß man die i. a. geometrisch
recht komplizierte Form der Einheitszelle des Kristalls durch eine
volumsgleiche Kugel mit dem Radius ersetzt. Innerhalb der Kugel
wird das Kristallpotential durch ein radialsymmetrisches Potential
angenähert:
In atomaren Einheiten (Längen in Bohr, Energien in
Rydberg) gilt außerdem
, und man erhält
Der Wigner-Seitz-Radius für das Natrium-Metall beträgt
3.9405 Bohr,
und als Potential innerhalb der WS-Zelle verwendeten
Wigner und Seitz in ihrer Originalarbeit das im folgenden
angegebene Prokofjew-Potential:
. FUNCTION POT(r: double): double; VAR q : double; BEGIN IF(r <= 0.01)THEN q:=11.0*r ELSE IF(r <= 0.15)THEN q:=-26.4*r*r+11.53*r-0.00264 ELSE IF(r <= 1.00)THEN q:=-2.84*r*r+4.46*r+0.5275 ELSE IF(r <= 1.55)THEN q:=1.508*r*r-4.236*r+4.876 ELSE IF(r <= 3.30)THEN q:=0.1196*r*r+0.2072*r+1.319 ELSE IF(r <= 6.74)THEN q:=0.0005*r*r+0.9933*r+0.0222 ELSE q:=r; pot:=-2.0*q/r/r END;
Die Bedingungen, die eine Eigenfunktion von (9.39) erfüllen muß, sind die folgenden:
Nun noch ein wichtiges Detail:
Die Funktion in (9.39) divergiert offenbar für wegen
Anmerkung: wie Sie sehen, braucht man zur Auswertung der Glg. (9.42) noch den Wert rechts neben dem Ende des eigentlichen Integrationsintervalls. Um diesen Wert zu erhalten, müssen Sie im Programm Numerov die Zählschleife von I=4(1)N+1 auf I=4(1)N+2 erweitern, und auch die Felder POTF und YF müssen dementsprechend um 1 Speicherplatz erweitert werden.
Die folgende Tabelle zeigt die 'Grobsuche' für das konkrete Problem:
. TESTRECHNUNG SHOOTING-NUMEROV 1999 Wigner-Seitz-Radius [Bohr] = 3.9405 Zahl der Subintervalle = 400 eigw y(b) b*y`(b) diff (soll Null sein!) -0.7500 0.3844 0.6288 -0.2444 -0.6500 0.3259 0.3886 -0.0627 -0.6400 0.3203 0.3664 -0.0461 -0.6300 0.3147 0.3445 -0.0298 -0.6200 0.3092 0.3229 -0.0138 -0.6100 0.3037 0.3017 0.0020 -0.6000 0.2982 0.2808 0.0175 -0.5900 0.2928 0.2601 0.0327 -0.5800 0.2875 0.2398 0.0476 -0.5700 0.2821 0.2198 0.0623 -0.5600 0.2769 0.2002 0.0767 -0.5500 0.2716 0.1808 0.0908 -0.4500 0.2215 0.0036 0.2179
Der gesuchte Energie-Eigenwert für das Natrium-Valenzelektron
liegt also im Bereich zwischen -0.62 und -0.61 Rydberg. Eine genauere
Rechnung führt zum Ergebnis
In bezug auf Randwertprobleme gewöhnlicher Differentialgleichungen hat die NAG-Bibliothek einiges zu bieten, wie die folgende Tabelle zeigt:
Sie finden in dieser Tabelle bekannte Begriffe wie shooting method und finite-difference method, wobei die Differenzenmethode nicht nur auf lineare, sondern auch auch nicht-lineare RWP angewendet werden kann. Die Kollokationsmethode nach Chebyshev ist für lineare Probleme ebenfalls sehr leistumgsfähig; es handelt sich dabei um eine Entwicklung der Lösungsfunktion nach Chebyshev-Polynomen und um eine numerische Berechnung der Polynom-Koeffizienten.
Obwohl in dieser LV. die Besprechung der numerischen Behandlung von
partiellen Differentialgleichungen (leider) keinen Platz hat, möchte
ich abschließend zumindest einige Hinweise auf das sehr umfangreiche
Software-Angebot in bezug auf diese Probleme geben:
Einen sehr guten Überblick über p.d. und kommerizielle Software finden
Sie in den schon häufig zitierten Büchern von Ch. Überhuber
([21],[22]). So ist z. B. in
[21], S. 331ff eine Fallstudie 'Software für elliptische
partielle Differentialgleichungen' enthalten, die einschlägige
Programme aus den (p.d.) Bibliotheken ELLPACK und TOMS sowie aus NAG und
IMSL vorstellt.