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
PR