F
ür mathematische Berechnungen sind Computer geradezu prädestiniert, man muß ihnen nur sagen, wie.Mit der Simpson-Formel lassen sich bestimmte Inte
grale der Form
|b b F(x) = J f(x)dx
|a a
näherungsweise berechnen.
Die Voraussetzungen hierbei sind:
- a,b sind reelle Zahlen und es gilt a < b.
- y = f(x ) sei eine im Intervall [a,b] stetige und differenzierbare Funktion.
Die Lösungsformel lautet:
n
F(a,b,2n) = (h/3) * [f(a) + f(b) + 4 E f(a + (2i-1)h) 1=1
n-1
+ 2 E f(a + 2ih>]
i= 1
wobei das Intervall [a,b] in 2n Teilintervalle der Breite h = (b - a)/2n
geteilt wird.
Das folgende Programm berechnet mit Hilfe der Funktion SIMPSON zu
x
der gegebenen Funktion f(x) = exp(-xz) J exp (t2)dt mit 0
dem Intervall [0,2]
und der Stellengenauigkeit m = 5
die Näherungswerte F(a,b,2i) mit 1 < = i < = k.
Die Berechnung wird abgebrochen, wenn folgendes Krite
rium erfüllt ist:
[F(a,b,2k) - F(a,b,2(k-1>)] < = IF(a,b,2k)l * 0 ,5 * 1 0 - m.
Programmbeschreibung:
In der Darstellung (Bild 1) sei die Unterteilung des Intervalls, sowie die Gewichtung der Funktionswerte an den Stützpunk
ten in der Simpson-Formel verdeutlicht:
Das Integral ergibt sich daraus folgendermaßen:
I = (h/3)*(Fktw ert(a) + Fktw ert(b)+4*(Fktw ert(S 1) + Fktwert(S3) + 2 * Fktwert(S2)
Da bei dieser Intervallunterteilung noch keine großen Ansprüche an die Rechengenauigkeit gestellt werden kön
nen, nimmt man eine feinere Intervalleinteilung vor.
Hierbei benutzen wir folgenden Trick:
Die Anzahl der Teilintervalle wird nicht um zwei erhöht, son
dern verdoppelt (Teilintervallhalbierung).
Der Schritt der Teilintervallhalbierung und seiner Auswir
kung ist in Bild 2 zu sehen.
Dies hat folgende Vorteile:
- Aus der Tatsache, daß die bisherigen Stützpunkte wieder benutzt werden, folgt eine wesentliche Verminderung des Rechenaufwandes.
- Da die Anzahl der Teilintervalle exponentiell ansteigt, kon
vergiert dieses Verfahren sehr schnell.
Dabei fällt folgendes auf:
- Alle bisherigen Stützpunkte, sowohl die mit der Gewich
tung 2 als auch die mit der Gewichtung 4, erhalten bei die
sem Schritt die neue Gewichtung 2.
- Die Anzahl der neuen Punkte (also die mit der Gewichtung 4) ist auf das Doppelte gewachsen.
- Die neuen Stützpunkte sind jeweils um die doppelte Inter
vallbreite voneinander entfernt.
Um diese Erscheinung nochmals zu verdeutlichen, ist in der Darstellung in Bild 3 eine weitere Intervallhalbierung vor
genommen worden.
Durch die stetige Intervallhalbierung wird irgendwann eine genügende Genauigkeit erreicht.
Zur Feststellung dieser Genauigkeit beziehungsweise als Abbruchkriterium für die Unterteilungsschleife eignet sich eine Formel, worin 2 k beziehungsweise 2 k_1 zwei aufeinan
derfolgende Intervallhalbierungen bedeuten. Dieses Abbruchkriterium verwendet schließlich auch M zur Angabe der Zahl der wesentlichen Stellen.
Die bisherigen Erläuterungen des Verfahrens lassen sich einfach in einen Algorithmus umsetzen.
Beim Funktionsaufruf werden eine reelle Funktion f, die untere und obere Integrationsgrenze (— > ug und og, reell), sowie die Rechengenauigkeit m (integer) angegeben. Das Ergebnis der Integration ist natürlich wieder reell.
Damit sieht der Funktionskopf wie folgt aus:
FUNCTION SIMPSON (FUNCTION F(X : REAL) : REAL;
OG, UG : REAL;
M : INTEGER)
: REAL;
In der Variablendeklaration sind folgende Bezeichner fest
gelegt:
SUMUO : REAL
- Summe der Fktwerte der Integrationsgrenzen
G e w i c h t u n g 1 4 . 2 4 1
M
I---- h — — I
(b)
Teilintervallbreite
Bild 1. Gewichtung In Teilintervallbreite: zwei wichtige Fakto
ren für die Simpson-Formel
bish e r i g e G e w i c h t u n g 1 4 2 4 1
neue Gewi c h t u n g *1 4 2 4 2 4 2 4 1
N u m m e r des n euen S t ü tzpunktes
(a)
1 2 3 4
(b)
Bild 2. Verkleinern wir die Teilintervalle, so wird das Ergebnis genauer
bisherige Gewichtung 1 4 2 4 2 4 2 4 1
neue Gewichtung 1 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 1
+— +-—t - f
(a) (b)
Nummer des neuen Stützpunktes
1 2 3 * 5 6 7 8
Bild 3. Jede weitere Halbierung macht das Ergebnis brauch
barer
SUMALT : REAL
- Summe der Fktwerte mit Gewichtung 2 in der Simpson- Formel
SUMNEU : REAL
- neuberechnete Summe der Fktwerte mit Gewichtung 4
INTALT, INTNEU : REAL
- letzter und aktueller Integralwert
DELTA : REAL
- aktuelle Teilintervallbreite (h)
ZNEUPKT : INTEGER
- aktuelle Anzahl der neuzuberechnenden Fktwerte
COUNT : INTEGER
- Schleifenzähler zur Summierung der Fktwerte
ZEHNHOCHM : INTEGER
- »Konstante«, enthält 10 für Abdruckkriterium
Nun zum eigentlichen Algorithmus. Beim Eintritt in das Pro
gramm sind mehrere Größen zu initialisieren:
SUMUO : = F(UG) + F(OG)
- SUMUO wird für den ganzen Fktaufruf so festgelegt.
DELTA : = (OG - UG)/2
- Startwert der Teilintervallbreite
ZEHNHOCHM : = TRUNC(EXP(M*LN(10)) ) = 1 0 IM SUMALT : = 0
INTNEU : = 0 SUMNEU : = 0
Startwerte für ersten Durchlauf ZNEUPKT : = 1
- beim ersten Durchlauf nur 1 neuer Fktwert
Nach der Initialisierung folgt nun die »Integrationsschleife«.
Ihre grundsätzliche Struktur stellt ein REPEAT-UNTIL- Konstrukt dar, wobei als Abbruchkriterium die eingangs erwähnte Formel dient.
In der Schleife werden folgende Aufgaben gelöst:
1 . ) Aufaddieren der neuen Funktionswerte (Abstand: dop
pelte Teilintervallbreite mit Beginn beim 1. Stützpunkt nach der Untergrenze).
2 . ) Wenn das Abbruchkriterium nicht erfüllt ist:
a. ) Verdoppeln der Anzahl der neu zu berechnenden Stütz
stellen = ZNEUPKT
b. ) Halbieren der Intervallbreite
c. ) Sichern des berechneten Integrals für den Vergleich beim nächsten Schleifendurchlauf.
d. ) Umbenennen der Funktionswerte mit Gewichtung 4 zu den Funktionswerten mit Gewichtung 2
e. ) Neuinitialisierung der Summe mit Gewichtung 4 auf 0 In der Funktion SIMPSON sieht die Schleife wie folgt aus:
REPEAT
INTALT := INTNEU;
SUMALT := SUMALT + SUMNEU;
SUMNEU := 0;
FOR COUNT := 1 TO ZNEUPKT DO
SUMNEU := SUMNEU + F(UG+(2*C0UNT-l)*DELTA);
ZNEUPKT := ZNEUPKT*2;
DELTA := DELTA/2;
UNTIL ( Abbruchkriterium e r f ü l l t ) ;
Den letzten Teil der Funktion bildet schließlich die Über
gabe des beim letzten Durchlauf berechneten Integrals, bei dem das Abbruchkriterium erfüllt wurde, an den Funktionsbe
zeichner.
Sollte der zur Verfügung stehende Compiler die Übergabe von Funktionen an eine Funktion nicht erlauben (zum Beispiel Turbo-Pascal), so muß man die Funktion SIMPSON im Haupt
programm realisieren oder man behilft sich, indem man die
_____________ ____________
---OMPJTE
75
Funktion als global vereinbart. Sie muß dann aber vor jedem Programmstart neu eingegeben werden - der Bedienungs
komfort sinkt.
Um die Funktion SIMPSON überhaupt nutzen zu können, benötigt man noch ein Rahmenprogramm, das die zusam
mengesetzte Funktion f(x) im Intervall [a,b] an n + 1 Stützstel
len auf die wesentliche Stellengenauigkeit m tabelliert.
Die Platzhalter der Daten für Intervallanfang, -ende, Stellen
genauigkeit und Stützpunkte wurden als Konstante darge
stellt. Sie könnten jedoch auch ohne weiteres vom Programm eingelesen werden. Somit sieht der Programmkopf wie folgt aus:
PROGRAM RAHMEN ( INPUT;OUTPUT);
CONST A = 0 . 0 ; (»T ab ellieru n gsan fan g ») B = 2 . 0 ; (»T abellierungsende * )
N = 4 0 ; (»N+1 S tü tz s te lle n zur T ab ellieru n g *)
M = 5 ; (» Anzahl der sig n ifik a n te n S te lle n zur In te g ra tio n »)
VAR X : REAL; (» S tü t z s t e l l e a l s Argument fü r die Funktion »)
H : REAL; (» S c h rittw e ite zur T ab ellieru n g ») I : 0 . . N; (» S ch le ife n z ä h le r fü r die
S tü tz s te lle n *)
Darauf folgt der Funktionsdeklarationsteil. Zum einen wird eine Funktion benötigt, die die angegebene, zusammenge
setzte Funktion beinhaltet (FUNCTION FKT). In dieser Funk
tion ist auch die Integrandenfunktion TEILFKT eingebettet, die dann als aktuelle Parameterfunktion für die Integrations
funktion SIMPSON verwendet wird. Hier wird mit dem linken Term EXP(-X*X) das Integral von TEILFKT EXP(X*X) in den Grenzen von UG bis OG auf eine Genauigkeit von M wesentli
chen Stellen berechnet.
Im abschließenden Hauptprogramm findet nur noch die Berechnung der Funktionswerte an den Stützstellen, sowie die tabellarische Ausgabe der Ergebnisse statt.
(Max Moser/Harry Paintner/hg)
Structogramm der Funktion Simpson
Initialisieren der Summe des Anfangs- und Endwerts der Intervallbreite des Integrals,
der Hilfsvariablen für die neuen bzw. alten Integralwerte und der Hilfsvariablen für die Summe der neuen bzw. alten Stützpunkte des Integrals
und des Anfangswerts der Anzahl der Teilintervalle des Integrals.
Teilintervallbreite des Integrals : * (Obergrenze-Untergrenze) / 2 Hilfsvariable für Zehnerpotenz := EXP (Stellengenauigkeit * LN(10))
Sicherung des letzten Integralwerts für die Abbruchent
scheidung
Sicherung der Summe der letztberechneten Stützpunkte Integralvariable auf Null setzen
Schleife von 1 bis Anzahl Stützpunkte Aufaddieren der neuen Stützpunkte Berechnung des neuen Integralwerts Verdoppeln der Anzahl der Stützpunkte
Wiederholen des Blocks, bis die gewünschte Stellengenauigkeit er
reicht ist.
Structogramm der Funktion FKT
Konstante für Integraluntergrenze
Funktionsaufruf SIMPSON, die Funktion TEILFKT, Integraluntergrenze, Argument und Stellengenauigkeit werden als Aktualparameter überge
ben
Funktionswert FKT := EXP ( - x * x ) * Funktionswert von SIMPSON
Structogramm des Hauptprogramms
Konstanten für Intervallanfang, Intervallende, Anzahl Teilintervalle und die gewünschte Stellengenauigkeit
Ausgabe der Überschrift für die Funktionsliste
Teilintervallbreite := (Intervallende -Intervallanfang) / Anzahl der Teil
intervalle
Schleife von 0 bis Anzahl der Teilintervalle
A rgum ent: = Intervallanfang + Schleifenindex * Teilinter
vallbreite
Funktionsaufruf von FKT, Argument und Stellengenauigkeit werden als Aktualparameter übergeben
Ausgabe des Arguments und des Funktionswerts in die Aus
gabeliste
p r o g r a m r a h m e n ( i n p u t » o u t p u t ) )
(»••••••••••••••«••••••«•••••••••••••••••»•••••••••••••••••••»••••••a
(••• ***
( * # * R a h m e n p r o g r a m m z u r N a e h e r u n g s w s i s » n I n t e g r a t i o n n a c h S i m p s o n • • •
( • • • * * *
(•••••«»•••••«••««••••••«»•••••••••••••••••••••••••••••••••••••••••a»
c o n s t a - O . O j
b - 2.0)
n ■ 4 0 |
m - 5)
v a r m( h i r e a l i i i O . . n j
<*** ***
( • • • N a z h i r u n g s H i i s e I n t e g r a t i o n u e b e r e i n e F u n k t i o n f ( x > * * *
<*** ***
f u nction S impson (function f < M I r e a l >t r e a l ) u g,og «real;
m t i n t e g e r )ireal|
v a r s u m u O f S u m a l t » B u m n e u i r e a l ) i n t a l t , i n t n e w , d e l t a i r e a l ) z n e u p k t » c o u n t » z m h n h o c h m s i n t e g e r )
b e g i n
s u m u o * “ f ( u g ) ♦ f ( o g ) ) d e l t a i - ( o g - u g ) / 2 )
z e h n h o c h m I ■ t r u n c ( e x p < m * l n ( 1 0 ) ) ) )
■ u m a l t i - O ) s u m n e u i “ O ) i n t n e u 1 * 0 ) % z n e u p k t i * 1 ) r e p e a t
i n t a l t i » i n t n e u ) s u m a l t t ■ s u m a l t + s u m n e u ) s u m n e u »■» O ) *
■ f o r c o u n t i » 1 t o z n e u p k t d o
s u m n e u » • s u m n e u ♦ f ( u g « - ( 2 * c o u n t - l ) * d e l t a ) ) i n t n e u * - d e l t a • ( s u m u o ♦ 4 • s u m n e u ♦ 2 * s u m a l t ) / 3 ) z n e u p k t - i - z n e u p k t * 2 )
d e l t a I « d e l t a / 2 )
u n t i l ( a b s ( i n t n e u - i n t a l t ) » z e h n h o c h m < ■ a b s ( i n t n e u / 2 ) ) ) S i m p s o n i « i n t n e u )
e n d ) ( ♦ o f f u n c t i o n S i m p s o n * )
(*****•*******************************************•*****■***********•)
(••• ***>
( • • • Z u t a b e l l i e r e n d e F u n k t i o n * • * )
(*•• ***)
(********•******************************************•***************>
f u n c t i o n f k t ( x i r e a l ) m i i n t e g e r ) i r e a l ) c o n s t u g = O . O )
f u n c t i o n t e i l f k t ( t i r e a l ) i r e a l ) b e g i n
t e i l f k t « - e x p ( t * t > ) e n d ) ( * o f f u n c t i o n t e i l f k t * ) b e g i n
f k t i - e x p ( - X * x > * s i m p s o n ( t e i I f k t » u g , x , m ) ) e n d ) ( * o f f u n c t i o n f k t * )
»••
•••
H a u p t p r o g r a m m
)) )>
)
b e g i n
w r i t e l n ) w r i t e l n ) w r i t e l n »
w r i t e l n ( ’ * * * x * * * f ( x > * * * ) )
w r i t e l n ) h I " ( b - a ) / n ) f o r i i = O t o n d o
b e g i n
x » - a i * h j
w r i t e l n ( * ’,xii2»* '»fkt( x , m ) i16))
e n d ) w r i t e l n ) e n d .
Listing. Integration nach Simpson