Turbo Toolbox

home *** CD-ROM | disk | FTP | other *** search

/ Turbo Toolbox / Turbo_Toolbox.iso / spezial / 15 / dok.asc next >

Wrap

Text File | 1989-06-15 | 79.6 KB | 2,022 lines

Numerische Mathematik anschaulich dargestellt von Heinz Hagemeyer Karsten Gieselmann und Peter Kurzweil Der Experimentator hat Daten gesammelt. Wie sind diese auszuwerten. Welche Gesetzmäßigkeiten liegen den Daten zugrunde. Was ist zu tun. Handelt es sich um eine lineare-, quadratische oder sonstige Funktion ? Liegen die Daten exakt vor oder sind sie mit Streuungen oder gar mit Ausreißern behaftet? Welche Methoden können angewendet werden. Gibt es Möglichkeiten die Daten aufzubereiten und sie dann graphisch darzustellen? Strategien zur Auswertung von Daten zu entwickeln ist u.a. die Aufgabe der numerischen Mathematik. Aber auch Integrale, die nicht geschlossen, d.h. formelmäßig gelöst werden können, können zumindest numerisch berechnet werden. Wie aber sieht es mit Differenzialgleichungen aus. Wer nicht weiß, was das ist, für den ist dieser Beitrag genau das Richtige! Zum Schluß lernen Sie noch eine andere Möglichkeit zum Lösen linearer Gleichungssysteme kennen. Damit ist das Gebiet der numerischen Mathematik noch lange nicht erschöft. Der Kenner hoft auf eine Fortsetzung. Inhalt : Wertetabellen graphisch darstellen······················· Seite 2 Beispiele : - Lösen von Differenzialgleichungen :·················Seite 5 - Freier Fall mit Luftwiderstand - - Aufladung eines Kondensators - - Gedämpfte Schwingung - - Differenzieren und Integrieren·-··················· Seite 8 Regression und Interpolation····························· Seite 10 Lineare Regression·····································Seite 10 - Potenzfunktionen - - Exponentialfunktionen - - Logarithmusfunktionen - Quadratische Regression······························· Seite 15 Interpolation mit kubischen Splines··················· Seite 16 Numerische Integration··································· Seite 19 Iterative Lösung eines linearen Gleichungssystems·········Seite 26 --------------------------------- Wertetabellen graphisch darstelen --------------------------------- Es kommt oft vor, daß Daten in Form von Werteta- bellen auszuwerten sind. Ein klassisches Beispiel hierfür ist, daß von einer Funktion eine Wer- tetabelle berechnet wurde mit dessen Hilfe an- schließend der Funktionsgraph zu zeichen ist. Klar, daß diese Wertetabelle einfach per Programm erstellt werden kann. Aber nicht immer sind die Daten so einfach zu gewinnen. Es gibt zahlreiche Fälle, wo es so einfach nicht geht. Man denke nur an empirisch ermittelte Daten, an numerische Lösungen von Differenzialgleichungen, an Ableitungen, die nicht ganz einfach zu gewinnen sind etc.. In all diesen Fällen ist es mühsam die Meßpunkte in ein Koordinatensystems einzuzeichnen um sie dann mit Hilfe von Kurvenlinealen zu verbinden. Wer hat bei dieser Arbeit noch nie geflucht ? -------- Die Idee -------- ist, diese Arbeit sich durch ein Programm abnehmen zu lassen. Diese Wertetabelle muß selbstver- ständlich ein bestimmtes Format besitzen. Um das Programm so universell wie möglich einsetzen zu können, besteht dieses Format aus ASCII - Zeichen. Als erster Wert einer Zeile muß der x - Wert (die unabhängige Variable) in beliebiger Darstellung (z.B. Gleitkommadarstellung wie 3.4E-02 oder in "normaler" Darstellung wie 123.76) eingetragen werden. Es folgt mindestens ein Leerzeichen und dann der y - Wert. Auf diese Weise kann die Wertetabelle in einer beliebigen Sprache wie BASIC, FORTRAN, C, PASCAL oder - mit Hilfe eines Editors - von Hand erzeugt werden. Wird dem Wert ein großes "P" vorangestellt, so wird dieser nicht als Punkt des Kurvenzuges sondern als singulärer Punkt interpretiert. Bei- spiel : P 0.0000000000E+00 0.0000000000E+00 P 1.0000000000E-01 8.3333333333E-01 P 2.0000000000E-01 1.5972222222E+00 0.0000000000E+00 0.0000000000E+00 1.0000000000E-01 8.3333333333E-01 2.0000000000E-01 1.5972222222E+00 3.0000000000E-01 2.2974537037E+00 4.0000000000E-01 2.9393325617E+00 5.0000000000E-01 3.5277215149E+00 6.0000000000E-01 4.0670780553E+00 7.0000000000E-01 4.5614882174E+00 8.0000000000E-01 5.0146975326E+00 9.0000000000E-01 5.4301394049E+00 1.0000000000E+00 5.8109611211E+00 Das Programm liest nun Zeile für Zeile diese Wertetabelle ein. Die mit einem "P" versehenen Wertepaare werden als Punkt (kleiner Kreis) in ein Koordinatensystem eingezeichnet, die anderen Werte werden durch Geraden verbunden. (Hier liegt auch die Schwäche des Algorithmus, die Punkte müssen dicht beieinander liegen sonst entsteht kein "glatter" Kurvenzug sondern eine "Fieberkurve"). Wird in die Tabelle eine Leerzeile eingefügt, so werden die nächsten Daten als neuer Kurvenzug interpretiert. Auf diese Art und Weise ist es möglich, mehrere Funktionsgraphen in einem Koordinatensytem darzustellen. Dies ist ebenfalls möglich durch die Eingabe mehrerer Dateien. Die Anzahl wird durch die Konstante MaxDateiAnzahl begrenzt werden müssen. Das Programm soll, kennt man die auszuwertende Datei(en) vorher, mit den(m) Dateinamen aufgerufen werden. Also so z.B.: GRAPHEN S.DAT V.DAT Dann werden durch das Programm die angegebenen Dateien ausgewertet. Ohne Parameter aufgerufen fordert das Programm die auszuwertenden Dateien "per Hand" an. Das Programm eröffnet sicherlich eine große Zahl von Anwendungen. Einige von ihnen werden desweiten beschrieben. Viel Spaß beim Ausdenken neuer Möglichkeiten. --------------- Algorithmisches --------------- Das Programm, welches diese Aufgabe übernimmt, finden Sie auf der Sonderdiskette unter dem Namen GRAPHEN.PAS. Es ist in der Sprache TURBO - Pascal 4.0 unter Verwendung der Unit GRAPH und der UNIT PGraph aus der Toolbox 12/89 entwickelt worden. Wenn das Programm ohne Paramer aufgerufen wird, wird zunächst abgefragt, welche Dateien ausgewertet werden sollen. Falls sich die Dateien nicht im aktuellen Inhaltsverzeichnis befinden oder nicht vorhanden sind, gibt die Procedure DateiEingabe eine entsprechende Fehlermeldung aus. Die eingegebenen Namen werden in der Variablen Datei_Name gespeichert, die Anzahl der eingebenen Dateien in der Variablen DateiAnzahl. Die Procedure wird Verlassen durch Eingabe von RETURN (Leerzeile) oder wenn MaxDateiAnzahl überschritten wird. Danach müssen die Parameter zum Zeichnen des Koordinatensystems eingegeben werden. Die Grafik wird initialisiert und das Koordinatensystem mit den vorher eingegebenen Werten gezeichnet. Die Da- teien werden nun nacheinander verarbeitet, wobei auf dem Bildschirm der Dateiname ausgegeben wird. Dabei wird zunächst der erste Punkt ( = erste Zeile der ASCII Datei) eingelesen. In der PROCEDURE Verwandle werden x und y - Wert voneinander getrennt und in reelle Zahlen (Datentyp Real) umgewandelt. Falls dabei ein Fehler auftritt, wird dies als Fehlermeldung ausgegeben. Die Umrechnung in die Bildschirmkoordinaten übernimmt die Procedure Umrechnung, die wiederum die Procedure Scale aus PGraph.pas aufruft. Dies ist so notwendig, weil PGraph ohne Clipping arbei- tet. Liegt ein Punkt außerhalb des Koordinatensystems, würde PGraph ihn auch dort darstellen. Die Anweisung SetViewPort in der Procedure KoordViewPort verhindert dieses durch ClipOn. Scale liefert dann leider falsche Werte zurück, die in der Procedure Umrechnung korrigiert werden. Wer will, mag Scale gleich passend ändern. Ist dem Punkt ein "P" vorangestellt, wird dieser durch die Procedure Mark durch einen kleinen normgerechten Kreis im Koordinatensystem darge- stellt. Ansonsten werden die x und y Werte den Variablen Xalt und Yalt zugewiesen. Sobald das "P" nicht mehr auftritt, wird der erste Punkt in den Variablen Xalt und Yalt und der zweite in Xneu und Yneu gespeichert. Mit Hilfe dieser 4 Variablen kann nun die erste Linie ge- zogen werden. Danach werden den "alten" Werten die "neuen" Werte zugewiesen und so wird Linie für Linie gezogen, bis EOF (End Of File) erreicht wird. ------------------------------- Das entscheidene Struktogramm : ------------------------------- ┌────────────────────────────────────────────────────────────┐ │ Solange Datei nicht zu Ende │ │ ┌────────────────────────────────────────────────────┤ │ │ Lese erste Zeile │ │ ├────────────────────────────────────────────────────┤ │ │ Verwandle die Zeile in die Zahlen Xalt, Yalt │ │ ├────────────────────────────────────────────────────┤ │ │ Solange Datei nicht zu Ende │ │ │ ┌───────────────────────────────────────────────┤ │ │ │ Lese nächste Zeile │ │ │ ├───────────────────────────────────────────────┤ │ │ │ Verwandle die Zeile in die Zahlen Xneu,Yneu │ │ │ ├───────────────────────────────────────────────┤ │ │ │ Ziehe Linie P(Xalt/Yalt) - P(Xneu/Yneu) │ │ │ ├───────────────────────────────────────────────┤ │ │ │ Xalt := Xneu ; Yalt := Yneu │ └───────┴────┴───────────────────────────────────────────────┘ Liest das Programm nun eine Leerzeile ein werden die anschliesend eingelesenen X und Y - Werte den Variblen Xalt (statt Xneu) und Yalt zugewiesen. Das bedeutet, daß ein neuer Kurvenzug gezeichnet wird ! Großer Wert wurde auf die Fehlerbehandlung gelegt, da die Dateien auch von Hand oder mit anderen Programmiersprachen erstellt werden können. Die erledigt die Procedure Fehler welche detalierte Fehlermeldungen an den Benutzer ausgibt. Man sieht, auch mit relativ einfachen Algorithmen lassen sich durchaus brauchbare Programme entwickeln. Es muß nicht immer gleich Backtracking oder ähnliches sein. -------------------------------------------------------------- - Lösen von Differenzialgleichungen mit numerischen Methoden - -------------------------------------------------------------- --------------- 1. Freier Fall --------------- Der freie Fall eines Körpers wird als klassisches Beispiel in jedem Physikunterricht behandelt. Dabei bleibt der Luftwiderstand des Körpers jedoch unberücksichtigt. Wie verhält sich nun ein Körper mit Berücksichtigung des Luftwiderstandes ? Diese Frage kann mit elementaren Methoden der Physik und Mathematik nicht beantwortet werden, oder doch? (vergl. Abb. 1) Betrachten wir den Körper zum Zeitpunkt t = 0. Dann wirkt auf ihn - da er ja noch nicht fällt - nur die Erdbeschleunigung g = 9,81 m / s2. Für einen "ganz kleinen" Zeitabschnitt t ändert sich daran so gut wie nichts, so daß wir nach diesem Zeitabschnitt folgende Geschwindigkeit v erhalten : v(t) = g ∙ t . Durch diese Geschwindigkeit (obschon sie noch recht klein ist) macht sich nun der Luftwiderstand bemerkbar. Durch Beobachtung hat man festgestellt, daß dieser proportional dem Quadrat der Geschwindigkeit ist. Es gilt also : FL = - kw ∙ v (t) 2 Das "-" Zeichen drückt aus, daß die durch den Luftwiderstand FL wirkende Kraft entgegengesetzt der Kraft ist, die durch die Erdbeschleunigung hervorgerufen wird. Teilt man diese Gleichung auf beiden Seiten durch die Masse m des Körpers, so erhält man nach dem dynamischen Grundgesetz : kw aL = - cw ∙ v (t)2 [ mit cw = -- ! ] m und somit wirkt auf den Körper nunmehr die Gesamtbeschleunigung: a = g - aL Mit dieser resultierenden Gesamtbeschleunigung erhalten wir für den nächsten Zeitabschnitt : v(t + t) = v(t) + a ∙ t Daraus läßt sich nun wieder die neue Gesamtbeschleunigung berechnen usw. usf.. Um bei diesem Verfahren brauchbare Ergebnisse zu erhalten, muß man die Zeitabschnitte t nur klein genug wählen. Das bedeutet, daß man sehr viele Berechnungen durchführen muß, um z.B die Endge- schwindigkeit zu ermitteln. Dies ist nur durch den Einsatz eines programmierbaren Rechners in einem vertretbaren Zeitaufwand möglich. Mit den elementaren Methoden der Mathematik, d.h. ohne Kenntnisse der Differenzialrechnung und speziell der Differenzialgleichungen ließ sich das oben beschriebene Problem vor dem Computerzeitalter in der Schule nicht lösen obschon ein freifallender Körper im luftleeren Raum doch ziemlich unwarscheinlich ist und weit an der Realität vorbei geht. Damit erhalten wir folgendes Struktogramm : ┌┬─────────────────────────────────────────────────────────────┬┐ ││ Dateien zum Schreiben vorbereiten ││ ├┼─────────────────────────────────────────────────────────────┼┤ ││ Anfangswerte definieren und Eingabe des Luftwiderstandes ││ ├┼─────────────────────────────────────────────────────────────┼┤ ││ Anfangswerte s(0) und v(0) in die Dateien schreiben ││ ├┴─────────────────────────────────────────────────────────────┴┤ │ Zähle mit ts von 0 bis 5 (* Zeit in Sekunden *) │ │ ┌──────────────────────────────────────────────────────┤ │ │ Zähle mit tn von 1 bis Round (1 / t) │ │ │ ┌──────────────────────────────────────────────┤ │ │ │ t := ts + t ∙ tn │ │ und │ ├──────────────────────────────────────────────┤ │ tue │ und │ a := g - cw ∙ v² │ │ │ tue ├──────────────────────────────────────────────┤ │ │ │ v := v + a ∙ t │ │ │ ├──────────────────────────────────────────────┤ │ │ │ s := s + v ∙ t │ │ │ ├┬────────────────────────────────────────────┬┤ │ │ ││ Schreibe s und v in die Dateien ││ ├────────┴───────┴┴────────────────────────────────────────────┴┤ │ Dateien schließen │ └───────────────────────────────────────────────────────────────┘ Die doppelte Zählschleife bedarf einer Erklärung: Selbstverständlich ließe sich statt dessen die benötigte Zeit t inerhalb einer REPEAT - Schleife folgendermaßen zu berechnen : t := t + t; Durch ständige Wiederholung dieser Anweisung bekäme man in der Berechnung von t große Rundungsfehler, die doch sehr unschön wirken. Dies verhindern die beiden Zählschleifen. Das Programm auf der Sonderdiskette (FALL.PAS) schreibt die berechneten Daten in die beiden Dateien FALL-S.DAT und FALL-V.DAT und gibt sie zur Kontrolle noch auf dem Bildschirm aus. ------------------------------ 2. Aufladung eines Kondensator ------------------------------ Die Aufladekurve eines Kondensators ( uC(t) , iC(t) ) dürfte jedem Elektriker / Elektroniker bekannt sein. Doch wie kommt diese Kurve, dessen Gleichung in jeder besseren Formelsammlung zu finden ist, zustande ? (Abb. 2) Auch hier gehen wie wieder wie beim freien Fall vor. Zum Zeitpunkt t = 0 beträgt die Kondensatorspannung : uC(0) = 0 ( Volt ) Demnach fließt ein Strom von :·············u0 - uC ···································iC(0) = ------- ··············································R Dieser Strom lädt den Kondensator auf. Für eine "sehr kleine" Zeit t kann er als konstant angenommen werden, so daß nach dieser Zeit die Spannung am Kondensator sich mit : Q = C * U errechnet : ····································iC(t - t) uC (t) = uC(t - t) + ---------- ········································C Wählt man die Zeit t klein genug, so erhält man durchaus brauchbare Ergebnisse. Auch hier gilt, daß sich diese Berechnungen nur dann durchführen lassen, wenn man einen Computer einsetzt. Das Program KOND.PAS auf der Sonderdiskette zeigt die Realisierung dieses Ansatzes in der Sprache Pas- cal. Auch hier wird eine Wertetabelle auf dem Bildschirm ausgegeben und gleichzeitig werden die Dateien U.DAT und I.DAT erzeugt, die man dann mit GRAPHEN.PAS auswerten kann. Hier ist die maßgebliche Schleife in PASCAL : FOR ts := 0 TO te DO { Fuer jede volle Sekunde } FOR tn := 1 TO Round (1/dt) DO { ein paar mal } BEGIN t := ts + t * tn ; { Neue Werte be- } Uc := Uc + Ic/C * t ; { rechnen und } ic := (U0 - Uc) / R ; Schreibe ; { in die Datei schreiben } END; Auch hier wurde wieder eine doppelte Zählschleife zur Berechnung der Werte für t realisiert, um Rundungfehler zu verhindern. Das Schreiben in die Dateien erledigt die PROCEDURE Schreibe. Das kom- plette Programm finden Sie wieder auf der Sonderdiskette. Es läßt sich bequem in andere Sprachen übertragen. Die Darstellung der mit dem Programm erzeugten Daten mit GRAPHEN.PAS zeigt Abb. 5 ----------------------- 3. Gedämpfte Schwingung ----------------------- Bei einer freien Schwingung wird ein schwingungsfähiges System nach einmaligem Anstoß sich selbst überlassen, d.h. es wirken keine wei- teren Kräfte von außen auf das System. Dies ist wiederum sehr unrealistisch, da Luftwiderstand und Reibung die Schwingung allmählich zum Abklingen führen. Um dies zu untersuchen nehmen wir das System wie in der Abb. 3 gezeigt an. Dabei wirken auf die Masse m des Körpers folgende Kräfte : 1.) Die rücktreibende Kraft FF der Feder. Diese ist proportional zur Auslenkung s der Feder, d.h. es gilt : FF = - fk ∙ s ( fk : Federkonstante ) 2.) Eine Dämpfungskraft FD, welche proportional zur Geschwindigkeit v der Masse m ist, d.h. es gilt : FD = - D ∙ v (D : Dämpfungskonstante) Wird die Masse der Feder nicht berücksichtigt, dann ergibt sich die Beschleunigungskraft mit : F = FD + FF Somit erhält man für die auf das System wirkende Beschleunigung : ·······························s(t) ∙ FK + D ∙ v(t) ·················a(t·+·t)·=·-·-------------------- ········································m Diese Beschleunigung ist wiederum für einen "sehr kleinen" Zeitabschnitt t konstant, so daß sich daraus analog zum freien Fall die Geschwindigkeit und der zurückgelegte Weg berechnen lassen. Das Programm SCHWINGUNG.pas auf der Sonderdiskette zeigt die Realisierung dieser Überlegungen. In der Abb. 4 sieht man sehr schön, wie diese Methode dem tatsächlichem Zeitverhalten entspricht. ------------------------ 4. Infinitisimalrechnung ------------------------ Unter der Infinitisimalrechnung versteht man die Zusammenfassung von Differenzial- und Integralrechnung kurz - die klassische Analysis. Hier wurde ein Algorithmus entwickelt, der folgendes leistet : - Eine beliebige Funktion tabellieren - Von dieser Funktion die 1. und 2. Ableitung tabellieren und - Die Stammfunktion dieser Funktion tabellieren. Die dabei entstehenden Wertepaare (x/y) werden wiederum zwecks Auswertung durch GRAPHEN.PAS in eine Datei (INFI.DAT) eingeschrieben. ------------------------ Verwendete Algorithmen : ------------------------ I Differenzation : Verwendet man bei der numerischen Berechnung die Definition der Ableitung einer Funktion : ·······························f(x+x) - f(x) ···············f'(x) = lim···· -------------- ······················x -> 0········x so erhält man nur ungenaue Werte. Besser ist es, wenn die Ableitung verwendet wird, die sich nach der Taylerreihenentwicklung ergibt : ·······························f(x+x) - f(x-x) ···············f'(x) = lim···· ----------------- ·····················x -> 0······· 2 ∙ x Dies kann man sich auch anschaulich klar machen. (Vergl. Abb. 10). Entsprechend erhält man für : ·································f'(x+x) - f'(x-x) ···············f''(x) = lim····· ------------------- ·······················x -> 0········· 2 ∙ x Setzt man in dieser Formel den o.a. Term ein, so erhält man die im Programm verwendeten Gleichungen. II. Integration : Die Integration wurde nach der Doppelstreifenregel berechnet. Diese ist die Grundform der Simpsonschen Regel, welche man in jedem guten Mathematikbuch findet. (U.a. Krohn - Rattay : Mathematik für technische Berufe, Verrlag Handwerk u. Technik). Näheres entnehme man auch dem Beitrag von Carsten Gieselmann : "Numerische Integration" . Abb. 6 zeigt die Auswertung der durch INFINI gewonnenen Daten am Beispiel der Funktion f : f(x) = x² Das Programm befindet sich unter dem Namen INFI.PAS auf der Sonderdiskette. ---------- REGRESSION & Interpolation ---------- - Auswertung von Meßwerten - Wer hat noch nie vor der Aufgabe gestanden, ein paar Meßpunkte - gleich wie diese entstanden sind - in ein Koordinatensystem zu bringen. Und wer hat noch nie den Fehler dabei gemacht, diese Meßpunkte mit Hilfe eines Lineals (schließlich will man ja sauber zeichnen und mit dem Kurvenlineal hat man zu wenig Übung) von Punkt zu Punkt zu verbinden. Dies ist zwar sauber, aber falsch. Da taucht natürlich die Frage auf, wenn nicht so, wie dann ? Welche Methode ist die richtige? Das kann eigent- lich nur dann beantworten werden, wenn der mathematische Zusammenhang der Meßergebnisse bekannt ist. Wie sind sie zustandegekommen? Welcher physikalische Zusammenhang ist gegeben. Sind die Meßergebnisse nahezu korrekt oder liegen Fehler (Streuungen) vor. Hat man es sogar mit Aus- reißern zu tun, die die Meßergebnisse total ver- fälscht haben, so daß sie nahezu unbrauchbar sind? Zur Analyse des des Problems werden verschiedene Methoden angeboten die zur Erstellung eines Datensatzes führen. Dieser wird dann wiederum, wie gehabt, durch das universelle Graphikprogramm GRAPHEN.PAS ausgewertet. -------------------- I Lineare Regression -------------------- Approximation durch eine Gerade (nach C.F. Gauß) Ist bei den Meßpunkten bekannt, daß die aufgenommene Funktion einer linearen Funktion der Form f(x) = mx + b entstammen, so sind bei einer optimalen Verbindung die beiden Koeffizienten m und b zu bestimmen. Geht man unbefangen an die Sache heran, so stellt man fest, daß es mehrere (genau : unendlich viele) Möglichkeiten gibt, eine Gerade durch die Reihe der Meßpunkte zu legen. Welche ist die optimale ? Hier hilft die Mathematik weiter. Die optimale Gerade erhält man durch folgende Definition : m und b sind so zu wählen, daß das Quadrat der Summe der Abstände zwischen den Meßwerten an den Stellen xi und den Funktionswerten f(xi) = mxi + b zu einem Minimum wird. Man spricht dabei von der Methode der kleinsten Fehlerquadrate. Mathematisch ausgedrückt : Σ (xi - f(xi))2 = Minimum Man kann zeigen, daß es nur eine Funktion f(x) mit der o.a. Eigenschaft gibt, diese also Eindeutigkeit ist. Doch wie sind das m und das b zu bestimmen ? Eine mathematische Herleitung ist nur mit Hilfe partieller Ableitungenden möglich. Diese ist in (fast) jedem Buch über numerische Mathematik bzw. Differenzialrechnung beschrieben. Eine gut lesbare Darstellung findet man in (1). Hier soll nur der Algorithmus zur Berechnung von m un b angegeben werden. 1. Man berechne aus den Meßwerten xi und yi die folgende Summen : Σ xi , Σ yi , Σ xi∙yi und Σ xi2 2. Aus Σ xi und Σ yi berechne man die Mittelwerte xm und ym indem man die Summen durch die Anzahl der Meßwerte dividiert. 3. Dann erhält man die Lösungen aus den Gleichungen : Σ xi∙yi - n ∙ xm ∙ ym i) m = ______________________ Σ xi2 - n ∙ xm2 ii) b = ym - m ∙ xm wobei das Summenzeichen Σ folgendermaßen definiert ist : Σ xi = x1 + x2 + ... + xn Beispiel : i x y x ∙ y x2 ____________________________________________ 1 1.0 4.2 4.2 1.0 2 2.0 8.5 17.0 4.0 3 3.0 12.7 38.1 9.0 4 4.0 17.2 68.8 16.0 ___________________________________________ Σ 10.0 42.6 128.1 30.0 10.0 42.6 ==> xm = ______── = 2.5 ym = _______ = 10.65 und damit 4 4 128.1 - 4 ∙ 2.5 ∙ 10.65 m = ________________________ = 4.32 30.0 - 4 ∙ 2.5 2 b = 10.65 - 4.32 ∙ 2.5 = - 0.15 Doch wie berechnet man eine solche Summe? Sicherlich ist diese Frage nur für den "Anfänger" interessant. Nachstehend der Algorithmus zur Berechnung einer beliebigen Summe in Form eines Struktogramms : ┌──────────────────────────────────────────────┐ │ Summe := 0 │ ├──────────────────────────────────────────────┤ │ Zähle mit i von 1 bis n │ │ ┌─────────────────────────────────────┤ │ und │ Berechnung Summand │ │ tue ├─────────────────────────────────────┤ │ │ Summe := Summe + Summand │ └────────┴─────────────────────────────────────┘ Die beiden Anweisungen Berechnung Summand und Summe := Summe + Summand können ggf. zu einer Anweisung zusammengefaßt werden. Das folgende Bei- spiel in Pascal zeigt die Berechnung der benötigten Summen : sumx := 0; sumy := 0; { Vorbesetzen der Summen } sumxy := 0; sumxx := 0; FOR i:=1 to n DO { Für jeden Meßpunkt } BEGIN xr := x[i]; { Zwischenwert, da einfacher } yr := y[i]; { zu schreiben } sumx := sumx + xr; { Berechnung der Summen } sumy := sumy + yr; { " } sumxy := sumxy + xr*yr; { " } sumxx := sumxx + sqr(xr) { " } END; Den vollständigen Algorithmus finden Sie auf der Sonderdiskette in der Funktion Berechne_m_b (File REGLIN.PAS) wieder, wobei diese noch zusätzlich die evtl. mögliche Division durch Null (Σ xi2 - n ∙ xm2 = 0) abprüft. Das ist z.B. dann der Fall, wenn nur ein Meßpunkt (evtl. auch mehrmals) einge- geben wurde. In diesem Fall ist keine eindeutige Gerade durch den singulären Punkt zu legen - die Funktion Lineare_Regression wird mit einer Feh- lermeldung abgebrochen. Doch wie gut ist eine solche Approximation, d.h. wie genau wird die berechnete Funktion den tatsächlichen Meßwerten gerecht ? Auskunft darüber gibt der sogenannte Korellationskoeffizient r : Σ (yi - ym)(xi - xm)2 r = _______________________________ √ [ Σ (yi - ym)2 Σ (xi - xm)2 ] den man nach der Berechnung von m und b bestimmen kann. Man kann zeigen: │r│ ≤ 1 Liegt der Betrag von r zwischen 0.4 und 0.6, so spricht man von einer mittleren, liegt er zwischen 0.6 und 0.8 von einer hohen und liegt er zwischen 0.8 und 1.0 von einer sehr hohen Korellation. Man beachte aber, daß ein sehr hoher Korellationskoeffizien nichts über einen etwaigen kausalen Zusammenhang zwischen der linearen Funktion und den Meßpunkten aussagen kann. Es gibt gute Gegenbeispiele dafür (Sogenannte Schein- oder Unsinnskorellationen z.B. der Zusammenhang zwischen der Kleidergröße und der Güte der Handschrift in den Klassen 1 bis 6). Diese sind gar nicht so selten und oft nur sehr schwer zu entlarven. Der Korrelationskoeffizient beweist gar nichts, er weist nur auf etwas hin, was einer genaueren Untersuchung bedarf. Im Programm wurde auf die Berechnung und Ausgabe von r verzichtet. An Hand des Graphens läßt sich beurteilen wie gut in ersten Näherung die Lineare Regression ist. Als aufmerksamer Leser müßten Sie in der Lage sein, diese Berechnung - wenn gewünscht - in das Programm einzufügen. Ist die Approximation Ihnen nicht gut genug, bietet wir Ihnen noch mehrere andere Möglichkeiten an, Ihre Meßwerte auszuwerten. ------------------- II Potenzfunktionen ------------------- Potenzfunktionen gehören zu der Sorte von Funktionen, die sich durch Umrechnen auf lineare Funktionen zurückführen lassen. Allgemein lautet die Gleichung einer Potenzfunktion : f : f(x) = b∙ax Durch beidseitiges Logarithmieren (hier zur Basis e) erhält man : f : ln (f(x)) = ln (b∙ax) = ln (b) + ln (ax) = ln (b) + x∙ln (a) Paßt man nun den Meßdatensatz entsprechend an, kann man ohne weiteres die bisher besprochenen Verfahren zur Auswertung von Daten heranziehen, von denen vermutet wird, daß deren mathematischer Zusammenhang mit einer Potenzfunktion beschrieben werden kann. Dazu berechnet man einen 2. Meßdatensatz durch : FOR i := 1 TO n DO BEGIN x[i] := Ln (x[i]); y[i] := Ln(y[i]); END; Mit diesem Satz führt man nun die lineare Regression durch und bekommt so die Werte von a und b. Bei der Tabelierung der Funktionswerte ist dann zu rechnen : y := b∙a^x Da PASCAL keinen Potenzoperator besitzt, muß dies umgerechnet werden in: y := b*EXP (x*Ln(a)); Das Verfahren ist nicht auf alle Potenzfunktionen anwendbar. Es gilt nur für x,y,a ε î+ ( x > 0, y > 0 und a > 0). ------------------------- III Exponentialfunktionen ------------------------- Auch Exponentialfunktionen der Form : ····························f··:··f(x)··=··a∙b1cx läßt sich zunächst umwandeln in :·····<=>······f··:··f(x)··=··a∙ecx∙ln(b) ···················<=>······f··:··f(x)··=··a∙ebx Diese Gleichung läßt sich durch logarithmieren in eine lineare Form umwandeln : ln (f(x)) = ln (a∙ebx) = ln (a) + ln (ebx) = ln (a) + b∙x Hier wird der 2 Meßdatensatz gewonnen durch : FOR i := 1 TO n DO y[i] := Ln (y[i]); Die Umrechnung fürs tabellieren geschieht dann durch die Anweisung : y := a*EXP(b*x); Auch hier gilt wieder, daß a,b,x ε î+ sein müssen. -------------------------- VI Die Logarithmusfunktion -------------------------- der Form f : f(x) = a + b∙ln(x) liegt schon linear vor. Hier muß nur noch jeder x - Wert durch logarithmieren umgewandelt werden und bei der Tabellierung die o.a. Funktionsvorschrift angewendet werden. Schon ist man fertig. Alle bisher besprochenen Verfahren lassen sich sowohl normal, d.h. nur mit Hilfe der Gaußschen Fehlerquadrate, als auch mit der im Heft be- schriebenen Methode des ROBUST FITTING anwenden. Wie gut dies funktioniert, zeigen die Abb. 8 und 9 wo jeweils einer Potenzfunktin und einer Exponentialfunktion "von Hand" ein einzelner Feh- ler aufgeprägt wurde. Man sieht deutlich, wie sich die Methode der Gaußquadratur aus der Bahn werfen läßt - also versagt - während ROBUST FITTING brauchbare Resultate liefert. Bei allen bisher beschriebenen Verfahren wurde der Meßdatensatz in zwei ARRAY's [1 .. Nmax] gehalten. Die Übergabe an Funktionen und Prozeduren erfolgte stets als CALL by Name. Dies hat einzig und allein den Grund, den Speicher klein zu halten, da Matrizen bekanntlich "Speicherfresser" sind. ------------------------- V Quadratische Regression ------------------------- Approximation durch eine Parabel (Funktion 2.ten Grades ) Genügt die lineare Approximation nicht, so kann man versuchen, die Meßpunkte durch eine Funktion 2.ten Grades zu approximieren. Bekannterweise ist diese gegeben durch : f(x) = ax2 + bx + c wobei jetzt die Koeffizienten a,b und c zu bestimmen sind. (Man beachte, dies ist keine Potenzfunktion, diese hätte die Gleichung f(x) = b∙x² !) Das ist natürlich nicht ganz so einfach wie vorhin, schließlich sind jetzt 3 Unbekannte zu berechnen. Geht man wieder nach der Methode der kleinsten Fehlerquadrate vor, so bekommt man fol- gendes lineares Gleichungssystem zur Bestimmung von a, b und c : a ∙ Σ xi2 + b ∙ Σ xi + c ∙ Σ 1 = Σ yi a ∙ Σ xi3 + b ∙ Σ xi2 + c ∙ Σ xi = Σ (xi ∙ yi) a ∙ Σ xi4 + b ∙ Σ xi3 + c ∙ Σ xi2 = Σ (yi ∙ xi2) (Man beachte : Σ 1 = n ) Die Koeffizienten dieses Gleichungssystems werden durch Summenbildung erzeugt. Das Gleichungssystem kann mit hinreichend bekannten Verfahren, wie z.B. der Cramerschen Regel (Determinantenverfahren) gelöst werden. Hier gibt es Fälle, wo durch unsinnige Eingaben das Lineare Gleichungssystem nicht eindeutig lösbar ist. Dann ist die Koeffi- zientendeterminate = 0 und die Berechnung wird üblicherweise mit einer Fehlermeldung abgebrochen. Hinweis : Die Berechnung der Ausgleichgeraden und der Ausgleichsparabel ist auch möglich bei mehrmaliger Eingabe des gleichen x - Wertes. Dies ist auch realistisch : bei einwandfreien Meßmethoden können durchaus die gleichen Werte mehrmals gemessen werden. Bedingt durch Meßungenauigkeiten etc. entstehen dabei in der Regel andere Meßergebnisse. Dies stört nicht weiter. Nur wenn nur ein Punkt (genauer nur ein x - Wert) eingegeben wurde, kann das Verfahren kein eindeutiges Ergebnis liefern. Die Folge wäre eine Division durch Null, welche im Programm durch ge- eignete Abfragen verhindert werden sollte. Ausgleichsfunktionen 3, 4, .. n-ten Grades lassen sich im Prinzip genau so aufstellen. Man kommt dann auf lineare Gleichungssysteme n-ten Grades. Deren Lösung macht im allgemeinen keine Schwierig- keit, jedoch haben die Funktionen die Eigenschaft, daß sie mitunter stark oszillieren (stark hin und her schwingen). Sie sind daher für eine Appriximation nicht so gut geeignet. ------------------ VI Interpolatition ------------------ - durch kubische Splines - Genügt auch die quadratische Regression nicht, so bleiben noch mehrere Wege offen : entweder man versucht eine Regression mit Funktionen noch hö- heren Grades bzw. nehme man einen anderen Funktionstyp oder man verbindet einfach die Punkte miteinander (Aber bitte nicht mit einem Lineal !). Dann erhält man eine Funktion f(x) mit der Eigenschaft f(xi) = yi, d.h. der entstehende Kurvenzug geht genau durch die Meßpunkte. In diesen Fällen spricht man von Interpolation. Dies könnte u.a. dadurch geschehen, daß man durch die n - Meßpunkte eine rationale Funktion vom Grade n - 1 legt. Bekanntlich ist das immer möglich : 2 Punkte lassen sich immer durch eine Gerade (= rationale Funktion 1.ten Grades) verbinden. Durch 3 Punkte läßt sich immer eine Parabel legen, bei 4 Punkten benötigt man eine rationale Funktion 3.ten Grades usw. usf.. Wo ist das Problem ? Geht man nach dieser Methode vor, so stört der Fundamentalsatz der Algebra ganz erheblich. Dieser besagt, daß eine rationale Funktion vom Grade n genau n Nullstellen besitzt. Dabei müssen diese nicht, können aber alle reell sein. Das bedeutet in der Praxis, daß die berechnete Funktion sehr oft die x - Achse schneidet. Der dabei entstehende Kurvenzug gibt dann nicht den erwarteten Funktionsverlauf wieder, da er stark oszilliert. Abhilfe schafft hier folgende Überlegung : Man interpoliere die gegebenen Meßpunkte mit Hilfe einer rationalen Funktion von möglichst niedrigem Grade intervallweise. Dabei sollen sich die Nahtstellen glatt aneinanderfügen, d.h. die Ableitungen der intervallweise definierten Funktionen müssen an den Schnittstellen gleichgroß sein (keine Ecken !). Weiter wird gefordert, daß die Länge der Kurve zu einem Minimum wird. Man sieht, ohne Kenntnisse der Differenzial- und Integralrechnung kommt man auch hier nicht weiter. Funktionen, welche die o.a. Forderungen erfüllen sind vom Grade 3 und heißen kubische Splines. Die mathematische Herleitung ist u.a. in (2) und (3) nachzulesen. Die Vorteile dieses Verfahrens sind : - Zur Berechnung sind kubische Splines wesentlich handlicher als rationale Funktionen höheren Grades. Dies gilt auch für die Ableitungen. - Kubische Splines geben im allgemeinen den Kurvenverlauf besser wieder, der entstehende Fehler wird kleiner. Sie wirken "glättend" haben "geringere Welligkeit" und "beruhigen" die Funktion. - Hervorragend geeignet sind sie, wenn die Meßpunkte in exakter Form vorliegen, wie das z.B bei durch n Punkte gegebenen Profilen der Fall ist. Ihren Namen haben diese Funktionen von dem gleichnamigen biegsamen Kurvenlineal (englisch : spline ). Algorithmus zur Berechnung der kubischen Splines : Die Funktion, welche wir suchen hat bekannterweise die Gleichung : f(x) = a + bx + cx2 + dx3 Einfacher wird das Verfahren jedoch, wenn wir statt dessen ansetzen : fi(x) = ai + bi∙(x - xi) + ci∙(x - xi)2 + di∙(x - xi)3 für jedes Intervall [xi,xi+1] und i = 0 .. n-1 (n Meßstellen) Vorgehensweise : 1) Man setze ai = yi···········für i = 0 .. n-1 2)··"····"···c0 = cn-1 = 0 3) Man löse das lineare Gleichungssystem zur Bestimmung der ci : 3 3 hi-1ci-1 + 2ci(hi-1 + hi) + hici+1 = -- (yi+1 - yi) - ---- (yi - yi-1) hi hi-1 für hi = xi+1 - xi und i = 1 .. n-2 ! 1 hi 4) bi = __ (yi+1 - yi) - __ (ci+1 + 2ci) für i = 0 .. n-2 hi 3 1 5) di = ____ (ci+1 - ci) für i = 0 .. n-2 3hi Bemerkungen : - Man kann zeigen, daß das lineare Gleichungssystem (3) eindeutig lösbar ist. Damit können dann die kubischen Splines berechnet wer- den. (2) bewirkt, daß die Splines für i = 0 und i = n-2 an der oberen bzw. unteren Intervallgrenze einen Wendepunkt besitzen, so daß sie über ihren Definitionsbereich hinaus linear fortgesetzt werden können. Die Realisierung des o.a. Algorithmus können Sie dem Quelltext (File SPLINE.PAS) auf der Sonderdsikette entnehmen. Zu beachten ist nur, daß die Meßwerte nicht wie o.a. von i = 0 .. n-1, sonder von i = 1 .. n in den beiden Arrays x und y gespeichert sind. Daher war eine Umrechnung der Indizes notwendig. Wegen der Klarheit des Pro- grammtextes wurde darauf verzichtet, die Berechnungen von i + 1 und i - 1 durch die beiden schnelleren Funktionen Succ (i) und Pred (i) vorzunehmen. Wer will, mags ändern. Aus dem gleichen Grund und wegen der Übertragbarkeit auf TURBO Pascal 3.0 wurden auch die Funktionen DEC und INC von TURBO PASCAL 4.0 nicht verwendet. Das ganze Verfahren kann nur angewandt werden, wenn die x - Werte in aufsteigender Reihenfolge vorliegen und kein x - Wert doppelt eingegeben wurde. Die Sortierung übernimmt die Function Bubblesort mit einem wenig effizienten aber dafür einfachen Algorithmus. In den meisten Fällen sind die Daten sortiert - dann ist Bubblesort schneller als hörere Verfahren (nur ein Durchlauf bei sortierten Daten). Sind x - Werte doppelt eingegeben worden, so können die Splines nicht berechnet werden, der Funktion Spline_Interpolation wird FALSE zuge- wiesen, die Berechnungen werden abgebrochen und das Programm gibt eine entsprechende Fehlermeldung aus. Das gleiche passiert, wenn die Anzahl der Punkte den Speicher sprengen würde. Immerhin müssen zur Berechnung des Gleichungssystems (3) ca. n² reelle Zahlen im Speicher gehalten werden. Da das Verfahren ursprünglich auf einem CPM - Rechner mit 64kB (ich kann mir heute nicht mehr vorstellen, wie man mit sowenig Speicherkapazität solche Programme schreiben kann) Hauptspeicher entwickelt wurde, sind die benötigten Matrizen dynamisch mit Hilfe von NEW und DISPOSE angelegt bzw. wieder entfernt worden. Wer will, mags ändern. Die Meßpunkt werden am Anfang des Programms mit einem vorangestellten "P" in die Datei REG.DAT geschrieben. GRAPHEN.PAS stellt diese Punkte im Koordinatensystem als normgerechte kleine Kreise dar. Beim Aufruf von REG.PAS meldet sich das Programm mit einer Auswahlmöglichkeit zur Gewinnung des Datensatzes. Zum einen kann eine Datei ausgewertet werden, die den gleichen Aufbau wie bei GRAPHEN.PAS - jedoch ohne Punkte - haben muß. Weiterhin können die Daten über die Tastatur eingegeben werden. Zum dritten ist es möglich, einen Testdatensatz mit Hilfe von Zufallszahlen zu erzeugen. Diese Verfahren ist ROBUSTLINFIT entnommen. Als besonderes Bonbon werden diesem Datensatz Ausreißer aufgeprägt. Mit Hilfe dieses Datensatzes kann man sehr anschaulich die Möglichkeiten von ROBUSTLINFIT (in diesem Heft) studieren. ROBUSTLINFIT ist in das Programm eingearbeitet worden. Ist man den Ausführungen gefolgt, dürfte das Nachvollziehen des Algorithmus keine großen Schwierigkeiten mehr bereiten. Zu bemerken ist, daß entgegen dem Handbuch von TURBO 4.0 (Band I S.251) Include - Dateien auch ohne vollständige Programmblöcke akzeptiert werden. Vermutlich müssen nur die Definitionen komplett abgeschlossen werden. Daher ist es nach wie vor möglich, ein Programm wie üblich aufzuteilen. Abb. 7 zeigt die lineare und quadratische Approximation eines Testdatensatzes und die Interpolation mit kubischen Splines. Man sieht sehr deutlich, daß die kubischen Splines im Gegensatz zu den Approximationen exakt durch die Meßpunkte gehen. Ich hoffe, daß Ihnen das Lesen etwas Spaß bereitet hat, auch wenn die Materie dem einen oder anderen sicherlich zu mathematisch ist. Vernünftige Programme bekommt man aber nur durch gründliche Problemanalyse - und dabei benötigt man eben mathematische Methoden. Ausprobieren ist eben doch nicht alles. Literatur : (1) W.Kroll : Grund- und Leistungskurs ANALYSIS Band 1 : Differenzialrechnung 1 Dümmlers Verlag ISBN 3-427-42811-7 Bonn 1975 (2) Dr. Hans Ade, Dr. Hugo Schell : Themenhefte Mathematik Numerische Mathematik Ernst Klett Verlag ISBN 3-12-739800-X Stuttgard 1975 (3) Prof. Dr. M.W. Müller : Höhere Mathematik IV Vorlesung gehalten im SS 1975 (nicht veröffentlicht) (4) Hans Jochen Bartsch : Taschenbuch mathematischer Formeln Verlag Harri Deutsch ISBN 3871442399 Frankf./Main 1975 ---------------------- Numerische Integration ---------------------- - Algorithmen zur computergestützten Berechnung von Integralen - Bei vielen Problemen in Wissenschaft und Technik gilt es Integrale von zu berechnen. Oftmals ist diesen Integralen mit dem erlernten mathematischen Handwerkszeug nicht beizukommen, eine geschlossene Auswertung ist entweder unmöglich oder zu aufwendig. In solchen Fällen ist man auf eine numerischen Bearbeitung des Problems angewiesen. Während sich die Wissenschaftler früherer Zeiten noch mit Rechenschieber und Tabellen abquälen mußten, steht ihren heutigen Kollegen ein wesentlich leistungsfähigeres Arbeitsmittel zur Verfügung: der Computer. ------------------- Der Integralbegriff ------------------- Wie die Differentialrechnung, so geht auch die Integralrechnung auf ein geometrisches Problem zurück, nämlich die Berechnung von Flächeninhal- ten. In einem x-y-Koordinatensystem bezeichnet man mit ⌠b A =···│ f(x) dx ⌡a (sprich: "Integral f(x) dx von a bis b") den von dem Funktionsgraphen y=f(x), der x-Achse und den Geraden x=a und x=b eingeschlossenen Flächeninhalt (Abb. 1). Da dieses Problem systematisch erstmals Mitte des 19. Jahrhunderts von dem deutschen Mathematiker Bernhard Riemann untersucht wurde, spricht man in Würdigung seiner Arbeiten auch häufig vom "Riemann'schen Integral". Neben diesem gibt es in der Mathematik noch weitere Integralbegriffe, welche jedoch meist nur von theoretischem Interesse sind. Wir wollen uns hier deshalb ausschließlich mit dem Riemann- Integral beschäftigen. --------------------------------------------------------- Trapezregel - ein erstes, einfaches Integrationsverfahren --------------------------------------------------------- Es gibt eine Reihe von Möglichkeiten, den Flächeninhalt unter einem Funktionsgraphen zu bestimmen. Die hieraus abgeleiteten Rechenvor- schriften und -regeln bezeichnet man als Quadraturformeln. Die wohl einfachste Formel dieser Art ist die sogenannte "Trapezregel". Die Eckpunkte [a,f(a)] und [b,f(b)] werden durch eine Gerade verbunden, die zu berechnende Fläche wird damit zu einem Trapez, dessen Inhalt nun einfach anzugeben ist (Abb.2, [1]). Wie man leicht sieht, stellt diese Annäherung nicht gerade das Nonplusultra an Genauigkeit dar, es besteht eine je nach "Glattheit" der Funktion mehr oder weniger große Abweichung vom tatsächlichen Wert des Integrals. Dies kann verbessert werden, indem man das zu integrierende Intervall in kleinere Teilintervalle unterteilt und auf jedes dieser Teilintervalle die Trapezregel anwendet. Die Trapezsehnen passen sich dann besser dem tatsächlichen Funktionsverlauf an; Aufsummieren der einzelnen Trapeze ergibt die sogenannte große Trapezregel (Abb. 2, [2]). ---------------------------- Darf es etwas genauer sein ? - Simpson's Formel ---------------------------- Ein Maß für die Güte einer Quadraturformel ist die sogenannte "Fehlerordnung O". Sie gibt für eine Unterteilung des gesamten Intervalls (a,b) in N Teilintervalle der Länge h = (b-a)/N an, zu wel- cher Potenz von h die Abweichung der gemachten Näherung vom tatsächlichen Wert proportional ist. Die Trapezregel besitzt zum Beispiel die Fehler- ordnung O(h2), der absolute Fehler einer Trapez- Näherung ist also proportional zum Quadrat von h. Eine Methode mit einer wesentlich besseren Fehlerordnung als die der Trapezregel ist die ebenfalls recht bekannte Simpson-Regel, auch Kep- lersche Faßregel genannt. Sie beruht auf einer Annäherung des Kurvenverlaufs einer Funktion f durch eine Parabel durch die Punkte [a,f(a)], [m,f(m)] und [b,f(b)], wobei der Parameter m = (b- a)/2 die Intervallmitte bezeichnet (Abb.3, [1]). Da eine Parabel zur Annäherung einer Funktion anschaulich gesprochen "biegsamer" als eine Gerade (beim Trapezverfahren) ist, lassen sich hier natürlich schon wesentlich bessere Näherungswerte für den Flächeninhalt erzielen. Mit einer Intervallunterteilung wie beim Trapezverfahren beschrieben gelangt man zur großen Simpson-Formel (Abb.3, [2]) mit der Fehlerordnung O(h4). Allerdings ist hier zu beachten, daß die große Simpson-Formel nur auf Intervalle mit ge- radzahliger Unterteilung anwendbar ist. ------------------- Gauß'sche Quadratur ------------------- Eine weitere Möglichkeit zur näherungsweisen Integralberechnung sind die sogenannten Gauß'schen Quadraturformeln. Die Ableitung dieser Rechenregln soll hier nicht vorgeführt werden, sie erfordert Einiges an Mathematik. Das Verfahren an sich soll jedoch nicht unerwähnt bleiben, da es gerade für den Einsatz in Verbindung mit elektronischen Datenverarbeitungsanlagen sehr geeignet ist. Die Grundidee der Gauß-Quadratur ist ähnlich der bei Trapez- und Simpson- Verfahren. Man konstruiert Regeln, die Polynome eines bestimmten Grades exakt integrieren (Trapez: Polynom 1.Grades (=Gerade), Simpson: Polynom 2.Grades (=Parabel)) und wendet diese auf die Funktion f an. Da sich jede Funktion lokal durch ein Polynom annähern läßt (je höher der Grad des Polynoms, umso besser ist meist die Annäherung), wird dieses Verfahren natürlich um so genauer, je höher der Grad der exakt zu integrierenden Polynome ist. Der Unterschied zu den bereits beschriebenen Verfahren ist nun, daß bei der Konstruktion der Formeln weder Stützstellen (die x-Werte, an denen ein Funktionswert ermittelt wird) noch Gewichte (die Zahlen, mit denen die einzelnen Funktionswerte multipliziert werden) von vorn herein festgelegt werden. Dies hat zur Folge, daß man Näherungsformeln noch höherer Genauigkeit, jedoch mit "krummen" Stützstellen und Gewichten erhält. Das ist auch der Grund dafür, warum die Gauß-Quadratur gerade für elektronische Rechenmaschinen geeignet ist: das Verfahren wird unrentabel, wenn man keine Möglichkeit hat, die Stützstellen und Gewichte dauerhaft abzuspeichern, um sie bei Bedarf in die Rechnung miteinzube- ziehen. Eine Gauß'sche Quadraturformel vom Grad n integriert Polynome vom Grad 2n+1 exakt, die Fehlerordnung der Formel ist O(h2n+1). Das neben den Prozeduren "Trapezoid" und "Simpson" ebenfalls in der Turbo-Pascal-Unit NUMINT.PAS implementierte Unterprogramm "Gauss" arbeitet mit einer Formel vom Grad 15, besitzt also die Fehlerordnung O(h31). ----------------------------- Ein erstes Anwendungsbeispiel ----------------------------- Ein Schraubenhersteller möchte gerne wissen, mit wieviel Prozent Produktionsausschuß er rechnen muß, wenn er seiner neusten Sorte Präzi- sionsschrauben einen Durchmesser von 14.95 bis 15.05 Milllimeter garantiert. Aus umfangreichen Messungen ist bekannt, daß der der mittlere Durchmesser dieser Schrauben µ=15 mm bei einer Streuung von σ=0.02 mm beträgt. Solche unvermeidlichen produktionsbedingten Abweichungen vom Sollwert entstehen zumeist durch zufällige Fehler und sind deshalb normalverteilt, d.h. ihre Verteilung (die Wahrscheinlichkeit für einen bestimmten Wert) wird durch die Gauß'sche Glockenkurve (Abb.4) beschrieben: 1··········┌ -(x-µ)2 ┐ p(x) = ───────── ·exp │─────────│ σ·(2π)½·······└···2σ2···┘ Den gesuchten Wert (besser gesagt das Komplement dazu) findet der Schraubenfabrikant nun, indem er alle Wahrscheinlichkeiten im Bereich von µ-δ bis µ+δ aufsummiert (δ=0.05), d.h. das Integral ⌠µ+δ P =···│ p(x) dx ⌡µ-δ berechnet. Da die Gauß-Funktion eine der Funktionen ist, die nicht geschlossen (also mit Papier und Bleistift) integrierbar sind, ist man bei der Berechnung auf eine numerische Auswertung angewiesen. Wendet man hierzu das Trapez- Verfahren an, kann man ausnutzen, daß der maximale Fehler f dieses Verfahrens durch die Fehlerordnung und das Maximum der 2.Ableitung der zu integrierenden Funktion im Integrationsintervall gegeben ist: f ≤ h2/12 · max (f"(x)) = g Aus dieser Gleichung kann für eine gegebene Genauigkeit g der Wert für h = (b-a)/N, damit also auch die Anzahl N der nötigen Intervallunter- teilungen ermittelt werden. Unter Benutzung des Programms PROBABIL.PAS und Eingabe der oben angegeben Wert ergibt sich das Ablaufprotokoll in Abb.5. Wie man sieht, besitzen 98.8% aller produzierten Schrauben den garantierten Durchmesser, es muß also mit 1.2% Ausschuß gerechnet werden. ------------------------------------------------------- Automatische Fehlerkontrolle durch Halbierungsverfahren ------------------------------------------------------- Die im letzten Beispiel verwendete Genauigkeitskontrolle ist zwar recht praktisch, doch oft ist man nicht im Stande bzw. es ist zu unwirtschaftlich, die Ableitungsmaxima zu berechnen (bei der Gauß-Integration müßte man beispielsweise das Maximum der 30. Ableitung von f berechnen, viel Spaß!). Aus diesem Grunde bedient man sich zur Genauigkeitskontrolle einer Methode, welche zwar etwas zeitintensiver ist, dafür dem Benutzer aber sämtliche Arbeit abnimmt: bei der Auswertung wird nicht nur eine Näherung, sondern durch sukzessive Unterteilung des Integrationsintervalls eine Folge von Näherungswerten berechnet. Sobald die Differenz zweier aufeinanderfolgender Näherungswerte die geforderte Genauigkeit unter- schreitet, kann man davon ausgehen, daß der letzte Näherungswert die geforderte Genauigkeit aufweist. Dieses Verfahren wollen wir sofort an einem Beispiel ausprobieren: der Umfang einer Ellipse soll (im Rahmen der Rechengenauigkeit eines Compu- ters) exakt berechnet werden. Nach einer kleinen Rechnung findet man, daß der exakte Umfang durch das Integral ⌠π ┌·····b2-a2·······┐½ U = 2b···│··│1 - ─────── sin2Θ│ dΘ ⌡0 └······b2·········┘ gegeben ist. Die Parameter "a" und "b" sind hierbei die Halbachsen der Ellipse (Abb.6). Das Integral ist ein sogenanntes elliptisches Inte- gral, es kann wie die Gauß-Funktion nicht geschlossen gelöst werden (aus diesem Grund findet man in Formelsammlungen auch nur Näherungsformeln für den Ellipsenumfang). Mit dem Programm ELLIPSE.PAS kann dem Problem auf den Leib gerückt werden: es berechnet den Umfang einer Ellipse mit Hilfe der Simpson-Regel und oben beschriebenem Kontrollverfahren, welches als Unterprogramm "Successive" in der NUMINT-Unit realisiert ist. ------------------------------- Adaptive Verfahren - Integration für Feinschmecker ------------------------------- Nichts ist so gut, daß es nicht noch zu verbessern wäre: auch die oben vorgestellte Methode hat einen Schwachpunkt. Naturgemäß sind Quadraturformeln dort ungenau, wo die zu integrierende Funktion starke Anstiege bzw. Abfälle besitzt. Integriert man nun eine Funktion, die an einer Stelle im Integrationsintervall einen solchen Peak besitzt, ansonsten aber recht glatt verläuft, so muß der Peak- Bereich in sehr kleine Intervalle aufgeteilt werden, damit eine gute Näherung möglich ist. Unser obiges Kontrollprogramm sieht jedoch eine Gleichbehandlung aller Bereiche des Integrationsintervalls vor, deshalb werden auch die glatten Teilstücke sehr fein aufgeteilt. Dies hat nun aber zur Folge, daß viele Funktionswerte umsonst berechnet werden, da an den glatten Stellen ja schon wenige Stützwerte zur erforderlichen Genauigkeit ausreichen. In den meisten Fällen führt das soweit, daß Teilstücke des Integrationsintervalls, die schon mit der ma- ximal überhaupt erreichbaren Rechengenauigkeit be- stimmt sind, immer noch weiter aufgespalten werden. Hierdurch verschenkt man natürlich jede Menge Rechenzeit, manche Integrale lassen sich so überhaupt nicht bestimmen. Abhilfe schafft hier ein "intelligentes" Steuerprogramm, welches sich die Bereiche merkt, in denen Schwierigkeiten enstehen und nur diese Intervalle besonders genau behandelt. Das intelligente Unterprogramm, es handelt sich um die Prozedur "Adaptive", welche sich ebenfalls in NUMINT.PAS befindet, läßt sich sehr schön rekursiv formulieren, der Rechner nimmt dem Programmierer so praktisch alle Verwaltungsarbeiten ab. Auch hierzu ein Anwendungsbeispiel: die Wurzelfunktion y = x½ soll im Intervall von 0 bis 1 integriert werden. Hierbei gibt es im Bereich nahe Null Probleme, denn dort wächst die Steigung der Wurzelfunktion gegen Unendlich (Abb.7). Das Programm SQROOT.PAS zeigt die numerische Behandlung des Wurzelintegrals mit dem oben erwähnten rekursiven Steuerprogramm. Der Leser möge mit diesem Programm einmal den Integralwert zwischen a=0 und b=1 berechnen (der exakte Wert ist 2/3) und bei fester Genauigkeit den Re- chenzeitunterschied zur Berechnung des Integrals mit "Successive" vergleichen! ----------------------------- Warum verschiedene Verfahren? ----------------------------- Nachdem wir jetzt einige Methoden und Verfahren kennengelernt haben, Integrale numerisch mit dem Computer zu "knacken", mag sich der ein oder andere Leser fragen, wieso hier so viele verschiedene Unterprogramme vorgestellt werden, letztendlich nimmt man ja doch das schnellste und leistungsfähigste. Doch gerade diese beiden Be- griffe können in verschiedenen Situationen auch völlig verschiedene Bedeutungen haben, bei der Auswahl des Verfahrens kommt es ganz auf die Erfordernisse an. Ein Beispiel soll dies verdeutlichen: die Integralfunktion ⌠∞ Γ(α) = │ exp(αx) dx ⌡0 soll auf dem Bildschirm in einem Γ-α-Koordinatensystem grafisch dargestellt werden. Für die Berechnung der vielen Integrale ist es sicherlich unzweckmäßig, das Gauß- Verfahren anzuwenden. Dieses zeigt seine Stärken erst bei höheren Genauigkeitsforderungen; da für die Ermittlung von Grafik-Koordinaten aber eine Genauigkeit von 4 Stellen schon vollkommen ausreicht, wird man sich lieber ein weniger genaueres, dafür aber sehr schnelles Verfahren (z.B. Simpson oder Trapez) herausgreifen, man will ja schließlich auch nicht stundenlang auf einen fertig aufgebauten Bildschirm warten. In anderen Anwendungen hingegegen wird wieder mehr die Genauigkeit gefragt sein, die Rechenzeit spielt dann keine so große Rolle. Unter den vorgestellten Programmen sollte eigentlich für jeden Zweck die passende Routine vorhanden sein, die Quadraturformeln ergeben in Verbindung mit den Kontrollprogrammen eine Vielfalt an Einsatzmöglichkeiten. ------------------------ Oberflächen und Volumina ------------------------ Es ist nun nicht so, daß die Integrationsprogramme lediglich Flächen unter irgendwelchen Funktionen berechnen können, mit ein Paar Handgriffen kann man sie auch dazu bewegen, Oberflächen und Volumi- na von dreidimensionalen Körpern zu berechnen. Vorraussetzung ist allerdings, daß diese Körper rotationssymmetrisch sind, d.h. sie werden von einer um eine Achse rotierenden Funktion f erzeugt. Beispiele für solche "Drehkörper": y = k (Zylinder) y = k·x (Kegel) y = k·x½ (Paraboloid) Abb.8 zeigt diese Funktionen, die bei Rotation um die x-Achse den entsprechenen Körper erzeugen. Der Parameter "k" soll eine beliebige reelle Zahl repräsentieren. Nach der sogenannten Guldin'schen Regel gilt für Volumen und Oberfläche (Mantel) solcher Körper ⌠b M = ··2π·│ y · (1 + y'2)½ dx ⌡a ⌠b V =··· π·│ y2 dx ⌡a Das Programm GULDIN.PAS berechnet nach diesen Formeln Volumen V und Oberfläche M eines Drehkörpers. ---------------------------- Mehrdimensionale Integration ---------------------------- Zum Abschluß wollen wir noch einen kurzen Abstecher in das überaus in- teressante Gebiet der Integration in mehreren Dimensionen machen. Analog zu der bisher besprochenen eindimensionalen Integration gibt es für Funktionen mit mehreren Variablen mehrdimensionale Integration, jede Variable entspricht einer Dimension. Geometrisch veranschaulicht werden hier nicht nur Flächen, sondern auch Oberflächen und Volumina (jedoch von beliebigen Körpern!) berechnet. Ein Beispiel für ein sol- ches Integral: ⌠π···⌠π···⌠π V =···│ dx │ dy │ dz sin (x+y+z) ⌡0···⌡0···⌡0 Im mehrdimensionalen Fall schreibt man oft die Differentiale (dx,dy,dz) direkt hinter das Integralzeichen, um die Zugehörigkeit von Integra- tionsvariable und Grenzen deutlich zu machen. Die mehrdimensionale Integration ist für den Computer ein ganz heißes Eisen: alle bisher in einer Dimension aufgetauchten Probleme vervielfa- chen sich hier mit der Potenz der Dimension. So hat man beispielsweise für ein dreidimensionales Integral (welches in allen drei Integrations- variablen numerisch gleich schwer zu behandeln sein soll), bei dem für eine Dimension n Funktionswerte berechnet werden müssen, letztlich we- nigstens n3 Funktionswerte zu berechnen. Zeitlich gesehen heißt das: benötigt der Rechner für eine Dimension beispielsweise 30 Sekunden, so wird man auf das Ergebnis in drei Dimensionen ganze 303 = 27000 Sekun- den, das sind 7½ Stunden warten müssen! Wie man sieht, sind in der Re- gel bei solchen Problemen, falls man sie wirklich numerisch angehen muß, extrem leistungsfähige Rechner gefragt. Beschränkt man sich jedoch auf Integrale mit maximal drei Dimensionen und gibt sich mit Genauigkeiten von wenigen Stellen (allgemein nicht mehr als fünf) zufrieden, so kann man auch auf seinem Hobby-Rechner schon ganz brauchbare Resultate erzielen. Das Unterprogramm "Romberg" in NUMINT.PAS implementiert ein relativ simples Verfahren, um mehrdimensionalen Integralen numerisch beizukom- men. Um die Funktionstüchtigkeit dieser Routine auszutesten, soll das Volumen einer Kugel berechnet werden. Dies ist zwar kein Paradebeispiel für den Einsatz des Programms - den Wert kann man in jeder Formelsamm- lung nachschlagen, falls man ihn ohnehin nicht schon aus dem Kopf kennt - doch läßt sich hiermit sehr schön der Gebrauch von "Romberg" demon- strieren. In sphärischen Polarkoordinaten ist das Volumen einer Kugel einfach das Integral über das Volumenelement in den Begrenzungen des Körpers: ⌠2π··⌠π···⌠R V =···│ dΦ │ dΘ │ dr r2·sinΘ ⌡0···⌡0···⌡0 Das Programm SPHERE.PAS berechnet dieses Integral, nach Eingabe des gewünschten Kugelradius wird das aus der Integration sowie zum Ver- gleich das aus der bekannten Formel V = 4/3·πR3 bestimmte Kugelvolumen angezeigt. Waren bisher alle Integrationsunterprogramme vom Aufruf her mehr oder weniger gleich aufgebaut, so gibt es hier nun einige gravierende Unter- schiede: die zu integrierende Funktion hat als Argument nicht mehr nur eine einzige REAL-Zahl "x", sondern einen Vektor von REAL-Zahlen. Für unser Kugel- Beispiel könnte das so aussehen, das x[1]=r, x[2]=Θ, x[3]=Φ ist. Entsprechend sind die Integrationsgrenzen keine einfachen REAL-Zahlen mehr, sondern ebenfalls Vektoren. Bei der Besetzung der Felder muß man natürlich darauf achten, daß die Indizes nicht durchein- ander wirft, also a[1],b[1] gehört zu x[1] usw. Desweiteren hat man die Möglichkeit, eine relative Unterteilung des Integrationsraumes in Bezug auf die verschiedenen Dimensionen vorzunehmen. Das ist sehr nützlich, wenn man ein Integral bearbeitet, welches in einer Richtung sehr schwierig zu integrieren ist. Man kann diese dann durch eine entspre- chende Unterteilung besonders genau behandeln, ohne das gleiche in den anderen Dimensionen auch tun zu müssen, wertvolle Rechenzeit wird ein- gespart. Literatur Wer bei diesem "Einstieg" in das Gebiet der Numerik mit Computern auf den Geschmack gekommen ist oder wem die hier angegebenen Erklärungen nicht ausreichen, dem sei folgende weiterführende Literatur empfohlen: Carsten Gieselmann (1) Engeln-Müllges/Reutter Numerische Mathematik für Ingenieure (2) Stoer Einführung in die Numerische Mathematik (3) Davis/Rabinowitz Methods of Numerical Integration ------------------------------------------------------------- Iterative Berechnung eines Linearen Gleichungssystems mit dem Einzelschritt - Verfahren von Gauss - Seidel. ------------------------------------------------------------- Die Berechnung linearer Gleichungssysteme bereitet in der Regel keine Schwierigkeiten, ist der dazugehörende Gauß'sche Algorithmus oft genug publiziert worden. So auch in der Zeitschrift PASCAL (Heft ???). Nun ist aber die Berechnung größerer Systeme recht problematisch und zwar in zweierlei Hinsicht. Erstens wächst mit steigendem n (n = Anzahl der Gleichungen) der Zeitbedarf proportional zu n3 (!). Dies fällt sicher nicht in's Gewicht, wenn man nur hin und wieder ein solches System zu lösen hat. Zweitens aber - und das ist viel wichtiger -, werden auch die Rundungsfehler beim Lösen mit dem Gauß'schen Algorithmus mit stei- gendem n immer größer. Das kann schließlich soweit führen, daß die Lö- sungen praktisch nicht mehr zu gebrauchen sind. Findige Mathematiker sind dieser Frage nachgegangen und haben schon sehr früh ein Verfahren entwickelt, welches diese beiden Nachteile nicht hat. Dafür hat dieser Algorithmus - wie könnte es anders sein, nichts auf der Welt ist Ideal - einen anderen Pferdefuß. Doch davon unten später mehr. ---------------------------- Der Banarsche Fixpunktsatz : ---------------------------- Bevor wir uns dem Verfahren zuwenden, müssen einige mathematische Be- griffe klargestellt werden. Aber keine Angst, ganz trocken und hochtra- bend wird die Sache hier nicht behandelt. Für eine exakte mathematische Beschreibung gibt es schließlich Bücher. Nehmen wir einmal an, wir hätten eine Funktion y = f(x), welche wir nicht nach x auflösen können (oder wollen). Beispielsweise die Funktion f(x) = x4 - 9x2 + 4x + 12 . Gesucht ist eine Nullstelle dieser Funk- tion, also f(x) = 0. Hat man überhaupt keine Ahnung vom Lösen von Glei- chungen und schreckt unbefangen vor nichts zurück, so kann man die Gleichung umformen in : 9x2 - 4x - 12 x = ─────────────── x3 Hat man nun eine Lösung, so müssen linke und rechte Seite übereinstim- men. Hat man keine Lösung, so kann man durch probieren versuchen, die Lösung zu finden. Wir setzen z.B. ein x = 3 und berechnen : 9∙32 - 4∙3 - 12 x = ──────────────── = 2.1111111111 33 Stimmt also nicht. Was hindert uns aber, den so errechneten Wert wieder in die Formel einzusetzen ? Dann erhalten wir x2 = 2.0902 ... x3 = 2.0762 ... x4 = 2.0660 ... . . . . x100 = 2.0050 ... Wir vermuten schließlich die Lösung bei x = 2, was wir durch Einsetzen bestätigt finden. Selbstverständlich rechnen wir dies nicht zu Fuß aus, sondern schreiben dazu ein kleines BASIC - Programm (warum eigentlich nicht, BASIC - Programme bis zu 10 Zeilen Länge haben schließlich ihre Vorteile). 10 X = 3 20 FOR I=1 TO 100 30 X = (9*X*X-4*X-12)/(X*X*X) 40 NEXT I 50 PRINT "X = ";X 60 END Die Lösung geschieht also in 2 Schritten : 1. Umwandeln der Gleichung in die Form x = f(x). 2. Erzeugen der Folge··················xi = f(xi-1) mit einem beliebigen Startwert x0 (in unserem Beispiel x0 = 3) Wenn der Fehler der so gefundenen Lösung beliebig klein gehalten kann, indem man das Verfahren genügend oft wiederholt, spricht der Mathemati- ker von Konvergenz. Es ist hoffentlich klar, daß das Verfahren seine Tücken hat. Die Kon- vergenz ist nicht bei jeder Umformung und nicht für jeden beliebigen Startwert x0 gewährleistet! Probieren Sie mal verschiedene Werte und Gleichungen aus. Der Banarch'sche Fixpunktsatz macht nun Aussagen darüber, wann eine solche Folge konvergiert, d.h. unter welchen Voraussetzungen nach end- lich vielen Schritten ein brauchbarer Wert für x gefunden wird. ----------------------------------------- Anwendung auf lineare Gleichungssysteme : ----------------------------------------- Wir wollen nun unser erworbenes Wissen auf ein lineares Gleichungssy- stem anwenden. Beispiel : 4x - y + z = 4 x - 2y = -1 3x + 4z = 7 Dazu stellen wir die 1. Gleichung nach x, die 2. nach y und die 3. nach z um. Wir erhalten : x = (4 + y - z)/4 y = (-1 - x)/(-2) z = (7 - 3x)/4 Beginnen wir mit einem beliebigen Startwert z.B. x = 4, y = -2 und z = 4, so erhalten wir nach je einem Schritt : ( 1) x = -0.5000 y = 0.2500 z = 2.1250 ( 2) x = 0.5312 y = 0.7656 z = 1.3515 ( 3) x = 0.8535 y = 0.9267 z = 1.1098 . . . . . . (10) x = 0.9999 y = 0.9999 z = 1.0000 Offensichtlich hat das Gleichungssystem die Lösung x = 1, y = 1 und z = 1, was sofort durch Einsetzen bestätigt wird. Das dazugehörende BASIC Programm : 10 X=1:Y=1:Z=1 20 FOR I=1 TO 10 30 X=(4+Y-Z)/4 40 Y=(-1-X)/(-2) 50 Z=(7-3*X)/4 60 PRINT "X= ";X;" Y= ";Y;" Z= ";Z 70 NEXT I 80 END Damit dürfte das Einzelschrittverfahren klar sein ! Bleibt zu untersu- chen, unter welchen Umständen das Verfahren konvergiert, also die ge- wünschte Lösung liefert. Die Anwendung des Banarch'schen Fixpunktsatzes liefert als Ergebnis das Zeilen- und Spaltensummenkriterium : Das Einzelschrittverfahren konvergiert genau dann, wenn 1. Das Gleichungssystem eine eindeutige Lösung besitzt und 2. Entweder das Zeilensummenkriterium oder das Spaltensummenkriterium erfüllt sind. Zeilensummenkriterium : Für jede (!) Zeile i (i = 1 .. n) der Matrix A muß gelten, daß die Summe der Beträge der Koeffizienten kleiner ist, als das doppelte Produkt des Elementes aii. Mit mathematischen Symbolen geschrieben : n Σ │aik│ < 2∙│aii│ (i = 1 .. n) k=1 Entsprechendes gilt für das Spaltensummenkriterium. ┌ ┐ │4 -1 1│ In unserem Beispiel ist A = │1 -2 0│ (Koeffizientenmatrix) │3 0 4│ └ ┘ => : 1. Zeile : 4 + 1 + 1 < 2∙4 o.k. 2. " : 1 + 2 + 0 < 2∙2 o.k. 3. " : 3 + 0 + 4 < 2∙4 o.k. In diesem Gleichungssystem ist also das Zeilensummenkriterium erfüllt. (Nebenbei, auch das Spaltensummenkriterium. Der Leser möge dieses durch Rechnung nachprüfen !). Damit konvergiert das Verfahren für jeden be- liebigen Startwert (x,y,z). Leider sagen diese Kriterien nichts darüber aus, wann ein Verfahren nicht konvergent ist! Das heißt, liefern Zeilen- oder Spaltensummenkri- terium keinen Beweis dafür, daß das Verfahren konvergiert, so heißt das noch lange nicht, daß das Verfahren keine brauchbaren Ergebnisse lie- fert. Auch kann man manchmal mit dem Vertauschen von Zeilen erreichen, daß ein System noch eines der Kriterien erfüllt, obschon dies vor der Vertauschen nicht der Fall war. Im Programm wurde darauf allerdings verzichtet. -------------- Das Programm : -------------- Um nicht jedesmal bei einem gegebenen Gleichungssystem ein neues Pro- gramm schreiben zu müssen, wurde das Program LINIT (LINeares Glei- chungssystem ITerativ) geschrieben, welches sich auf der Sonderdiskette befindet. Hier konnte die Zielsprache nur eine Sprache sein, die struk- toriertes Programmieren zuläßt. Primitive BASIC - Dialekte scheiden daher aus. Die Sprache PASCAL wurde ausgewählt. Um eine möglichst leichte Übertragbarkeit zu gewährleisten, wurde das Programm so geschrieben, daß es sowohl unter TURBO - PASCAL 3.0 wie auch unter TURBO - PASCAL 4.0 compeliert werden kann. Die dazu minima- len notwendigen Änderungen sind im Quell - Text angegeben. Matritzen sind Speicherfresser. Daher wurden alle Matritzen in den Pro- zeduren und Funktion als Call by Value übergeben (VAR !). Nur so war es mir möglich, auch Gleichungssysteme mit großem n (getestet wurde bis n = 31) direkt in den Speicher (640k) zu compelieren. --------------------------------------------------------- Die wichtigsten Funktionen und Proceduren : --------------------------------------------------------- ------------------------ I. FUNCTION Konvergent : ------------------------ Diese Funktion überprüft wie oben angegeben das Zeilen- und Spaltensum- menkriterium. Ist keines der beiden Kriterien erfüllt, so wird der Be- nutzer nach einer Warnung aufgefordert einzugeben, ob er den Versuch der Lösung trotzdem waagen will. Entsprechend seiner Eingabe wird Kon- vergent auf TRUE oder FALSE gesetzt. ----------------------- II. FUCTION Maxfehler : ----------------------- Die Beurteilung der Genauigkeit einer Lösung eines linearen Gleichungs- systems ist außerordendlich schwierig. Hier wird man, je nach Problem, eine entsprechendes Anpassung vornehmen müssen. -------------------------- III. PROCEDURE Iteration : -------------------------- Zwecks schnellerer Verarbeitung durch den Rechner wurde der Iterations- schritt folgendermaßen umgestellt : Für jede Zeile i = 1 .. n gilt : n Σ aikxk = bi k=1 n <=> aiixi + Σ aikxk - aiixi = bi k=1 n <=> aiixi = bi - (aiixi + Σ aikxk) k=1 n bi + aiixi - Σ aikxk k=1 <=> xi = ────────────────────── aii ··············n Das··Produkt Σ aikxk wird durch die Funktion Produktsumme ·············k=1 ermittelt. ------------------- IV PROCEDURE Eingabe : ------------------- Die Eingabe der Koeffizienten des Gleichungssystems kann sowohl über die Tastatur (umständlich und nur für Testzwecke geeignet) als auch mit Hilfe einer Textdatei erfolgen. Diese hat dann folgenden Aufbau : (vergl. dazu das obige Beispiel) 3 4 -1 1 4 1 - 2 0 -1 3 0 4 7 Die Datei TEST.lin enthält ein Gleichungssystem mit 31 Unbekannten zum ausprobieren. Die Arbeitsweisen der anderen Proceduren und Funktionen gehen aus dem Quelltext hervor und bedürfen sicherlich keinerlei Erklärungen mehr. -------------------------------------- Anwendungen und Grenzen des Verfahrens. -------------------------------------- Die Zeit, welche der Gauß'sche Algorithmus benötigt, ist wie oben be- reits angegeben proportional zu n3. Dagegen benötigt das Gesamtschritt- verfahren für einen Schritt eine Zeit, die proportional zu n2 ist. Al- lerdings kommt man mit einem Schritt in den seltesten Fällen aus, so daß der Zeitgewinn nicht alzu groß in's Gewicht fällt. Man erreicht jedoch eine erhebliche Steigerung der Genauigkeit. Das Iterationsver- fahren ist bei Konvergenz numerisch stabil und daher dem Elimentations- verfahren nach Gauß vorzuziehen. Dies gilt insbesondere bei "großen" Systemen (n > 50), bei denen das Elimentationsverfahren zu total un- brauchbaren Ergebnissen führen kann. Besonders vorteilhaft wirkt sich das Iterationsverfahren dann aus, wenn die Matrix A viele Nullen ent- hällt, wie das Beispielsweise bei der Berechnung kubischer Splines oder auch bei der Lösung mancher Differenzialgleichungen der Fall ist. Eine "blinde" Anwendung des Gaußalgorithmus würde diese Nullen "zerstören", was sich dann nachteilig auf die Genauigkeit auswirkt. Ich möchte aber davor warnen dieses Verfahren "blind" einzusetzen. Ist die Konvergenz nicht gewährleistet, so kann es passieren, daß die be- rechneten xi mit der exakten Lösung nichts zu tun haben. In der Regel wird man dieses daran merken, daß die xi bei jedem Iterationsschritt immer größer werden. Dies führt letztlich zu einem Overflow - Error und damit zu einem Abbruch des Programms. Bei linear abhängigen Systemen gibt es 2 Möglichkeiten : entweder das Verfahren konvergiert gegen eine Lösung (obschon es deren unendlich viele gibt) oder es divergiert. Man sollte sich daher die Werte stets ausgeben lassen um nicht die Kontrol- le über die Berechnung zu verlieren. Das oben gesagte gilt auch für Gleichungen oder Gleichungssysteme, die nicht linear sind und mit der oben geschilderter Methode berechnet werden sollen. Und nun viel Spaß beim Austesten und Ausprobieren und ggf. selber schreiben. Heinz Hagemeyer