Numeric Library / Referenz

 Omikron Basic im Internet: http://www.berkhan.de

Allgemeines -blättern- Inhaltsverzeichnis


2. Numeric Library Referenz

2.1 An- und Abmelden der Library
2.2 Spezielle Funktionen
2.3 Berechnung von Nullstellen und Extremwerten
2.4 Berechnung von Eigenwerten und Lösen von Gleichungssystemen.
2.5 Interpolation
2.6 Approximation von Funktionen durch Funktionensysteme
2.7 Numerische Differentiation
2.8 Numerische Integration
2.9 Lösung von Differentialgleichungen
2.10 Graphische Ausgabe
 
 
 
 


In diesem Abschnitt sollen die Prozeduren und Funktionen der Numeric Library erläutert werden. Der theoretische Hintergrund der einzelnen Befehle kann jedoch hier nicht erläutert werden. Für den sinnvollen Einsatz der Library ist es daher empfehlenswert, zusätzliche Literatur zu verwenden. Dabei möchten wir besonders auf das Buch "Numerical Recipes" hinweisen, das viele Aspekte der numerischen Mathematik sehr ausführlich und praxisbezogen darstellt.
 
 




2.1 An- und Abmelden der Library
 
Num_Init
Rufen Sie diese Prozedur zu Beginn Ihres Programms einmal auf. Vorher kann die Numeric Library nicht benutzt werden.
Wichtig: Diese Prozedur verändert den DATA-Zeiger. Wenn Sie also Ihre eigenen Daten einlesen wollen nachdem Sie Num_Init aufgerufen haben, müssen Sie vorher den DATA-Zeiger mit RESTORE auf die gewünschten Daten setzen.
 
 
Num_Exit
Rufen Sie diese Prozedur am Ende Ihres Programms einmal auf. Danach können Sie die Numeric Library nicht mehr benutzen.
 
 
Numeric
Es wird nur eine Copyright-Meldung der Numeric Library ausgegeben.

 
 
Num_Mode M
M Betriebsart:
M=0 : Fehlermeldungen anzeigen.
M=1 : Fehlermeldungen unterdrücken.
Mit dieser Prozedur können Sie die Betriebsart der Numeric Library einstellen. Wenn ein Fehler auftritt, wird dieser normalerweise in einer Alertbox angezeigt. Mit M=1 können Sie die Anzeige unterdrücken, z.B. um selbst eine Fehlerbehandlung durchzuführen.

 
 
FN Num_Error
Liefert die Nummer des zuletzt aufgetretenen Fehlers. Wenn Sie die automatischen Fehlermeldungen mit Num_Mode 1 abgeschaltet haben, können Sie diese Funktion verwenden, um herauszufinden, welcher Fehler aufgetreten ist.

Im einzelnen sind folgende Fehlernummern möglich:

 0 : Kein Fehler.
 1 : FN Li#(N,X#) arbeitet nur mit N=1 oder N=2.
 2 : FN Bessel#(N,X#) arbeitet nur mit N>=0.
 3 : FN Weber#(N,X#) arbeitet nur mit N>=0.
 4 : FN Mod_Bessel#(N,X#) arbeitet nur mit N>=0.
 5 : FN Macdo#(N,X#) arbeitet nur mit N>=0.
 6 : FN Zero# konnte keine Nullstelle finden.
 7 : FN Minimum# konnte kein Minimum finden.
 8 : FN Maximum# konnte kein Maximum finden.
 9 : PROC Eigen arbeitet nur mit symmetrischen Matrizen.
10 : Eigenwertberechnung in PROC Eigen konvergiert nicht.
11 : Das Gleichungssystem in PROC Gauss hat keine Lösung, da die Matrix singulär
   ist.
12 : PROC Gauss_Seidel konvergiert nicht.
13 : Die Anzahl der Gleichungen ist nicht sinnvoll in PROC Band.
14 : PROC Newton_Sys konvergiert nicht.
15 : Zu wenig Stützpunkte in PROC Spline_Int.
16 : PROC Gauss_Fit konvergiert nicht.
17 : Fehler in PROC Fft. Verwenden Sie für N nur gerade, positive Zahlen!
18 : PROC Mean_Approx konvergiert nicht.
19 : Zu wenig Daten in PROC Derive (N>=7).
20 : Zu wenig Daten in FN Newton_Cotes (N>=5).
21 : FN Romberg konvergiert nicht.
22 : FN Romberg_2 konvergiert nicht.
23 : Ungültiger Parameter in FN Gauss_Leg# (2<=N<=15).
 
 


2.2 Spezielle Funktionen
 
Achtung: Nicht alle Funktionen können exakt berechnet werden. Viele Funktionen haben an einigen Stellen Singularitäten, müssen über Polynom-Interpolationen oder rekursiv berechnet werden. Da die Genauigkeit von Fließkomma-Zahlen begrenzt ist, kann es zur Aufschaukelung von Fehlern kommen, so daß die Ergebnisse völlig falsch werden. Z.B. liefert die modifizierte Bessel-Funktion für N=100 in der Umgebung von X#=0 völlig falsche Ergebnisse. Für normale Anwendung sind die Funktionen jedoch ohne Probleme zu verwenden und liefern recht genaue Funktionswerte.

Tip: Setzen Sie in der Testphase das Steuerwort COMPILER "DEBUG ON". Dann werden Sie auf problematische Stellen aufmerksam gemacht, da alle Prozessor-Exceptions gemeldet werden.
 
 

FN Atn2#(X#,Y#)
X# Fließkommazahl -INF<X#<INF
Y# Fließkommazahl -INF<Y#<INF AND Y#<>0
Berechnet den zweiwertigen Arcus Tangens (ARCTAN(X#/Y#)). Dabei wird je nach Vorzeichen von X# und Y# der Quadrant so ausgewählt, daß das Ergebnis im Intervall [-PI,PI] liegt. Diese Funktion eignet sich daher z.B. sehr gut, um kartesische Koordinaten in Polarkoordinaten umzurechnen.

Hinweis: Die BASIC Funktion ARCTAN(X#/Y#) liefert nur Ergebnisse im Intervall [-PI/2,PI/2] und stimmt daher nur für positive X# mit FN Atn2# überein.
 
 
FN Hypot#(X#,Y#)
X# Fließkommazahl -INF<X#<INF
Y# Fließkommazahl -INF<Y#<INF
Berechnet SQR(X#^2+Y#^2). Die Berechnung erfolgt direkt, also ohne erst die Summe der Quadrate zu bilden. Dadurch kommt es nicht so leicht zu einem Overflow, wie bei einer direkten Auswertung der Formel.
 
 
FN Compound#(X#,Y#)
X# Fließkommazahl 0<X#<INF
Y# Fließkommazahl 0<Y#<INF
Berechnet (1+X#)^Y#. Für kleine X# liefert diese Funktion genauere Werte als eine direkte Auswertung der Formel. Diese Funktion ist besonders für kaufmännische Anwendungen wichtig.
 
 
FN Annuity#(X#,Y#)
X# Fließkommazahl 0<X#<INF
Y# Fließkommazahl 0<Y#<INF
Berechnet (1-(1+X#)^(-Y#)/X#. Für kleine X# liefert diese Funktion genauere Werte als eine direkte Auswertung der Formel. Diese Funktion ist besonders für kaufmännische Anwendungen wichtig.
 
 
FN Exp2#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet 2^X#.
 
 
FN Expm1#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet EXP(X#)-1. Für kleine X# liefert diese Funktion genauere Werte als eine direkte Auswertung der Formel.
 
 
FN Lnp1#(X#)
X# Fließkommazahl -1<X#<INF
Berechnet LN(X#+1). Für kleine X# liefert diese Funktion genauere Werte als eine direkte Auswertung der Formel.
 
 
FN Erf#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet (2/PI)*INTEGRAL[0,X#]EXP((-t)^2)dt. Dieses Integral ist auch als Fehlerfunktion bekannt und hat in der Statistic und Wahrscheinlichkeitsrechnung eine große Bedeutung.
 
 
FN Erfc#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet 1-FN Erf#(X#). Diese Funktion ist auch als komplementäre Fehlerfunktion bekannt und hat in der Statistic und Wahrscheinlichkeitsrechnung eine große Bedeutung. Für große positive Zahlen (ab 10) liefert diese Funktion genauere Ergebnisse als eine direkte Auswertung der Formel.
 
 
FN Gamma#(X#)
X# Fließkommazahl -INF<X#<INF AND X#<>0,-1,-2,-3 ...
Berechnet INTEGRAL[0,INF]EXP(-t)*t^(X#-1)dt. Dieses Integral ist auch als Gammafunktion bekannt. Diese Funktion hat Singularitäten bei null und allen negativen ganzen Zahlen. Mit der Fakultätsfunktion des Omikron Basics gibt es den Zusammmenhang: FN Gamma#(X#)=FACT(X#-1).
 
 
FN Lgamma#(X#)
X# Fließkommazahl -INF<X#<INF AND X#<>0,-1,-2,-3 ...
Berechnet LN(ABS(FN Gamma#(X#)). Weil die Gammafunktion sehr schnell gegen unendlich strebt, ist es häufig von Vorteil, mit dem Logarithmus der Gammafunktion zu arbeiten, um einen Overflow zu vermeiden.
 
 
FN Hermite#(N,X#)
N Grad des Polynoms 0<=N<INF
X# Fließkommazahl -INF<X#<INF
Berechnet den Wert des Hermite-Polynoms vom Grade N an der Stelle X#.
 
 
FN Cheby#(N,X#)
N Grad des Polynoms 0<=N<INF
X# Fließkommazahl -INF<X#<INF
Berechnet den Wert des Tschebyscheff-Polynoms vom Grade N an der Stelle X#. Verwechseln Sie diese Funktion nicht mit der zur Tschebyscheff-Approximation, die weiter unten besprochen wird.
 
 
FN Legendre#(N,X#)
N Grad des Polynoms 0<=N<INF
X# Fließkommazahl -INF<X#<INF
Berechnet den Wert des Legendre-Polynoms vom Grade N an der Stelle X#.
 
 
FN Legendre_D#(N,X#)
N Grad des Polynoms 0<=N<INF
X# Fließkommazahl -INF<X#<INF AND X#<>-1 AND X#<>1
Berechnet den Wert der Ableitung des Legendre-Polynoms vom Grade N an der Stelle X#. Diese wird auch als Gegenbauer-Polynom bezeichnet.
 
 
FN Horner#(&A#(),N,X#)
A#(0:N) Koeffizienten der ganzrationalen Funktion.
N Grad der Funktion 0<=N<INF
X# Fließkommazahl -INF<X#<INF
Berechnet den Wert der ganzrationalen Funktion
A#(N)*X#^N+A#(N-1)*X#^(N-1)+...+A#(1)*X#+A#(0) an der Stelle X#.
 
 
FN Horner_D#(&A#(),N,X#)
A#(0:N) Koeffizienten der ganzrationalen Funktion.
N Grad der Funktion 0<=N<INF
X# Fließkommazahl -INF<X#<INF
Berechnet den Wert der Ableitung der ganzrationalen Funktion
A#(N)*X#^N+A#(N-1)*X#^(N-1)+...+A#(1)*X#+A#(0)an der Stelle X#.
 
 
FN Beta#(X#,Y#)
X# Fließkommazahl -INF<X#<INF AND X#<>0,-1,-2,-3 ...
Y# Fließkommazahl -INF<Y#<INF AND X#<>0,-1,-2,-3 ...
Berechnet die Beta-Funktion. Es gilt die Beziehung:
FN Beta#(X#,Y#)=FN Gamma#(X#)*FN Gamma#(Y#)/FN Gamma#(X#+Y#)
 
 
FN Zeta#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet die Riemannsche Zeta-Funktion.
 
 
FN Kappa#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet die Kappa-Funktion.
 
 
FN Eta#(X#)
X# Fließkommazahl -INF<X#<INF
Berechnet die Eta-Funktion.
 
 
FN Li#(N,X#)
N Ordnung des Polylogarithmus 1<=N<=2
X# Fließkommazahl -1<=X#<INF
Berechnet den Wert des Polylogarithmus der Ordnung N an der Stelle X#.
 
 
FN Chi#(N,X#)
N Ordnung der Chi-Funktion 1<=N<=2
X# Fließkommazahl -1<=X#<INF
Berechnet den Wert der Chi-Funktion der Ordnung N an der Stelle X#.
 
 
FN Bessel#(N,X#)
N Ordnung der Bessel-Funktion 0<=N
X# Fließkommazahl -INF<=X#<INF
Berechnet den Wert der Bessel-Funktion der Ordnung N an der Stelle X#.
Im Ordner DEMO finden Sie das Programm "Bessel.BAS", mit dem Sie alle Besselfunktionen graphisch darstellen können.
 
 
FN Weber#(N,X#)
N Ordnung der Weber-Funktion 0<=N
X# Fließkommazahl 0<X#<INF
Berechnet den Wert der Weber-Funktion der Ordnung N an der Stelle X#. Die Weber-Funktion wird in der Literatur auch häufig als Bessel-Funktion 2. Art bezeichnet.
Im Ordner DEMO finden Sie das Programm "Bessel.BAS", mit dem Sie alle Besselfunktionen graphisch darstellen können.
 
 
FN Mod_Bessel#(N,X#)
N Ordnung der modifizierten Bessel-Funktion 0<=N
X# Fließkommazahl -INF<=X#<INF
Berechnet den Wert der modifizierten Bessel-Funktion der Ordnung N an der Stelle X#.
Im Ordner DEMO finden Sie das Programm "Bessel.BAS", mit dem Sie alle Besselfunktionen graphisch darstellen können.
 
 
FN Macdo#(N,X#)
N Ordnung der MacDonald-Funktion 0<=N
X# Fließkommazahl 0<X#<INF
Berechnet den Wert der MacDonald-Funktion der Ordnung N an der Stelle X#. Die MacDonald-Funktion wird in der Literatur auch häufig als modifizierte Bessel-Funktion 2. Art bezeichnet.
Im Ordner DEMO finden Sie das Programm "Bessel.BAS", mit dem Sie alle Besselfunktionen graphisch darstellen können.
 
 





2.3 Berechnung von Nullstellen und Extremwerten
 
Nullstellen und Extremwerte sind in den meisten Fällen nur schwer analytisch zu finden. Darum haben wir einige Funktionen in die Numeric Library implementiert, die derartige Probleme mit iterativen Methoden lösen.
 
 
FN Zero#(&FN Fun#(0,0),I,X0#,X1#)
Fun#(I,X#) Dies ist die Funktion, deren Nullstelle berechnet werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X0# Linke Grenze des Intervalls, in dem die Nullstelle liegt.
X1# Rechte Grenze des Intervalls, in dem die Nullstelle liegt.
Diese Funktion sucht eine Nullstelle einer beliebigen Funktion mit Hilfe der Sekantenmethode. Dabei werden die Funktionswerte von X0# und X1# als Startwerte benutzt. Wenn nach 100 Schritten keine Nullstelle gefunden wurde, erzeugt FN Zero# eine Fehlermeldung und bricht ab. Als Funktionswert wird in diesem Fall NAN (Not A Number) zurückgegeben.

Beispiel:
Es soll die Stelle berechnet werden, an der SIN(3*X#)=COS(5*X#) ist. Diese Gleichung formen wir zuerst in die Form SIN(3*X#)-COS(5*X#)=0 um. Als Grenzen des Suchintervalls nehmen wir X0#=0 und X1#=PI. Damit sieht das Programm wie folgt aus:

Num_Init
PRINT FN Zero#(&FN Trig#(0,0),0,0,PI)
INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Trig#(I,X#)=SIN(3*X#)-COS(5*X#)

Wenn Sie alles richtig gemacht haben, müssten Sie als Wert 2.356... erhalten.
 
 
FN Minimum#(&FN Fun#(0,0),I,E#,S#)
Fun#(I,X#) Dies ist die Funktion, deren Minimum berechnet werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
E# Schätzwert für die Lage des Minimums.
S# Schrittweite, mit der die Suche beginnen soll. Ein typischer Wert ist z.B. S#=5.
Diese Funktion sucht ein Minimum nach der Methode der parabolischen inversen Interpolation. Dabei wird mit dem übergebenen Schätzwert E# und der Schrittweite S# begonnen. Wenn nach 200 Schritten kein Minimum gefunden wurde, erzeugt FN Minimum# eine Fehlermeldung und bricht ab. Als Funktionswert wird in diesem Fall NAN (Not A Number) zurückgegeben.

Beispiel:
Es soll das Minimum der Funktion X#^4-10*X#^2+5*X# in der Nähe von 2 gefunden werden:

Num_Init
PRINT FN Minimum#(&FN Poly#(0,0),0,2,1)
INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Poly#(I,X#)=X#^4-10*X#^2+5*X#

Wenn Sie alles richtig gemacht haben, müssten Sie als Wert 2.098... erhalten.
 
 
FN Maximum#(&FN Fun#(0,0),I,E#,S#)
Fun#(I,X#) Dies ist die Funktion, deren Maximum berechnet werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
E# Schätzwert für die Lage des Maximums.
S# Schrittweite, mit der die Suche beginnen soll. Ein typischer Wert ist z.B. S#=5.
Diese Funktion sucht ein Maximum nach der Methode der parabolischen inversen Interpolation. Dabei wird mit dem übergebenen Schätzwert E# und der Schrittweite S# begonnen. Wenn nach 200 Schritten kein Maximum gefunden wurde, erzeugt FN Maximum# eine Fehlermeldung und bricht ab. Als Funktionswert wird in diesem Fall NAN (Not A Number) zurückgegeben.
 
 





2.4 Berechnung von Eigenwerten und Lösen von Gleichungssystemen.
 
Gleichungen und Gleichungssysteme können in den meisten Fällen nur schwer analytisch gelöst werden. Entweder ist es überhaupt nicht möglich, eine analytische Lösung zu finden, oder die auftretenden Gleichungen sind sehr komplex und nur noch schwer zu handhaben. Darum haben wir einige Funktionen in die Numeric Library implementiert, die derartige Gleichungen zumindest näherungsweise lösen.
 
 
Eigen &M#(,),N,&E#()
M#(0:N,0:N) In diesem Feld müssen die Elemente der Matrix übergeben werden. Nach der Rückkehr enthält dieses Feld die Eigenvektoren. Wenn Sie Ihre Matrix später im Programm noch benötigen, müssen Sie diese zuvor mit dem BASIC-Befehl MAT in ein anderes Feld kopieren.
N Höchster Index der Matrix (Dimension minus eins).
E#(0:N) Nach der Rückkehr enthält dieses Feld die Eigenwerte.
Diese Prozedur berechnet die Eigenwerte einer reellen symmetrischen Matrix. Wenn die Matrix nicht symmetrisch ist oder der Algorithmus nicht konvergiert, erhalten Sie eine Fehlermeldung.


Beispiel:
Es sollen die Eigenwerte der Matrix berechnet werden, deren Elemente in der DATA Zeile angegeben sind:

Num_Init
-Matrix
DATA 2,-1,2,-1,2,-2,2,-2,5
DIM M#(2,2),E#(2)
RESTORE Matrix
PRINT "Matrix:"
FOR J=0 TO 2
 FOR I=0 TO 2
  READ M#(J,I):PRINT M#(J,I),
 NEXT I:PRINT
 NEXT J
Eigen &M#(,),2,&E#()
PRINT:PRINT "Eigenwerte:"
FOR I=0 TO 2:PRINT E#(I):NEXT I
PRINT
INPUT "Ende mit [Return]";Dummy
Num_Exit
END
 
 

Gauss &M#(,),&C#(),N,&X#(),R D#
M#(0:N,0:N) In diesem Feld müssen die Koeffizienten des Gleichungssystems in Form einer Matrix übergeben werden. Wenn Sie Ihre Matrix später im Programm noch benötigen, müssen Sie diese zuvor mit dem BASIC-Befehl MAT in ein anderes Feld kopieren.
C#(0:N) Dieses Feld muß die konstanten Werte auf der rechten Seite des Gleichungssystems enthalten. Wenn Sie diese Feld später im Programm noch benötigen, müssen Sie es zuvor mit dem BASIC-Befehl MAT in ein anderes Feld kopieren.
N Höchster Index der Felder (Anzahl der Gleichungen minus eins).
X#(0:N) Nach der Rückkehr enthält dieses Feld den Lösungsvektor.
D# In dieser Variablen wird die Determinante zurückgegeben.
Diese Prozedur berechnet den Lösungsvektor X#(), der das Gleichungssystem A#(,)*X#()=C#() erfüllt. Dabei wird das Verfahren von Gauss mit Pivotisierung verwendet, um den rechnerbedingten Rundungsfehler gering zu halten.

Beispiel:
Es soll das Gleichungssystem

3*X0#+4*X1#=11
2*X0#-X1#=0


gelöst werden:
Num_Init
-Equation
DATA 3,4,11,2,-1,0
DIM M#(1,1),C#(1),X#(1)
RESTORE Equation
FOR J=0 TO 1
 FOR I=0 TO 1
  READ M#(J,I)
 NEXT I
 READ C#(J)
NEXT J
Gauss &M#(,),&C#(),1,&X#(),D#
PRINT:PRINT "Die Lösung lautet:"
PRINT "X0#=";X#(0)
PRINT "X1#=";X#(1)
PRINT "Determinante D#=";D#
PRINT
INPUT "Ende mit [Return]";Dummy
Num_Exit
END
 
 

Gauss_Seidel &M#(,),&C#(),N,&X#()
M#(0:N,0:N) In diesem Feld müssen die Koeffizienten des Gleichungssystems in Form einer Matrix übergeben werden.
C#(0:N) Dieses Feld muß die konstanten Werte auf der rechten Seite des Gleichungssystems enthalten.
N Höchster Index der Felder (Anzahl der Gleichungen minus eins).
X#(0:N) In diesem Feld müssen Sie eine genäherte Lösung übergeben. Nach der Rückkehr enthält dieses Feld den verbesserten Lösungsvektor.
Diese Prozedur berechnet den Lösungsvektor X#(), der das Gleichungssystem A#(,)*X#()=C#() erfüllt. Bei sehr großen Gleichungssystemen (100 oder mehr Variablen) benötigt der Gauss-Algorithmus eine zu lange Rechenzeit und der Rundungsfehler wird zu groß. In diesen Fällen ist es besser, ein Iterationsverfahren wie das von Gauss und Seidel zu verwenden, besonders wenn man schon ungefähre Werte für den Lösungsvektor kennt.
Natürlich ist es auch möglich, dieses Verfahren dem Gauss-Algorithmus nachzuschalten und dadurch den Rundungsfehler zu verringern.

Achtung: Da das Gauss-Seidel-Verfahren iterativ arbeitet, kann es sein, daß es nicht konvergiert. Dieser Effekt tritt besonders dann auf, wenn kein genäherter Lösungsvektor übergeben wird. Wenn der Algorithmus nicht konvergiert, erhalten Sie eine Fehlermeldung.
 
 

Band &M#(,),N,K,&X#()
M#(0:N,0:2*K+1) In diesem Feld müssen die Koeffizienten des Gleichungssystems übergeben werden. Bitte beachten Sie, daß in der letzten Spalte die konstanten Werte von der rechten Seite des Gleichungssystems übergeben werden müssen. Wenn Sie dieses Feld später im Programm noch benötigen, müssen Sie es zuvor mit dem BASIC-Befehl MAT in ein anderes Feld kopieren.
N Höchster Index der Felder (Anzahl der Gleichungen minus eins).
K Abstand des am weitesten von der Hauptdiagonalen entfernten Elements, das ungleich null ist.
X#(0:N) Nach der Rückkehr enthält dieses Feld den Lösungsvektor.
Diese Prozedur berechnet den Lösungsvektor X#(), der das Gleichungssystem A#(0:N,0:N)*X#(0:N)=A#(0:N,N+1) erfüllt.
In der Praxis treten häufig Gleichungssysteme auf, bei denen die Koeffizienten-Matrix nur in der Nähe der Hauptdiagonalen von null verschiedene Elemente enthält. Solche Matrizen bezeichnet man als Bandmatrizen mit der Breite
K. Ein solches Gleichungssystem läßt sich mit der Prozedur Band viel schneller lösen als mit dem Gauss-Verfahren, da überflüssige Rechenoperationen mit den Nullen der Matrix vermieden werden. Der Vorteil ist umso größer, je größer die Matrix und je schmaler das Band im Verhältnis zur Dimension der Matrix ist.

Beispiel:
Es soll das folgende Gleichungssystem gelost werden:

2*X0#+9*X1#=1
2*X0#-3*X1#+12*x2#=2
8*X1#+3*X2#-5*X3#=3
6*X2#+4*X3#+X4#=4
X3#-2*X4#=-5

Wie man leicht erkennt, sind nur die Elemente auf der Hauptdiagonalen, sowie Elemente rechts und links davon ungleich null. Die Koeffizienten-Matrix kann daher auf eine Band-Matrix mit 3 Spalten reduziert werden. In der vierten Spalte werden dann noch die konstanten Elemente von der rechten Seite des Gleichungssystems eingetragen. Daher ergibt sich mit
K=1 folgendes Programm:

Num_Init
-Band
DATA 0,2,9,1
DATA 2,-3,12,2
DATA 8,3,-5,3
DATA 6,4,1,4
DATA 1,-2,0,-5
DIM M#(4,3),X#(4)
RESTORE Band
FOR J=0 TO 4
 FOR I=0 TO 3
  READ M#(J,I)
 NEXT I
NEXT J
Band &M#(,),4,1,&X#()
PRINT:PRINT "Die Lösung lautet:"
FOR I=0 TO 4
 PRINT "X"+MID$(STR$(I),2)+"#=";X#(I)
NEXT I
PRINT
INPUT "Ende mit [Return]";Dummy
Num_Exit
END
 
 

Newton_Sys &FN Fun#(0),&X#(),N
Fun#(I) Durch diese Funktion definieren Sie das Gleichungssystem, das gelöst werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, der die Nummer der Gleichung identifiziert. Die Gleichung muß so umgeformt sein, daß auf der rechten Seite eine 0 steht. Ihre Funktion muß dann den Funktionswert der linken Seite zurückliefern.
X#(0:N) In diesem Feld sollten Sie eine genäherte Lösung übergeben, falls diese bekannt ist. Nach der Rückkehr enthält dieses Feld den Lösungsvektor.
N Höchster Index des Feldes (Anzahl der Gleichungen minus eins).
Diese Prozedur löst ein nichtlineares Gleichungssystem nach dem Newton-Verfahren. Wenn nach 100 Schritten keine Lösung gefunden wurde, erzeugt Newton_Sys eine Fehlermeldung und bricht ab.

Beispiel:
Es soll das nichtlineare Gleichungssystem

X0#*X1#-5*X1#-5=0
X0#-2+X1#*COS(X2#)=0
SIN(X0#+3*X2#)=0


gelöst werden:

Num_Init
DIM X#(2)
Newton_Sys &FN Equation#(0),&X#(),2
PRINT:PRINT "Die Lösung lautet:"
FOR I=0 TO 2
 PRINT "X"+MID$(STR$(I),2)+"#=";X#(I)
NEXT I
PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Equation#(I):LOCAL X#
 SELECT I
  CASE 0:X#=X#(0)*X#(1)-5*X#(1)-5
  CASE 1:X#=X#(0)-2+X#(1)*COS(X#(2))
  CASE 2:X#=SIN(X#(0)+3*X#(2))
 END_SELECT
RETURN X#
 
 




2.5 Interpolation
 
Oft kommt es vor, daß von einer benötigten Funktion nur tabellierte Werte bekannt sind. In solchen Fällen bietet es sich an, die Zwischenwerte durch Interpolation mit einer bekannten Funktion zu ermitteln. Die Numeric Library bietet Ihnen diverse Prozeduren und Funktionen, die für fast alle Interpolationsprobleme geeignet sind.
Bitte denken Sie aber immer daran, daß eine Interpolation nicht notwendigerweise die richtigen Zwischenwerte liefert. Z. B. können nicht differenzierbare oder gar unstetige Funktionen dadurch überhaupt nicht korrekt dargestellt werden.
 
 
Lagrange_Int &X#(,),N,&C#()
X#(0:N,0:1) In diesem Feld übergeben Sie die Stützstellen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
N Höchster Index der Felder (Anzahl der Stützpunkte minus eins).
C#(0:N) Nach der Rückkehr enthält dieses Feld die Koeffizienten des Interpolationspolynoms. Zur Berechnung der Interpolationsfunktion muß dieses Feld an FN Lagrange# übergeben werden.
Die Lagrange-Interpolation ist ein sehr gutes Verfahren, weil es im Gegensatz zu anderen auch nicht äquidistante Werte zuläßt. Bei sehr vielen Stützpunkten beginnt das Interpoaltionspolynom allerdings an den Intervallenden zu schwingen, wodurch die Qualität der Interpolation verschlechtert wird.

Im Ordner DEMO finden Sie das Programm "Interpolation.BAS", mit dem Sie selbstgewählte Stützpunkte mit den verschiedenen Verfahren interpolieren können.
 
 
FN Lagrange# X#,&X#(,),N,&C#()
X# An dieser Stelle wird der Funktionswert der Lagrange-Interpolation berechnet. X# muß natürlich innerhalb des Intervalls liegen, in dem die Funktionswerte bekannt sind, sonst erhalten Sie unsinnige Ergebnisse.
X#(0:N,0:1) Hier müssen Sie das gleiche Feld mit den Stützstellen wie bei der Prozedur Lagrange_Int übergeben.
N Höchster Index der Felder (Anzahl der Stützpunkte minus eins).
C#(0:N) Dieses Feld muß die Koeffizienten des Interpolationspolynoms enthalten. Es ist genau das Feld, das Sie von Lagrange_Int zurückerhalten haben.
Diese Funktion dient dazu, die Werte des Interpolationspolynoms zu berechnen. In der Praxis rufen Sie also zuerst die Prozedur Lagrange_Int auf, um die Koeffizienten des Interpolationspolynoms zu ermitteln und berechnen dann mit FN Lagrange# die Zwischenwerte.
 
 

Spline_Int &X#(,),N,&C#(,)
X#(0:N,0:1) In diesem Feld übergeben Sie die Stützstellen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
N Höchster Index der Felder (Anzahl der Stützpunkte minus eins).
C#(0:N,0:3) Nach der Rückkehr enthält dieses Feld die Koeffizienten der Interpolationspolynome (4 pro Stützstelle). Zur Berechnung der Interpolationsfunktion muß dieses Feld an FN Spline# übergeben werden.
Diese Prozedur berechnet eine kubische Spline-Interpolation. Dadurch, daß am Ende jedes Interpolationsintervalls auch die Ableitung berücksichtigt wird, liefert die Spline-Interpolation einen sehr glatten Interpolationsverlauf. Auch die Spline-Interpolation erlaubt unterschiedliche Abstände zwischen den Stützpunkten.
Im Unterschied zur Lagrange-Interpolation wird durch die Verwendung von Polynomen dritten Grades, die an den Stützpunkten zusammengesetzt werden, das starke Schwingen der Polynome bei einer hohen Anzahl von Stützstellen vermieden.
Die Spline-Interpolation hat die Eigenschaft, daß die Spannung in der Kurve minimal wird, die Kurve also einem durch die Stützpunkte gelegten elastischen Lineal (engl. spline) entspricht.

Im Ordner DEMO finden Sie das Programm "Interpolation.BAS", mit dem Sie selbstgewählte Stützpunkte mit den verschiedenen Verfahren interpolieren können.
 
 
FN Spline# X#,&X#(,),N,&C#()
X# An dieser Stelle wird der Funktionswert der Spline-Interpolation berechnet. X# muß natürlich innerhalb des Intervalls liegen, in dem die Funktionswerte bekannt sind, sonst erhalten Sie unter Umständen falsche Ergebnisse.
X#(0:N,0:1) Hier müssen Sie das gleiche Feld mit den Stützstellen wie bei der Prozedur Spline_Int übergeben.
N Höchster Index der Felder (Anzahl der Stützpunkte minus eins).
C#(0:N,0:3) Dieses Feld muß die Koeffizienten des Interpolationspolynoms enthalten (4 pro Stützstelle). Es ist genau das Feld, das Sie von Spline_Int zurückerhalten haben.
Diese Funktion dient dazu, die Werte der Interpolation zu berechnen. In der Praxis rufen Sie also zuerst die Prozedur Spline_Int auf, um die Koeffizienten der Interpolationspolynome zu ermitteln und berechnen dann mit FN Spline# die Zwischenwerte.


 
 

Gauss_Int &FN Fun#(0,0),I,&X#(,),N,&P#(),M
Fun#(I,X#) Dies ist eine Funktion mit M frei wählbaren Parametern, die im Feld P#() stehen. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X#(0:N,0:2) In diesem Feld übergeben Sie die Stützstellen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
In
X#(0:N,2) können Sie noch eine Fehlertoleranz für die einzelnen Stützpunkte angeben. Im allgemeinen kann X#(0:N,2)=1 gesetzt werden.
N Höchster Index des Feldes X#(,) (Anzahl der Stützpunkte minus eins).
P#(0:M) Nach der Rückkehr enthält dieses Feld die Parameter Ihrer Funktion.
M Höchster Index des Feldes P#() (Anzahl der frei wählbaren Parameter der Funktion FN Fun minus eins).
Diese Prozedur passt eine von Ihnen definierte Funktion mit M+1 frei wählbaren Parametern an N+1 Stützpunkte an. Dabei wird die Methode der kleinsten Fehlerquadrate nach Gauss verwendet, d.h. die Summe aus den Quadraten der Fehler an den Stützstellen wird minimiert. Auch die Gauss-Interpolation erlaubt unterschiedliche Abstände zwischen den Stützpunkten.

Hinweis: Wenn Ihre Funktion weniger frei wählbare Parameter hat, als Stützpunkte übergeben werden, läuft die Funktion im allgemeinen nicht durch die Stützpunkte.

Im Ordner DEMO finden Sie das Programm "Interpolation.BAS", mit dem Sie selbstgewählte Stützpunkte mit den verschiedenen Verfahren interpolieren können.

 
 

Fourier_Int &X#(,),N,X1#,X2#,&C#(),M
X#(0:N,0:1) In diesem Feld übergeben Sie die Stützstellen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
N Höchster Index des Feldes X#(,) (Anzahl der Stützpunkte minus eins).
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
C#(0:M+M) Nach der Rückkehr enthält dieses Feld die Koeffizienten der trigonometrischen Interpolationsfunktion. Zur Berechnung der Funktionswerte muß dieses Feld an FN Fourier# übergeben werden.
M Die Anzahl der trigonometrischen Terme.
Diese Prozedur berechnet die Koeffizienten einer Interpolationsfunktion, die aus Kosinus- und Sinus-Funktionen besteht. Dieser Interpolationstyp eignet sich daher besonders zur Interpolation von Werten, die sich periodisch wiederholen, wie zum Beispiel Töne.
Die Interpolationsfunktion hat dabei die Form:

C#(0)+C#(1)*COS(X#)+C#(2)*COS(2*X#)+ ... +C#(M)*COS(M*X#)+ ... +C#(M+1)*SIN(X#)+C#(M+2)*SIN(2*X#)+C#(M+M)*SIN(M*X#)

mit X#=2*PI*X#/(X2#-X1#)-PI

Achtung: Die Stützstellen müssen äquidistant sein. Bei periodischen Funktionen muß sich das Intervall über mindestens eine oder mehrere ganze Perioden erstrecken.

Im Ordner DEMO finden Sie das Programm "FourierInterpolation.BAS", das die Verwendung dieser Interpolationsmethode demonstriert.
 
 

FN Fourier# X#,X1#,X2#,&C#(),M
X# An dieser Stelle wird der Funktionswert der Fourier-Interpolation berechnet. Wenn eine nichtperiodische Punktreihe interpoliert werden soll, muß X# natürlich innerhalb des Intervalls liegen, in dem die Funktionswerte bekannt sind, sonst erhalten Sie unter Umständen falsche Ergebnisse.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
C#(0:M+M) Dieses Feld muß die Koeffizienten der trigonometrischen Interpolationsfunktion enthalten. Es ist genau das Feld, das Sie von Fourier_Int zurückerhalten haben.
M Die Anzahl der trigonometrischen Terme.
Diese Funktion dient dazu, die Werte der Interpolation zu berechnen. In der Praxis rufen Sie also zuerst die Prozedur Fourier_Int auf, um die Koeffizienten der trigonometrischen Terme zu ermitteln und berechnen dann mit FN Fourier# die Zwischenwerte. Bei periodischen Meßreihen gilt die Interpolation auch außerhalb des interpolierten Intervalls.
 
 

Fourier_Int_Sin &X#(,),N,X1#,X2#,&C#(),M
X#(0:N,0:1) In diesem Feld übergeben Sie die Stützstellen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
N Höchster Index des Feldes X#(,) (Anzahl der Stützpunkte minus eins).
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
C#(0:M+M) Nach der Rückkehr enthält dieses Feld die Koeffizienten der trigonometrischen Interpolationsfunktion. Zur Berechnung der Funktionswerte muß dieses Feld an FN Fourier_Sin# übergeben werden.
M Die Anzahl der Sinus-Funktionen.
Diese Prozedur berechnet die Koeffizienten einer Interpolationsfunktion, die nur aus Sinus-Funktionen und einem konstanten Term C#(0) besteht. Dieser Interpolationstyp eignet sich daher besonders zur Interpolation von Werten, die sich periodisch wiederholen, wie zum Beispiel Töne.
Die Interpolationsfunktion hat dabei die Form:

C#(0)+C#(1)*SIN(X#+C#(M+1))+C#(2)*SIN(2*X#+C#(M+2))+ ... +C#(M)*SIN(M*X#+C#(M+M)

mit X#=2*PI*X#/(X2#-X1#)-PI

Im Unterschied zu Fourier_Int enthält die Interpolationsfunktion nur die halbe Anzahl an trigonometrischen Funktionen; die Berechnug dürfte im allgemeinen also schneller erfolgen.

Achtung: Die Stützstellen müssen äquidistant sein. Bei periodischen Funktionen muß sich das Intervall über mindestens eine oder mehrere ganze Perioden erstrecken.

Im Ordner DEMO finden Sie das Programm "FourierInterpolation.BAS", das die Verwendung dieser Interpolationsmethode demonstriert.
 
 

FN Fourier_Sin# X#,X1#,X2#,&C#(),M
X# An dieser Stelle wird der Funktionswert der Fourier-Sinus-Interpolation berechnet. Wenn eine nichtperiodische Punktreihe interpoliert werden soll, muß X# natürlich innerhalb des Intervalls liegen, in dem die Funktionswerte bekannt sind, sonst erhalten Sie unter Umständen falsche Ergebnisse.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
C#(0:M+M) Dieses Feld muß die Koeffizienten der trigonometrischen Interpolationsfunktion enthalten. Es ist genau das Feld, das Sie von Fourier_Int_Sin zurückerhalten haben.
M Die Anzahl der Sinus-Funktionen.
Diese Funktion dient dazu, die Werte der Interpolation zu berechnen. In der Praxis rufen Sie also zuerst die Prozedur Fourier_Int_Sin auf, um die Koeffizienten und Phasen der Sinus-Funktionen zu ermitteln und berechnen dann mit FN Fourier_Sin# die Zwischenwerte. Bei periodischen Meßreihen gilt die Interpolation auch außerhalb des interpolierten Intervalls.
 
 

Fft &R#(),&I#(),N,Flag
R#(0:N) Realteile der komplexen Werte. Nach der Rückkehr enthält dieses Feld die Realteile der transformierten Werte. Wenn Sie die ursprünglichen Werte später im Programm noch benötigen, müssen Sie diese zuvor mit dem BASIC-Befehl MAT in ein anderes Feld kopieren.
I#(0:N) Imaginärteile der komplexen Werte. Nach der Rückkehr enthält dieses Feld die Realteile der transformierten Werte. Wenn Sie die ursprünglichen Werte später im Programm noch benötigen, müssen Sie diese zuvor mit dem BASIC-Befehl MAT in ein anderes Feld kopieren.
N Anzahl der Werte. Dies muß eine gerade Zahl sein.

Hinweis: Die Prozedur arbeitet besonders schnell, wenn N eine Zweierpotenz ist, aber auch Zahlen, die sich gut in Faktoren von 2,3,4,5 zerlegen lassen, ergeben kurze Rechenzeiten.
Flag Flag = 1 : Die normale Fourier-Transformation wird durchgeführt.
Flag = -1 : Die inverse Fourier-Transformation wird durchgeführt.
Diese Prozedur führt eine diskrete Fourier-Transformation durch. Dabei werden im Prinzip folgende Summen gebildet.

Normale Fourier-Transformation:

R#(I)=SUM[J=0,N-1]R#(J)*COS(2*PI*I/N)-I#(J)*SIN(2*PI*I/N)
I#(I)=SUM[J=0,N-1]I#(J)*COS(2*PI*I/N)+R#(J)*SIN(2*PI*I/N)
R#(N)=R#(0):I#(N)=I#(0)


Inverse Fourier-Transformation:

R#(I)=(SUM[J=0,N-1]R#(J)*COS(2*PI*I/N)+I#(J)*SIN(2*PI*I/N))/N
I#(I)=(SUM[J=0,N-1]I#(J)*COS(2*PI*I/N)-R#(J)*SIN(2*PI*I/N))/N
R#(N)=R#(0):I#(N)=I#(0)


Die Prozedur arbeitet allerdings nach der Methode der schnellen Fourier-Transformation, wodurch die Berechnungszeit wesentlich kürzer ausfällt, als wenn man die obigen Summen direkt ausrechnen würde.

Im Ordner DEMO finden Sie das Programm "FFT.BAS", das die Verwendung dieser Prozedur demonstriert.
 
 


2.6 Approximation von Funktionen durch Funktionensysteme
 
In vielen Fällen steht man vor dem Problem, daß eine Funktion zwar analytisch bekannt, die Auswertung jedoch schwierig ist. Dann bietet es sich an, die Funktion durch andere, einfachere Funktionen anzunähern. Die gängigsten Verfahren für diesem Zweck sind in der Numeric Library enthalten.
 
 
Cheby_App &FN Fun#(0,0),I,X1#,X2#,&C#(),N
Fun#(I,X#) Dies ist die Funktion, die approximiert werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
C#(0:N) Nach der Rückkehr enthält dieses Feld die Koeffizienten für die Tschebyscheff-Approximation. Zur Berechnung der Funktionswerte muß dieses Feld an FN Cheby# übergeben werden.
N Höchster Index des Feldes C#() (Anzahl der Tschebyscheff-Polynome minus eins).

Diese Prozedur passt eine von Ihnen definierte analytische Funktion FN Fun# an das orthogonale Funktionensystem aus Tschebyscheff-Polynomen an.
Der Vorteil der Tschebyscheff-Approximation besteht darin, daß normalerweise nur wenige Koeffizienten benötigt werden, um eine sehr gute Approximation zu erreichen und das es sehr effiziente Verfahren gibt, die Tschebyscheff-Polynome schnell zu berechnen. Dazu können Sie die Funktion
FN Cheby# verwenden.

Hinweis: Da, anders als bei der weiter unten besprochenen Prozedur Mean_App, zur Berechnung der Koeffizienten keine Integrale ausgewertet werden müssen, sind für die Anzahl der Polynome auch größere Zahlen möglich (N > 100). Zwar dauert dann die Berechnung schon etwas länger, man erhält aber recht gute Ergebnisse.

Im Ordner DEMO finden Sie das Programm "Approximation.BAS", das die Verwendung dieser Prozedur zeigt.
 
 
FN Cheby# X#,X1#,X2#,&C#(),N
X# An dieser Stelle wird der Funktionswert der Tschebyscheff-Approximation berechnet. X# muß natürlich innerhalb des Intervalls liegen, in dem die Funktionswerte bekannt sind, sonst erhalten Sie unter Umständen falsche Ergebnisse.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
C#(0:N) Dieses Feld muß die Koeffizienten der Tschebyscheff-Approximationsfunktion enthalten. Es ist genau das Feld, das Sie von Cheby_App zurückerhalten haben.
N Höchster Index des Feldes C#() (Anzahl der Tschebyscheff-Polynome minus eins).
Diese Funktion dient dazu, die Werte der Approximation zu berechnen. In der Praxis rufen Sie also zuerst die Prozedur Cheby_App auf, um die Koeffizienten der Tschebyscheff-Polynome zu ermitteln und berechnen dann mit FN Cheby# die Approximationswerte.
Verwechseln Sie diese Funktion nicht mit der oben besprochenen Funktion
FN Cheby#(I,X#), die nur den Funktionswert eines ganz bestimmten Tschebyscheff-Polynoms berechnet.

Im Ordner DEMO finden Sie das Programm "Approximation.BAS", das die Verwendung dieser Prozedur zeigt.
 
 
Mean_App &FN Fun#(0,0),I,X1#,X2#,&FN Sys#(0,0),&C#(),N
Fun#(I,X#) Dies ist die Funktion, die approximiert werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
Sys#(I,X#) Dies ist das Funktionen-System, mit dem FN Fun# approximiert werden soll. Es muß aus N+1 Elementen bestehen und die Koeffizienten C#(0:N) verwenden.
Dieses Funktionensystem müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
C#(0:N) Nach der Rückkehr enthält dieses Feld die Koeffizienten für Ihr Funktionen-System.
N Höchster Index des Feldes C#() (Anzahl der von Ihnen definierten Funktionen minus eins).

Diese Prozedur passt eine von Ihnen definierte analytische Funktion FN Fun# an ein ebenfalls von Ihnen definiertes Funktionen-System FN Sys# dergestalt an, daß das Fehlerquadrat-Integral minimal wird. Die angepasste Funktion berechnet sich aus:

C#(0)*FN Sys#(0,X#)+C#(1)*FN Sys#(1,X#)+ ...
+C#(N)*FN Sys#(N,X#)

Hinweis: Da zur Berechnung der Koeffizienten sehr viele Integrale ausgewertet werden müssen, kann die Rechenzeit bei größeren Funktionensystemen (N >> 10) schon recht lang werden. Außerdem können Zahlenbereichsüberschreitungen auftreten, die das Ergebnis völlig verfälschen. Sie sollten darum zumindest in der Testphase das Steuerwort COMPILER "DEBUG ON" setzen, auch wenn die Rechengeschwindigkeit dadurch etwas abnimmt.

Beispiel:
Es soll die Eulersche Funktion (EXP(X#)) durch eine Potenzreihe approximiert werden:

Num_Init
N=5
DIM C#(N)
Mean_App &FN Euler#(0,0),0,-1,1,&FN Potenz_Sys#(0,0),&C#(),N
PRINT:PRINT " X#";TAB(7);"Ausgangsfunktion";
PRINT TAB(27);"Approximation"
FOR X#=-1 TO 1 STEP 0.1
 PRINT USING "##.##";X#;TAB(6);USING "";EXP(X#);
 PRINT TAB(25);FN Potenz_App#(X#)
NEXT X#
PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Euler#(I,X#)=EXP(X#)
DEF FN Potenz_Sys#(I,X#)=X#^I
DEF FN Potenz_App#(X#):LOCAL I,Y#=0
 FOR I=0 TO N
  Y#+=C#(I)*X#^I
 NEXT I
RETURN Y#


Wie Sie mit diesem Programm leicht ausprobieren können, ergibt N=5 bessere Approximationswerte als N=10. Bei N=15 ergeben sich schon flasche Werte, weil Rundungsfehler sich bei der Berechnung der Koeffizienten akkumulieren und bei N=20 kommt es bereits zu einem Underflow.


Im Ordner DEMO finden Sie das Programm "Approximation.BAS", das ebenfalls die Verwendung dieser Prozedur zeigt.
 
 





2.7 Numerische Differentiation
 
Dieses Kapitel beschäftigt sich mit der Ableitung von Funktionen. Die numerische Differentiation ist sehr empfindlich gegen Rundungsfehler. Die Ergebnisse können daher immer nur Näherungswerte sein.
 
 
FN Derivative &FN Fun#(0,0),I,X#[,D#]
Fun#(I,X#) Dies ist die Funktion, deren Ableitung berechnet werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X# Die Stelle, an der die Ableitung berechnet werden soll.
D# In diesem Parameter kann optional eine Intervallgrenze angegeben werden, bei der die Differentiation abgebrochen wird. Der Standardwert ist 1D-10. Größere Werte können bei ungenauen Funktionen von Vorteil sein, da der Algorithmus die feinen Ungenauigkeiten dann nicht "sieht".
Diese Funktion berechnet die erste Ableitung einer analytischen Funktion. Durch wiederholte Anwendung auf sich selbst können auch höhere Ableitungen erhalten werden. Dann muß man im allgemeinen allerdings schrittweise die Intervallgrenze erhöhen.

Im Ordner DEMO finden Sie das Programm "Derivative.BAS", das die Verwendung dieser Funktion zeigt.
 
 
Derivative &X#(,),N
X#(0:N,0:3) In diesem Feld übergeben Sie die Funktionswerte. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
Nach der Rückkehr enthält
X#(0:N,1) die 1. Ableitung, X#(0:N,2) die 2. Ableitung und X#(0:N,3) die 3. Ableitung and der Stelle X#(0:N,0).
N Höchster Index des Feldes X#(,) (Anzahl der Funktionswerte minus eins).
Diese Prozedur berechnet die ersten drei Ableitungen einer tabellarisch vorliegenden Funktion. Bedenken Sie bitte, daß die Qualität der Ergebnisse stark von der Qualität der Eingaben in X#(0:N,0:1) abhängt. Dafür gilt nämlich die alte Weißheit: "Garbage in, garbage out (Fehlerhafte Eingabewerte liefern fehlerhafte Ausgabewerte)".

Beispiel:
Es soll die ersten drei Ableitungen der Funktion X#^3 berechnet werden:

Num_Init
N=20:I=0
DIM X#(N,4)
FOR X#=-1 TO 1 STEP 0.1
 X#(I,0)=X#
 X#(I,1)=X#^3
 I+=1
NEXT X#
Derivative &X#(,),N
PRINT:PRINT " X#";TAB(7);"Funktionswert";
PRINT TAB(27);"1. Ableitung";TAB(47);"2. Ableitung";
PRINT TAB(67);"3. Ableitung"
I=0
FOR X#=-1 TO 1 STEP 0.1
 PRINT USING "##.##";X#;TAB(7);USING "##.#########";X#(I,1);
 PRINT TAB(27);USING "##.#########";X#(I,2);TAB(47);
 PRINT TAB(27);USING "##.#########";X#(I,3);TAB(67);X#(I,4)
 I+=1
NEXT X#
PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END

Da die Funktionswerte von einer glatten analytischen Funktion stammen, sind die Ergebnisse in diesem Fall auch recht gut, wie Sie leicht nachrechnen können.

 
 





2.8. Numerische Integration
 
Im Unterschied zur Differentiation gibt es für die Integration von Funktionen schon keine allgemeingültigen Algorithmen mehr, die es immer erlauben, den Wert des Integrals immer in Form eines geschlossenen Ausdrucks anzugeben. So können schon relativ einfach aussehende Funktionen wie SIN(X#)/X# nicht elementar integriert werden.
Um so wichtiger ist es, daß man numerische Methoden zur Verfügung hat, mit denen man solche Integrale mit guter Genauigkeit berechnen kann.
 
 
FN Newton_Cotes#(&X#(,),N)
X#(0:N,0:1) In diesem Feld übergeben Sie die Funktionswerte. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.

Achtung: Die X-Werte müssen äquidistant sein.
N Höchster Index des Feldes X#(,) (Anzahl der Funktionswerte minus eins). N muß größer oder gleich 5 sein, sonst erhalten Sie eine Fehlermeldung.
Diese Funktion berechnet das Integral der im Feld X#(,) übergebenen diskreten Funktionswerte.

Hinweis: Wenn Sie keine äquidistanten Werte haben, können Sie auch zunächst eine der zuvor besprochenen Interpolationsfunktionen verwenden, um sich äquidistante Werte zu verschaffen und diese dann an FN Newton_Cotes# übergeben. Alternativ können Sie auch an eine der weiter unten besprochenen Integrationsfunktionen verwenden, an die Sie dann die Adresse der Interpolationsfunktion übergeben.

Beispiel:
Bei einer Messung haben sich die den DATA-Zeilen aufgeführten Meßwerte ergeben. Es soll jetzt das Integral, also die Fläche unter der Meßkurve, berechnet werden. Da die X-Werte nicht äquidistant sind, wird zunächst interpoliert:

Num_Init
-Values
DATA 0,1,0.5,1.649,1,2.718,2,1.4,4,0.5,5,2
N=5
DIM X#(N,1),Y#(N,1),C#(N)
RESTORE Values
FOR I=0 TO 5:READ X#(I,0),X#(I,1):NEXT I
'Erstmal äquidistante Meßwerte durch
'Interpolation erzeugen.
Lagrange_Int &X#(,),N,&C#()
FOR I=0 TO 5
 Y#(I,0)=I
 Y#(I,1)=FN Lagrange#(I,&X#(,),N,&C#())
NEXT I
PRINT
PRINT "Das Integral hat den Wert: ";
PRINT FN Newton_Cotes#(&Y#(,),N)
PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END
 
 
FN Gauss_Leg#(&FN Fun#(0,0),I,X1#,X2#,N)
Fun#(I,X#) Dies ist die Funktion, die integriert werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
N Mit diesem Parameter können Sie die Genauigkeit beeinflussen. N kann Werte zwischen 2 und 15 annehmen. Je höher der Wert für N ist, desto genauer ist das Ergebnis, aber desto länger dauert auch die Berechnung.
Die Integration nach Gauss-Legendre ist eine der schnellsten Methoden, Integrale numerisch zu berechnen. Sie funktioniert allerdings nur gut für relativ harmlose Funktionen wie SIN(X#) oder COS(X#) auf einem begrenzten Intervall. Schon die Berechnung des Integrals von 1/X# im Intervall [1,1000] führt selbst bei N=15 zu völlig falschen Ergebnissen. Verwenden Sie daher für solche Fälle besser die weiter unten besprochene Funktion FN Romberg#.

Beispiel:
Es soll das Integral von SIN(X#) in Intervall [0,PI] berechnet werden:

Num_Init
PRINT
PRINT "Das Integral hat den Wert: ";
PRINT FN Gauss_Leg#(&FN Sinus#(0,0),0,0,PI,5)
PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Sinus#(I,X#)=SIN(X#)

Das Ergebnis von 2.00000011020611 liegt ziemlich dicht bei dem theoretischen Ergebnis von 2, obwohl für N nur 5 gewählt wurde. Probieren Sie die Berechnung mal selbst mit anderen Werten für N und anderen Intervallgrenzen.

Im Ordner DEMO finden Sie das Programm Integral.BAS, das verschiedene Integrationsverfahren miteinander vergleicht.
 
 
FN Romberg#(&FN Fun#(0,0),I,X1#,X2#)
Fun#(I,X#) Dies ist die Funktion, die integriert werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
Diese Funktion berechnet das Integral nach der Methode von Romberg. Die Rechengeschwindigkeit ist deutlich geringer als bei FN Gauss_Leg#, dafür sind die Ergebnisse aber auch bei problematischen Funktionen recht genau. Dies wird dadurch erreicht, daß ein rekursiver Algorithmus die Funktion selbst wieder aufruft, bis ein relativer Fehler von 1D-10 unterschritten ist. Dazu wird Speicher auf dem BASIC-Stack benötigt, den Sie gegebenenfalls mit COMPILER "STACK X" vergrößern müssen.
Trotz der guten Konvergenz von
FN Romberg# können Sie natürlich nicht ohne weiteres über Singularitäten hinwegintegrieren. Sie müssen dafür sorgen, daß Ihre Funktion im gesamten Integrationsintervall nur Werte innerhalb des darstellbaren Zahlenbereichs abzüglich eines Sicherheitsintervalls annimmt.
Falls die Berechnung nicht konvergiert, erhalten Sie eine Fehlermeldung.

Bitte beachten Sie auch das Beispielprogramm Romberg.BAS im Ordner DEMO. Es berechnet den Integral-Sinus, den Integral-Cosinus, die Integral-Exponentialfunktion und den Integral-Logarithmus. Alle 4 Integrale können nicht in geschlossener Form angegeben werden und die Integranden der letzten 3 haben obendrein Singularitäten, die durch endliche Werte ersetzt werden müssen.
Bei der zweiten und dritten Funktion ist als Integrationsgrenze 100 bzw. -100 angenommen. Theoretisch müßte hier natürlich unendlich bzw. minus unendlich stehen. Da der Integrand jedoch gegen 0 strebt, bekommt man nur Rundungsfehler, die das gesamte Ergebnis verfälschen, wenn man das Intervall zu groß wählt.
 
 
FN Romberg_2#(&FN Fun1#(0,0),I1,&FN Fun2#(0,0),I2,X1#,X2#)
Fun1#(I1,X#)
Fun2#(I2,X#)
Dies sind die beiden Funktionen, deren Produkt integriert werden soll. Diese Funktionen müssen Sie natürlich selbst definieren und dann die Adressen übergeben.
In
I1 und I2 wird Ihren Funktionen ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswerte Ihre Funktionen zurückliefern muß.
I1,I2 Ein Index, der an FN Fun1# bzw. FN Fun2# übergeben wird.
X1# Linke Intervallgrenze.
X2# Rechte Intervallgrenze.
Bei vielen Aufgaben, wie z.B. der Entwicklung von Funktionen nach Funktionensystemen, stellt sich das Problem, daß man nicht einzelne Funktionen, sondern die Produkte von unterschiedlichen Funktionen integrieren muß. Dafür haben wir FN Romberg_2# in die Numeric Library mit aufgenommen. Es wird nach der Methode von Romberg integriert. Die Ergebnisse sind auch bei problematischen Funktionen recht genau. Dies wird dadurch erreicht, daß ein rekursiver Algorithmus die Funktion selbst wieder aufruft, bis ein relativer Fehler von 1D-10 unterschritten ist. Dazu wird Speicher auf dem BASIC-Stack benötigt, den Sie gegebenenfalls mit COMPILER "STACK X" vergrößern müssen.
Trotz der guten Konvergenz von
FN Romberg_2# können Sie natürlich nicht ohne weiteres über Singularitäten hinwegintegrieren. Sie müssen dafür sorgen, daß Ihre Funktionen im gesamten Integrationsintervall nur Werte innerhalb des darstellbaren Zahlenbereichs abzüglich eines Sicherheitsintervalls annehmen.
Falls die Berechnung nicht konvergiert, erhalten Sie eine Fehlermeldung.
 
 




2.9 Lösung von Differentialgleichungen
 
Das Lösen von Differentialgleichungen und Differentialgleichungssystemen ist ein sehr wichtiger Bereich in der numerischen Mathematik, da die Lösungsfunktionen in vielen Fällen nicht anders gefunden und häufig auch nicht in geschlossener Form dargestellt werden können. Wir haben daher die gängigsten Verfahren in die Numeric Library integriert.


 
 
Runge_Kutta &FN Fun#(0,0,0),I,Xs#,Xe#,R Y#,R S#,E#
Fun#(I,X#,Y#) Diese Funktion definiert die Differentialgleichung über die Beziehung Y#'=FN Fun#(I,X#,Y#). Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
I Ein Index, der an FN Fun# übergeben wird.
Xs# Anfangswert für X#.
Xe# Endwert für X#.
Y# Vor dem Aufruf müssen Sie hier den Anfangswert für Y# übergeben. Nach der Rückkehr erhalten Sie in dieser Variablen den Endwert von Y# zurück.
S# In diesem Parameter müssen Sie eine Empfehlung für die erste Schrittweite übergeben (S#>=0.01). Nach der Rückkehr enthält diese Variable die zuletzt verwendete Schrittweite.
E# Hiermit legen Sie fest, wie groß der lokale Fehler maximal sein darf (1D-14 <= E# <= 1D-4). Unter dem lokalen Fehler versteht man den verfahrensimplizierten Fehler, der bei jedem Integrationsschritt noch zulässig sein soll. Dadurch wirkt sich der Wert von E# direkt auf die Schrittweitensteuerung aus.
Diese Prozedur löst das Anfangswertproblem einer Differentialgleichung erster Ordnung nach dem Runge-Kutta-Verfahren. Die Prozedur übernimmt dabei die komplette Schrittweitensteuerung.

Im Ordner DEMO finden Sie das Programm "DiffEquation.BAS", welches das Runge-Kutta-Verfahren anhand eines Beispiels mit dem Adam-Bashforth-Verfahren vergleicht.
 
 
FN Runge_Step#(&FN Fun#(0,0,0),I,X#,Y#,S#)
Fun#(I,X#,Y#) Diese Funktion definiert die Differentialgleichung über die Beziehung Y#'=FN Fun#(I,X#,Y#). Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
I Ein Index, der an FN Fun# übergeben wird.
X# X-Wert, für den der Schritt ausgeführt werden soll.
Y# Y-Wert, für den der Schritt ausgeführt werden soll.
S# In diesem Parameter müssen Sie die Schrittweite übergeben.
Diese Funktion führt einen einzelnen Schritt nach dem Runge-Kutta-Verfahren durch. Sie können diese Funktion für ein eigenes Schrittweiten-Steuerungsprogramm verwenden oder einfacher gleich die Prozedur Runge_Kutta verwenden.
 
 
Adam_Bash &FN Fun#(0,0,0),I,Xs#,Xe#,Y#,&X#(,),N
Fun#(I,X#,Y#) Diese Funktion definiert die Differentialgleichung über die Beziehung Y#'=FN Fun#(I,X#,Y#). Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
I Ein Index, der an FN Fun# übergeben wird.
Xs# Anfangswert für X#.
Xe# Endwert für X#.
Y# Vor dem Aufruf müssen Sie hier den Anfangswert für Y# übergeben.
X#(0:N,0:1) In diesem Feld erhalten Sie die Lösungswerte Ihrer Differentialgleichung an den äquidistanten Stützstellen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
N Höchster Index des Feldes X#(,) (Anzahl der Stützstellen minus eins).
Diese Prozedur löst das Anfangswertproblem einer Differentialgleichung erster Ordnung nach dem Adams-Bashforth-Vierschritt-Verfahren. Die Prozedur arbeitet etwas schneller als Runge_Kutta, funktioniert aber nur mit einer festen Schrittweite. Deshalb ist sie besonders gut zu verwenden,wenn Sie die Lösung einer Differentialgleichung zwischen Anfangs- und Endwert an äquidistanten Stützstellen erhalten wollen.

Im Ordner DEMO finden Sie das Programm "DiffEquation.BAS", welches das Runge-Kutta-Verfahren anhand eines Beispiels mit dem Adam-Bashforth-Verfahren vergleicht.
 
 
FN Adam_Step#(&FN Fun#(0,0,0),I,X#,Y#,Y0#,Y1#,Y2#,Y3#,S#)
Fun#(I,X#,Y#) Diese Funktion definiert die Differentialgleichung über die Beziehung Y#'=FN Fun#(I,X#,Y#). Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
I Ein Index, der an FN Fun# übergeben wird.
X# X-Wert, für den der Schritt ausgeführt werden soll.
Y# Y-Wert, für den der Schritt ausgeführt werden soll.
Y0#,Y1#,
Y2#,Y3#
Hier müssen Sie die Lösungswerte Ihrer Differentialgleichung an den letzten vier Stützstellen übergeben.
S# In diesem Parameter müssen Sie die Schrittweite übergeben.
Diese Funktion führt einen einzelnen Schritt nach dem Adams-Bashforth-Vierschritt-Verfahren durch. Sie können diese Funktion für ein eigenes Lösungsprogramm verwenden oder einfacher gleich die Prozedur Adam_Bash verwenden.
 
 
Runge_Step_Sys &FN Fun#(0,0,0,0),I,X#,&Y#(),N,S#
Fun#(I,X#,Y,Yd) Diese Funktion definiert das Differentialgleichungssystem. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
X# ist die unabhängige Variable.
Y ist die Adresse eines Fließkommafeldes, das die abhängigen Variablen representiert.
Yd ist die Adresse eines Fließkommafeldes, in dem Ihre Funktion die ersten Ableitungen der abhängigen Variablen zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
X# X-Wert, für den der Schritt ausgeführt werden soll.
Y#(0:N) In diesem Feld müssen die Anfangswerte der abhängigen Variablen übergeben werden.
N Höchster Index des Feldes Y#() (Anzahl der Differentialgleichungen minus eins).
S# In diesem Parameter müssen Sie die Schrittweite übergeben.
Diese Prozedur führt einen einzelnen Schritt nach dem Runge-Kutta-Verfahren für Differentialgleichungssysteme durch. Sie können diese Prozedur für ein eigenes Lösungsprogramm verwenden und dabei mit einer konstanten Schrittweite arbeiten oder eine an das Gleichungssystem angepaßte Schrittweitensteuerung verwenden.

Im Ordner DEMO finden Sie das Programm "RungeKuttaSys.BAS", welches das Runge-Kutta-Verfahren für Differentialgleichungssysteme anhand eines Beispiels erläutert und dabei auch auf den Einfluß unterschiedlicher Schrittweiten eingeht.
 
 
 
 
De &FN Fun#(0,0,0,0),I,R Xs#,Xe#,&Y#(),N,R Re#,R Ae#,R F
Fun#(I,X#,Y,Yd) Diese Funktion definiert das Differentialgleichungssystem. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
X# ist die unabhängige Variable.
Y ist die Adresse eines Fließkommafeldes, das die abhängigen Variablen representiert.
Yd ist die Adresse eines Fließkommafeldes, in dem Ihre Funktion die ersten Ableitungen der abhängigen Variablen zurückliefern muß.
I Ein Index, der an FN Fun# übergeben wird.
Xs#,Xe# Diese beiden Parameter definieren das Intervall, über das integriert werden soll. Dabei kann Xs# auch größer als Xe# sein, da die Prozedur auch rückwärts integrieren kann.
Nach der Rückkehr enthält
Xs# den zuletzt verwendeten Wert. Diese Information kann für Sie wichtig sein, falls die Integration abgebrochen werden mußte.
Y#(0:N) In diesem Feld müssen die Anfangswerte der abhängigen Variablen übergeben werden. Nach der Rückkehr enthält dieses Feld den Lösungsvektor an der Stelle Xe#.
N Höchster Index des Feldes Y#() (Anzahl der Differentialgleichungen minus eins).
Re#,Ae# Mit diesen Parametern legen Sie den erlaubten relativen Fehler (Re#) und den absuluten Fehler (Ae#) fest. Die Prozedur verwendet dabei folgende Ungleichung:

ABS(Absolute_Error#)<=Re#*ABS(Y#(0:N))+Ae#

Durch Nullsetzen von Re# bzw. Ae# ist es möglich, sowohl ein reines absolutes als auch ein reines relatives Fehlerkriterium zu erhalten. Das wird jedoch nicht oft sinnvoll sein, da man bei den meisten Differentialgleichungen davon ausgehen muß, daß kaum Informationen über das Verhalten der Lösung bekannt sind. Aus diesem Grunde ist es normalerweise sinnvoll, sowohl einen relativen als auch einen absoluten Fehler zu definieren. Bei einer betragsmäßig großen Schwankung der Lösung wird dann hauptsächlich Re# für die Steuerung des Integrators verwendet, bei einer betragsmäßig kleinen Lösung Ae#. Damit sind Sie also immer auf der sicheren Seite.

F Über diese Variable kommuniziert die Prozedur mit Ihrem Program. Im einzelnen haben die Werte von F (-1 oder 1 bis 6) folgende Bedeutung:

F=-1,1:
Diesen Wert muß der Benutzer setzen, wenn er die Integration neu starten will, wenn also eine Integration neu begonnen werden soll oder die Differentialgleichung z.B. eine Unstetigkeitsstelle besitzt, über die man nicht hinwegintegrieren kann. Im Fall der Unstetigkeit ist es auch sinnvoll, Mit
F=-1 zu beginnen. Dadurch wird der Integrator angewiesen, nicht über Xe# hinauszuintegrieren, die Unstetigkeit also nicht zu überschreiten. Aus Gründen der optimalen Schrittweitenwahl versucht der Algorithmus bei positiven Werten von F nämlich nicht, Xe# exakt zu erreichen. Vielmehr wird der Lösungsvektor des Systems an der Stelle Xe# normalerweise durch Interpolation ermittelt. Das hat aber keinen Genauigkeitsverlust zur Folge.

Während des Integrationsprozesses kann F verschiedene Werte annehmen, die Ihr Programm über den Verlauf der Integration informieren. Diese Werte dürfen Sie nicht ändern.

F F=2:
Dieser Wert signalisiert einen erfolgreichen Abschluß der Integration.
Xs# ist gleich Xe#, so daß Sie Xe# vor dem nächsten Aufruf nur noch auf einen neuen Wert setzen müssen, um die Integration fortzusetzen. Dafür sind alle internen Felder des Integrators schon initialisiert. In Y#(0:N) steht der Lösungsvektor des Differentialgleichungssystems an der Stelle Xe#.

F=3:
Re# und Ae# sind für die vorhandene Rechengenauigkeit (16 Stellen) zu klein. Xs# wird auf den Wert gesetzt, der Xe# während der Integration am nächsten lag und Y#(0:N) enthält den Lösungsvektor an dieser Stelle. Re# und Ae# wurden von De auf annehmbare Werte gesetzt, so daß der Integrator nur noch neu aufgerufen werden muß, um die Integration mit diesen höheren Fehlerschranken fortzusetzen.

F F=4:
Die Integration wurde abgebrochen, weil mehr als 500 interne Schritte (d.h. 1000 Auswertungen des Differentialgleichungssystems) erforderlich sind, um Xe# zu erreichen.
Xs# wird auf den Wert gesetzt, der Xe# während der Integration am nächsten lag. Ihr Programm kann jetzt entscheiden, ob es die Integration trotz der längeren Rechenzeit fortsetzen oder abbrechen will. Zur Fortsetzung der Integration braucht De nur wieder neu aufgerufen zu werden.

F=5:
Dies hat die gleiche Bedeutung wie F=4, es wurde jedoch zusätzlich festgestellt, daß die Differentialgleichungen steif zu sein scheinen. Das kann dazu führen, daß weitere Rechnungen sehr lange dauern, auch wenn die Ergebnisse immer noch gut sein können. Zur Fortsetzung der Integration braucht
De nur wieder neu aufgerufen zu werden.

F F=6:
Bei jedem Aufruf von
De wird überprüft, ob die Eingabeparameter zulässig sind. Demzufolge wird F=6 gesetzt, wenn Xs#=Xe#, Re# oder Ae# kleiner null sind oder F betragsmäßig nicht zwischen 1 und 5 liegt.

Diese Prozedur ist ein sehr häufig verwendeter Integrator zur Lösung von Differentialgleichungssystemen erster Ordnung. Es handelt sich dabei ursprünglich um ein FORTRAN-Programm, das hier in einer Omikron Basic Variante zur Verfügung gestellt wird.
Es wird nach dem Adams-Verfahren variabler Ordnung und variabler Schrittweite integriert. Die Prozedur steuert die Ordnung und die Schrittweite selbsständig und versucht dabei, den verfahrensimpliziten Fehler unterhalb eines vom Anwender definierten relativen und absoluten Fehlers zu halten.
Die Ordnung kann zwischen 1 und 13 variieren, d. h. der Algorithmus wird Differentialgleichungen, deren Lösungen Polynome maximal 13. Grades sind, noch exakt integrieren (von Rundungsfehlern einmal abgesehen).
Außerdem kann eine Steifheit der Differentialgleichungen in den meisten Fällen erkannt werden. In diesen Fällen wird die Integration abgebrochen.
Diese Prozedur ist langsamer als die weiter oben besprochenen, liefert aber im allgemeinen bessere Ergebnisse.
Wenn Sie sich sehr oft mit der numerischen Lösung von Differentialgleichungen beschäftigen, sollten Sie mal einen Blick in das Buch von Shampine & Gordon werfen. Dort wird die Herleitung und Verwendung der Prozedur
De in der FORTRAN-Variante ausführlich diskutiert.

Hinweis: Bei der Integration kann es zu Problemen kommen, wenn eine der Differentialgleichungen im Laufe der Integration genau gleich null wird, da dann der relative Fehler nicht definiert ist. Das wird aber aufgrund der internen Rechenfehler so gut wie nie vorkommen, außer die Differentialgleichungen sind schon am Anfangswert gleich null.
In diesem Fall können Sie sich damit behelfen, die entsprechenden Werte statt auf null auf sehr kleine Zahlen zu setzen (z.B 1D-10). Das erzeugt praktisch keine weiteren Rechenfehler und vermeidet diese Situation.
In kritischen Fällen kann es auch sinnvoll sein, das Steuerwort
COMPILER "DEBUG ON" zu setzen, um auf unerlaubte Fließkommaoperationen (z.B. Division durch 0) aufmerksam gemacht zu werden.

Für die folgenden Beispiele wurden leicht durchschaubare Differentialgleichungen verwendet, die auch analytische Lösungen haben, mit denen die numerischen Ergebnisse verglichen werden können:

Beispiel 1:

Num_Init
'Lösung von Y'=-Y.
DIM Y#(0),Yd#(0)
'Anfangswerte.
Xs#=0:Xe#=10
Y#(0)=1
Re#=1D-5:Ae#=0
F=1
De &FN Diff_Equation#(0,0,0,0),0,Xs#,Xe#,&Y#(),0,Re#,Ae#,F
PRINT
PRINT "Die numerische Lösung lautet:.. ";
PRINT "X=";Xs#,"Y=";Y#(0)
PRINT
PRINT "Die theoretische Lösung lautet: ";
PRINT "X=";Xs#,"Y=";EXP(-10)
PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Diff_Equation#(I,X#,Y,Yd)
 *Yd#(0)=-*Y#(0)
RETURN 0


Beispiel 2:
Die Kreisbewegung eines masselosen Körpers um ein Massenzentrum wird durch die Differentialgleichungen
Y0''=-Y0*R^-3 und Y1''=-Y1*R^-3 beschrieben. Hierbei handelt es sich um ein System zweiten Grades, das ja bekanntlich immer in ein System ersten Grades umgeformt werden kann. In diesem Fall ergeben sich vier Differentialgleichungen ersten Grades:

Num_Init
'Lösung eines himmelsmechanischen Problems.
DIM Y#(3),Yd#(3)
'Anfangswerte.
Xs#=0:Xe#=2*PI:'Entspricht einem Umlauf.
Y#(0)=1:Y#(1)=0:'Anfangsposition.
Y#(2)=0:Y#(3)=1:'Anfangsgeschwindigkeit.
Re#=1D-6:Ae#=1D-6
F=1
REPEAT
 De &FN Orbit#(0,0,0,0),0,Xs#,Xe#,&Y#(),3,Re#,Ae#,F
UNTIL F=2 or F=6
If F=6 THEN PRINT "Integration wurde abgebrochen":GOTO Ende
PRINT
PRINT "Numerische Lösung:","Theoretische Lösung:"
PRINT
PRINT Y#(0),COS(Xe#)
PRINT Y#(1),SIN(Xe#)
PRINT Y#(2),-SIN(Xe#)
PRINT Y#(3),COS(Xe#)
-Ende:PRINT:INPUT "Ende mit [Return]";Dummy
Num_Exit
END

DEF FN Orbit#(I,X#,Y,Yd)
 LOCAL R#=(*Y#(0)^2+*Y#(1)^2)^-1.5:'Abstand der Körper.
 *Yd#(0)=*Y#(2):'Geschwindigkeit ist die
 *Yd#(1)=*Y#(3):'Ableitung des Ortes.
 *Yd#(2)=-*Y#(0)/R#:'Beschleunigung ist die
 *Yd#(3)=-*Y#(1)/R#:'Ableitung der Geschwindigkeit.
RETURN 0


Im Ordner DEMO finden Sie das Programm "Riccati.BAS", das an einem weiteren Beispiel die vielseitige Verwendbarkeit dieser Prozedur zeigt.
 
 




2.10 Graphische Ausgabe
 
Die grapgische Darstellung von Funktionen oder Meßwerten ist in den meisten Fällen viel aussagekräftiger als eine tabellarische Auflistung. Darum gibt es in der Numeric Library einige Prozeduren, die Funktionen oder einzelne Werte in Koordinatenkreuze einzeichnen. Wir haben uns dabei bemüht, diese Prozeduren so flexibel wie möglich zu halten, damit sie für möglichst viele verschiedene Fälle eingesetzt werden können. Dennoch stellen vorgefertigte Routinen immer einen Kompromiß dar zwischen einer einfachen Handhabung und einer möglichst vielfältigen Einsatzmöglichkeit.
 
 
Plot_Coordinates X,Y,W,H,X1#,X2#,Y1#,Y2#,F
X X-Wert der linken oberen Ecke des Koordinatensystems in Pixel.
Y Y-Wert der linken oberen Ecke des Koordinatensystems in Pixel.
W Breite des Koordinatensystems in Pixel.
H Höhe des Koordinatensystems in Pixel.
X1# Minimaler X - Wert.
X2# Maximaler X - Wert.
Y1# Minimaler Y - Wert.
Y2# Maximaler Y - Wert.
F Mit diesem Parameter können Sie festlegen, welche Elemente gezeichnet werden sollen. Jedes Bit dieser Integer-Zahl steht dabei für ein bestimmtes Element, das nur dann gezeichnet wird, wenn das entsprechende Bit gesetzt ist. Die einzelnen Bits haben die folgende Bedeutung:

Bit 0 : Bereich des Koordinatensystems füllen.
Bit 1 : Rahmen um das Koordinatensystem zeichnen.
Bit 2 : Koordinatenkreuz einzeichnen.
Bit 3 : Koordinatenachsen beschriften.
Bit 4 : Hilfslinien einzeichnen.

Diese Prozedur zeichnet ein Koordinatenkreuz mit verschiedenen Attributen. Dabei werden für die Zeichenelemente die aktuellen Einstellungen für Linienbreite, Linienfarbe, Füllfarbe usw. verwendet. Sie müssen also selbst dafür sorgen, daß diese Attribute geeignet eingestellt sind. Andererseits haben Sie dadurch auch die Möglichkeit, geziehlt Einfluß zu nehmen auf das Aussehen Ihres Koordinatenkreuzes.
Wenn Sie z.B. erreichen möchten, daß die Koordinatenachsen mit dickeren Linien gezeichnet werden als die Hilfslienen, stellen Sie die Linienbreite zunächst auf 2, zeichen das Koordinatensystem ohne die Hilfslinien, stellen dann auf eine Linienbreite von 1 um und rufen die Prozedur nochmal auf, wobei beim zweiten mal die Flags so gesetzt werden müssen, daß nur die Hilfslinien gezeichnet werden.

Achtung: Diese Prozedur setzt das Clipping-Rechteck zurück. Falls Sie also selbst ein Clipping-Rechteck verwenden, müssen Sie es nach dem Aufruf neu einstellen.

Im Ordner DEMO finden Sie diverse Programme, die diese Prozedur verwenden.
 
 

Plot_Function &FN Fun#(0,0),I,X,Y,W,H,X1#,X2#,Y1#,Y2#,F[,S]
Fun#(I,X#) Dies ist die Funktion, deren Graph gezeichnet werden soll. Diese Funktion müssen Sie natürlich selbst definieren und dann die Adresse übergeben.
In
I wird Ihrer Funktion ein Parameter übergeben, den Sie für eigene Zwecke verwenden können, z.B. um verschiedene Varianten einer Funktion zu identifizieren.
In
X# wird das Argument übergeben, dessen Funktionswert Ihre Funktion zurückliefern muß.

I Ein Index, der an FN Fun# übergeben wird.
X X-Wert der linken oberen Ecke des Koordinatensystems in Pixel.
Y Y-Wert der linken oberen Ecke des Koordinatensystems in Pixel.
W Breite des Koordinatensystems in Pixel.
H Höhe des Koordinatensystems in Pixel.

X1# Minimaler X - Wert.
X2# Maximaler X - Wert.
Y1# Minimaler Y - Wert.
Y2# Maximaler Y - Wert.

F Mit diesem Parameter können Sie festlegen, welche Elemente gezeichnet werden sollen. Jedes Bit dieser Integer-Zahl steht dabei für ein bestimmtes Element, das nur dann gezeichnet wird, wenn das entsprechende Bit gesetzt ist. Die einzelnen Bits haben die folgende Bedeutung:

Bit 0 : Bereich des Koordinatensystems füllen.
Bit 1 : Rahmen um das Koordinatensystem zeichnen.
Bit 2 : Koordinatenkreuz einzeichnen.
Bit 3 : Koordinatenachsen beschriften.
Bit 4 : Hilfslinien einzeichnen.
S Mit diesem Parameter können Sie festlegen, alle wieviel Pixel in X-Richtung ein Funktionswert berechnet wird. Wenn Sie diesen Parameter weglassen oder einen Wert <=1 übergeben wird S=1 gesetzt, es wird also ein Funktionswert pro Pixel berechnet.
In Fällen, bei denen die Berechnung der Funktionswerte sehr aufwendig ist, kann es sinnvoll sein, für
S einen größeren Wert zu wählen. Siehe z.B. das Programm Riccati.BAS im DEMO-Ordner.

Diese Prozedur zeichnet den Graphen einer Funktion. Dabei kann auch gleich noch ein Koordinatenkreuz mitgezeichnet werden. Für die graphischen Elemente werden die aktuellen Einstellungen für Linienbreite, Linienfarbe, Füllfarbe usw. verwendet. Sie müssen also selbst dafür sorgen, daß diese Attribute geeignet eingestellt sind.
Wenn man graphisch anspruchsvollere Darstellungen erhalten möchte, sollte man zunächst mit
Plot_Coordinates nur das Koordinatensystem zeichnen, um dann mit geänderter Linienfarbe etc. den eigentlichen Graphen bei abgeschalteten Attributen (F=0) einzuzeichnen.

Achtung: Diese Prozedur setzt das Clipping-Rechteck zurück. Falls Sie also selbst ein Clipping-Rechteck verwenden, müssen Sie es nach dem Aufruf neu einstellen.

Im Ordner DEMO finden Sie diverse Programme, die diese Prozedur verwenden.
 
 

Plot_Array &X#(0,0),N,X,Y,W,H,X1#,X2#,Y1#,Y2#,F
X#(0:N,0:1) In diesem Feld übergeben Sie die Punkte, die gezeichnet werden sollen. Dabei stehen in X#(0:N,0) die X-Werte und in X#(0:N,1) die Y-Werte.
N Höchster Index des Feldes X#(,) (Anzahl der Punkte minus eins).
X X-Wert der linken oberen Ecke des Koordinatensystems in Pixel.
Y Y-Wert der linken oberen Ecke des Koordinatensystems in Pixel.
W Breite des Koordinatensystems in Pixel.
H Höhe des Koordinatensystems in Pixel.
X1# Minimaler X# - Wert.
X2# Maximaler X# - Wert.
Y1# Minimaler Y# - Wert.
Y2# Maximaler Y# - Wert.

F Mit diesem Parameter können Sie festlegen, welche Elemente gezeichnet werden sollen. Jedes Bit dieser Integer-Zahl steht dabei für ein bestimmtes Element, das nur dann gezeichnet wird, wenn das entsprechende Bit gesetzt ist. Die einzelnen Bits haben die folgende Bedeutung:

Bit 0 : Bereich des Koordinatensystems füllen.
Bit 1 : Rahmen um das Koordinatensystem zeichnen.
Bit 2 : Koordinatenkreuz einzeichnen.
Bit 3 : Koordinatenachsen beschriften.
Bit 4 : Hilfslinien einzeichnen.

Die Bits 8 - 10 legen fest, ob die einzelnen Punkte mit Linien verbunden oder als geometrische Einzelobjekte gezeichnet werden sollen.

Bit 10 9 8
    0 0 0 : Punkte durch Linien verbinden.
    0 0 1 : x-förmige Kreuze.
    0 1 0 : normale Kreuze.
    0 1 1 : Auf der Basislinie stehende Dreiecke.
    1 0 0 : Auf der Spitze stehende Dreiecke.
    1 0 1 : Auf der Spitze stehende Quadrate.
    1 1 0 : Auf der Seitenlinie stehende Quadrate.
    1 1 1 : Kreise.

Die Möglichkeit, verschiedene Punktformen definieren zu können, ist besonders dann von Vorteil, wenn man mehrere Kurven in ein Koordinatenkreuz einzeichnen möchte, zur Unterscheidung aber bei einer Schwarz-Weiß-Darstellung keine Farben zur Verfügung hat. Dann kann man die einzelnen Kurven durch verschiedene Punktformen unterscheiden.

Tip: Durch vergrößern der Linienbreite können Sie auch ausgefüllte Objekte erzeugen.

Diese Prozedur zeichnet diskrete Punkte (z.B. Meßwerte), die entweder durch Linien verbunden oder als geometrische Objekte dargestellt werden können. Dabei kann auch gleich noch ein Koordinatenkreuz mitgezeichnet werden. Für die graphischen Elemente werden die aktuellen Einstellungen für Linienbreite, Linienfarbe, Füllfarbe usw. verwendet. Sie müssen also selbst dafür sorgen, daß diese Attribute geeignet eingestellt sind.
Wenn man graphisch anspruchsvollere Darstellungen erhalten möchte, sollte man zunächst mit
Plot_Coordinates nur das Koordinatensystem zeichnen, um dann mit geänderter Linienfarbe etc. die eigentlichen Punkte bei abgeschalteten Attributen (F=0) einzuzeichnen.

Achtung: Diese Prozedur setzt das Clipping-Rechteck zurück. Falls Sie also selbst ein Clipping-Rechteck verwenden, müssen Sie es nach dem Aufruf neu einstellen.

Im Ordner DEMO finden Sie diverse Programme, die diese Prozedur verwenden.
 
 
 

Allgemeines -blättern- Inhaltsverzeichnis

Omikron Basic im Internet: http://www.berkhan.de

Copyright 1998-1999 by Berkhan-Software