aki heaeia

Einführung in die Numerik

Olaf Ippish

Universität Heidelberg

Interdisziplinäres Zentrum für Wissenshaftlihes Rehnen

Im Neuenheimer Feld 368, D-69120 Heidelberg

email: olaf.ippishiwr.uni-heidelberg.de

12. Juli 2011

Dieses Dokument enthält Zusatzmaterial inder Form von Beispielen und

AlgorithmensowiederenpraktisherImplementierunginC++zurVorlesung

EinführungindieNumerik gehaltenimSommersemester2011.Esbasiertim

wesentlihen aufArbeitenvonPeter Bastian(IWR,UniversitätHeidelberg).

Inhaltsverzeihnis

1 Motivation 2

1.1 Modellbildungund Simulation . . . 2

1.2 Ein einfahes Beispiel: Das Fadenpendel . . . 4

1.3 Inhaltsübersiht der Vorlesung . . . 12

2 Gleitpunktzahlen 13

2.1 Beispiel zur Fehlerakkumulation . . . 14

3 Gauÿelimination 16

4 Interpolation und Approximation 21

4.1 Motivation . . . 21

4.2 Polynominterpolation . . . 21

4.3 Spline Interpolation . . . 26

4.4 Praktishes zur Diskreten Fourier Analyse . . . 31

4.5 Gauÿ-Approximation . . . 31

5 Quadratur 32 5.1 Newton-Cotes Formeln . . . 32

5.2 Fehlerkontrolle . . . 35

5.3 Gauÿ- und adaptiveQuadratur . . . 36

5.4 Mehrdimensionale Quadratur . . . 39

6 Ein kleiner Programmierkurs 42 6.1 Hallo Welt . . . 42

6.2 Variablenund Typen . . . 44

6.3 Entsheidung . . . 46

6.4 Wiederholung . . . 47

6.5 Funktionen . . . 52

6.6 Funktionsshablonen . . . 54

6.7 HDNUM . . . 56

Lehrbüher Numerik 63

Lehrbüher C++ 63

Weiterführende Literatur 63

1 Motivation

1.1 Modellbildung und Simulation

Die Wissenshaftlihe Methode

•

Experiment: Beobahtewas passiert.

•

^{Theorie: Versuhe die Beobahtung mitHilfe vonModellen zu erklären.}

•

^{Theorieund Experimentwerden sukzessive verfeinert und verglihen, biseine ak-}

zeptable Übereinstimmung vorliegt.

•

^InNaturwissenshaftundTehnikliegenModelleoftinFormmathematisherGlei- hungen vor.

•

^{Oftkönnendie}Modellgleihungennihtgeshlossen(mit Papierund Bleistiftoder Mathematia...) gelöst werden.

Simulation

•

^{Simulation: Gleihungennumerishlösen.}

UndurhführbareExperimentewerden möglih(z. B. Galaxienkollisionen).

Teuere Experimentewerden eingespart (z. B. Modelle imWindkanal).

Parameterstudienshneller durhführbar.

(Automatishe) Optimierungvon Prozessen.

•

Vielfältiger Einsatz in Naturwissenshaft, Tehnik und Industrie: Strömungsbe- rehnung (Klimasimulation),Festigkeit von Bauwerken ...

•

^{Grundlagefür allediese} Anwendungen sind numerishe Algorithmen!

mathematisches Modell konzeptionelles Modell

numerisches Modell Computerprogramm

Realität

wesentliche Prozesse

Wellenausbreitung Transport von Materie Reaktion

Phasenübergänge ...

algebraische Gleichungen Differentialgleichungen Wahrscheinlichkeiten

Funktionen, ...

Objekte: reelle Zahlen,

Näherungsverfahren zur Lösung oben genannter Gleichungen Komplexe SW

SW−Engineering, Qualität High Performance Comp.

Visualisierung

?

Simulation

Effizienz ("Teraflop")

Abbildung1: Prinzipielles Vorgehen imWissenshaftlihen Rehnen.

Die prinzipielle Herangehensweise im Wissenshaftlihen Rehnen zeigt Abbildung 1.

Die erfolgreihe Durhführung einer Simulation erfordert die interdisziplinäre Zusam-

menarbeitvonPhysikern oder Ingenieuren mitMathematikern und Informatikern.

Fehlerquellen

Untershiede zwishen Experiment und Simulation haben vershiedene Gründe:

•

Modellfehler:Einrelevanter Prozess wurdeniht oder ungenau modelliert(Temp.

konstant, Luftwiderstandvernahlässigt, ...)

•

Datenfehler:Messungen von Anfangsbedingungen,Randbedingungen,Wertenfür Parameter sind fehlerbehaftet.

•

Abshneidefehler: Abbruh von Reihen oder Iterationsverfahren, Approximation von Funktionen

•

Rundungsfehler: Reelle Zahlen werden imRehner genähert dargestellt.

Untersuhung vonRundungsfehlernund AbshneidefehleristeinzentralerAspektder

Vorlesung!

1.2 Ein einfahes Beispiel: Das Fadenpendel

Pisa, 1582

Der Student Galileo Galileisitzt in der Kirhe und ihmist langweilig.Er beobahtet

den langsam über ihm pendelnden Kerzenleuhter und denkt: Wie kann ih nur die

Bewegung dieses Leuhters beshreiben?.

Konzeptionelles Modell

Welhe Eigenshaften (physikalishen Prozesse) sind für die gestellteFrage relevant?

•

^{Leuhter istein} Massenpunkt mitder Masse

m

^.

•

^{DerFaden der Länge}

ℓ

^{wird alsrigide und masselos} angenommen.

•

^DerLuftwiderstand wird vernahlässigt.

Nunentwikle mathematishes Modell.

Abbildung 2 zeigt das Fadenpendel, welhes aus dem sogenannten konzeptionellen

Modellresultiert.

Kräfte

•

^{Pendel läuftauf Kreisbahn: Nur}Tangentialkraft istrelevant.

•

^TangentialkraftbeiAuslenkung

φ

^:

F ~ T (φ) = − m g sin(φ)

cos(φ) sin(φ)

.

PSfrag replaements

(0, 0)

ℓ φ

F ~ F ~ N

F ~ T

m

Abbildung 2:Das Fadenpendel.

•

^Alsoetwa:

F ~ T (0) = − mg 0

0

, F ~ T (π/2) = − mg

0 1

.

•

^{Vorzeihen kodiert Rihtung.}

Dies überlegt man sihso. Die Gewihtskraft zeigt immer nahunten, also

F ~ (φ) = mg 0

− 1

.

DieNormalkomponentezeigtimmer inRihtung

~n(φ) = (sin φ , − cos φ) ^T

^{und damit ist}

dieKraftin Normalenrihtung

F ~ _N (φ) = ( F ~ (φ) · ~n(φ))~n =

mg 0

− 1

·

sin φ

− cos φ

sin φ

− cos φ

= mg cos φ

sin φ

− cos φ

.

Damitrehnet man dieTangentialkraftaus

F ~ _T (φ) + F ~ _N (φ) = F ~ (φ)

^aus:

F ~ T (φ) = F ~ (φ) − F ~ N (φ) = mg 0

− 1

− mg cos φ

sin φ

− cos φ

= − mg

cos φ sin φ 1 − cos ² φ

= − mg sin φ

cos φ sin φ

.

Weg, Geshwindigkeit, Beshleunigung

•

^Weg

s(t)

^,Geshwindigkeit

v (t)

^,Beshleunigung

a(t)

^erfüllen:

v(t) = ds(t)

dt , a(t) = dv(t) dt .

•

^Fürden zurükgelegtenWeg (mit Vorzeihen!)gilt

s(t) = ℓφ(t)

^.

•

^{Alsofür die}Geshwindigkeit

v (t) = d s(φ(t))

dt = d ℓφ(t)

dt = ℓ dφ(t) dt

•

^{und die}Beshleunigung

a(t) = d v(φ(t))

dt = ℓ d ² φ(t) dt ² .

Bewegungsgleihung

•

^{Einsetzen indas 2.Newton'she Gesetz}

m a(t) = F (t)

^{liefert nun:}

mℓ d ² φ(t)

dt ² = − mg sin(φ(t)) ∀ t > t 0 .

•

^{DieKraftisthierskalar}(vorzeihenbehafteterBetrag der Tangentialkraft),dawir nurden zurükgelegtenWegbetrahten.

•

^{Ergibtgewöhnlihe} Dierentialgleihung 2.Ordnung für dieAuslenkung

φ(t)

^:

d ² φ(t)

dt ² = − g

ℓ sin(φ(t)) ∀ t > t ₀ .

⁽¹⁾

•

^{Eindeutige Lösungerfordert zwei} Anfangsbedingungen (

t 0 = 0

^):

φ(0) = φ 0 , dφ

(0) = u 0 .

Lösung bei kleiner Auslenkung

•

^{AllgemeineGleihung fürdas Pendel ist shwer analytish zu lösen.}

•

^{Fürkleine Winkel}

φ

^gilt

sin(φ) ≈ φ,

z.B.

sin(0, 1) = 0, 099833417

^.

•

^{DieseNäherung reduziert die Gleihung auf}

d ² φ(t)

dt ² = − g ℓ φ(t).

•

^Ansatz

φ(t) = A cos(ωt)

^liefertmit

φ(0) = φ 0

^,

dφ

dt (0) = 0

^{dann dieaus der Shule}

bekannteFormel

φ(t) = φ ₀ cos r g

ℓ t

(3)

Volles Modell; Verfahren 1

•

^{Löse das volleModell mitzweinumerishen Verfahren.}

•

^{ErsetzeGleihung zweiter Ordnung durh zwei Gleihungenerster Ordnung:}

dφ(t)

dt = u(t), d ² φ(t)

dt ² = du(t) dt = − g

ℓ sin(φ(t)).

•

^ErsetzeAbleitungen durh Dierenzenquotienten:

φ(t + ∆t) − φ(t)

∆t ≈ dφ(t)

dt = u(t), u(t + ∆t) − u(t)

∆t ≈ du(t)

dt = − g

ℓ sin(φ(t)).

•

^Mit

φ ⁿ = φ(n∆t)

^,

u ⁿ = u(n∆t)

^{erhält man Rekursion(Euler):}

φ ⁿ⁺¹ = φ ⁿ + ∆t u ⁿ φ ⁰ = φ 0

⁽⁴⁾

u ⁿ⁺¹ = u ⁿ − ∆t (g/ℓ) sin(φ ⁿ ) u ⁰ = u ₀

⁽⁵⁾

-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2

0 0.5 1 1.5 2 2.5 3 3.5 4

Auslenkung

Zeit

Konvergenz Differenzenquotient (Euler), phi=0.1 vereinfachtes Modell

dt=0.2 dt=0.1 dt=0.01 dt=0.001

Abbildung3:Simulation des Fadenpendels (volles Modell) bei

φ 0 = 0.1 ≈ 5.7 ^◦

^mitdem

Eulerverfahren.

Volles Modell; Verfahren 2

•

^Nutze Näherungsformel für diezweite Ableitung, sog. Zentraler Dierenzenquoti- ent:

φ(t + ∆t) − 2φ(t) + φ(t − ∆t)

∆t ² ≈ d ² φ(t) dt ² = − g

ℓ sin(φ(t)).

•

^Auösennah

φ(t + ∆t)

^ergibtRekursionsformel (

n ≥ 2

^):

φ ⁿ⁺¹ = 2φ ⁿ − φ ^{n − 1} − ∆t ² (g/ℓ) sin(φ ⁿ )

⁽⁶⁾

mitder Anfangsbedingung

φ ⁰ = φ 0 , φ ¹ = φ 0 + ∆t u 0 .

⁽⁷⁾

(Diezweite Bedingung kommt aus dem Eulerverfahren oben).

Abbildung 3 zeigt das Eulerverfahren inAktion. Für festen Zeitpunkt

t

^und

∆t → 0

konvergiert dasVerfahren.Fürfestes

∆t

^und

t → ∞

^{nimmtdasVerfahrenimmergröÿere}

Werte an.

Abbildung4zeigtzum Vergleihdas zentrale Verfahrenfürdiegleihe Anfangsbedin-

gung.ImUntershied zum explizitenEulersheint derFehlerbeifestem

∆t

^und

t → ∞

-0.15 -0.1 -0.05 0 0.05 0.1 0.15

0 0.5 1 1.5 2 2.5 3 3.5 4

Auslenkung

Zeit

Konvergenz zentriertes Verfahren, phi=0.1 vereinfachtes Modell

dt=0.2 dt=0.1 dt=0.01 dt=0.001

Abbildung4:Simulation des Fadenpendels (volles Modell) bei

φ 0 = 0.1 ≈ 5.7 ^◦

^mitdem

zentralen Verfahren.

Nunkönnenwir dasvolleModellmitdemvereinfahten Modellvergleihenundsehen

welhe Auswirkungen die Annahme

sin φ ≈ φ

^{auf das Ergebnis hat. Abbildung 5 zeigt}

dienumerishe Simulation.Selbstbei

28.6 ^◦

^istdieÜbereinstimmungnoheinigermaÿen passabel.Für gröÿere Auslenkungen istdas vereinfahteModell völligunbrauhbar, die

Form der Shwingung ist keinKosinusmehr.

Eine Geothermieanlage

Wir betrahten die Modellierung und Simulation einer Geothermieanlage. Den she-

matishen Aufbau zeigt die Abbildung 6. Kaltes Wasser ieÿt in einer Bohrung nah

unten, wird erwärmt und ineinem isoliertenInnenrohr wiedernah oben geführt.

•

Grundwasserströmung gekoppelt mitWärmetransport.

•

^{Welhe Leistung erzielt soeine Anlage?}

Modell für eine Geothermieanlage

•

^{Strömung des Wassers inund um das Bohrloh}

∇ · u = f, u = − K

µ ( ∇ p − ̺ w G)

-0.6 -0.4 -0.2 0 0.2 0.4 0.6

0 0.5 1 1.5 2 2.5 3 3.5 4

Auslenkung

Zeit

Zentriertes Verfahren, phi=0.5 vereinfachtes Modell

dt=0.01 dt=0.0001

-4 -3 -2 -1 0 1 2 3 4

0 0.5 1 1.5 2 2.5 3 3.5 4

Auslenkung

Zeit

Zentriertes Verfahren, phi=3.0

vereinfachtes Modell dt=0.0001

-4 -3 -2 -1 0 1 2 3 4

0 2 4 6 8 10

Auslenkung

Zeit

Zentriertes Verfahren, phi=3.14

vereinfachtes Modell dt=0.0001

Abbildung5:Vergleih von vollem und vereinfahtem Modell (jeweils in rot) bei den

Winkeln

φ = 0.5, 3.0, 3.14

^{gerehnet mit dem zentralen Verfahren.}

PSfrag replaements

Vorlauf Rüklauf

Wasserstrom

Wärmestrom

Abbildung6:Shematisher Aufbaueiner Geothermieanlage.

•

^{Transport der Wärme durh Konvektion und} Wärmeleitung

∂(c e ̺ e T )

∂t + ∇ · q + c _w ̺ _w f T = g,

q = c w ̺ w uT − λ ∇ T

inAbhängigkeitdiverserParameter:Bodendurhlässigkeit,Wärmekapazität,Dih-

te, Wärmeleitfähigkeit, Pumprate sowie Rand- und Anfangsbedingungen.

Abbildung 6 zeigt die Entzugsleistung so einer Anlage für die ersten hundert Tage

Betriebszeit. Dabei wurden plausibleWerte für dieModellparametergewählt.

Shadstoausbreitung

Als weiteres Beispiel nennen wir die Modellierung und Simulation einer Shadsto-

ausbreitungimBoden,siehe Abbildung8.DerShadstowird hieralsnihtmitWasser

mishbar und shwerer als Wasser angenommen (Bestimmte Reinigungsmittel haben

diese Eigenshaften).

•

^{Wo erreiht der Shadsto welhe} Konzentrationen?

•

^{Wiebekommt man den Shadsto wieder weg?}

•

^{Wohin bewegt sihgelöster Shadsto?}

300000 400000 500000 600000 700000 800000 900000 1e+06 1.1e+06 1.2e+06 1.3e+06

0 10 20 30 40 50 60 70 80 90 100

Entzugsleistung [Watt]

Zeit [Tage]

Entzugsleistung im kalibrierten System

Abbildung7:Entzugsleistung der Geothermieanlage über die Zeit für bestimmte

Systemparameter.

1.3 Inhaltsübersiht der Vorlesung

Wirwerden in dieser Vorlesung diefolgendenThemengebiete behandeln:

•

Grundbegrie, Gleitpunktzahlen, Gleitpunktarithmetik

•

^{DirekteMethoden zur Lösung linearer} Gleihungssysteme

•

Interpolationund Approximation

•

^Numerishe Integration

•

Iterationsverfahren zur Lösung linearerGleihungssysteme

•

Iterationsverfahren zur Lösung nihtlinearerGleihungssysteme

•

^{Eigenwerte und} Eigenvektoren

Ausblikauf weiterführende Vorlesungen

•

^GewöhnliheDierentialgleihungenbeshreibenvielezeitabhängigeProzesse(Nu- merik1).

Abbildung 8:Ausbreitung eines niht Wasser niht mishbaren Shadstoes im Boden.

•

^Praktishe Anwendungenführen oftauhaufpartielleDierentialgleihungen.Ge- suht sind dannFunktionen inmehreren Raumdimensionen(Numerik 2).

•

^Tehnishe Anwendungen erfordern Optimierung (Numerik 3).

•

^Eziente Programmierung und Visualisierung der Ergebnisse sind sehr wihtig (z.B.Objektorientiertes ProgrammierenimWissenshaftlihen Rehnen).

2 Gleitpunktzahlen

Zahlen im Computer

Alle Programmiersprahen stellen elementare Datentypen zur Repräsentation von

Zahlenzur Verfügung. In C/C++:

unsigned int , unsigned short , unsigned long

N ₀

int , short , long

Z

float, double

R

omplex< float >, omplex<double>

C

Diese sind Idealisierungen der Zahlenmengen

N ₀ , Z , R , C

.

Beiunsigned int ,int...bestehtdieIdealisierungdarin,dasseseinegröÿte(bzw.kleinste)

darstellbareZahl gibt. Ansonsten sind dieErgebnisse exakt.

Bei float und double kommt hinzu, dass die meisten innerhalb des erlaubten Bereihs

liegenden Zahlennur näherungsweise dargestellt werden können.

2.1 Beispiel zur Fehlerakkumulation

Potenzreihe für

e ^x

• e ^x

^{lässt sihmiteiner} Potenzreihe berehnen:

e ^x = 1 + X ∞

n=1

x ⁿ

n! = 1 + X ∞

n=1

y n .

•

^Diesformulierenwir rekursiv:

y 1 = x, S 1 = 1 + y 1 ,

(Anfangswerte).

Rekursion:

y _n = x

n y _{n − 1} , S _n = S _{n − 1} + y _n .

•

^Probierevershiedene Genauigkeiten und Werte von

x

^.

Positives Argument

•

^Für

x = 1

^und float-Genauigkeiterhalten wir:

1.000000000000000e+00 1 2.000000000000000e+00

5.000000000000000e-01 2 2.500000000000000e+00

1.666666716337204e-01 3 2.666666746139526e+00

4.166666790843010e-02 4 2.708333492279053e+00

8.333333767950535e-03 5 2.716666936874390e+00

1.388888922519982e-03 6 2.718055725097656e+00

1.984127011382952e-04 7 2.718254089355469e+00

2.480158764228690e-05 8 2.718278884887695e+00

2.755731884462875e-06 9 2.718281745910645e+00

2.755731998149713e-07 10 2.718281984329224e+00

0.000000000000000e+00 100 2.718281984329224e+00

ex 2.718281828459045e+00

...also7 gültigeZiern.

•

^Für

x = 5

^...

9.333108209830243e-06 21 1.484131774902344e+02

ex 1.484131591025766e+02

...dito.

Negatives Argument

•

^Für

x = − 1

^und float -Genauigkeiterhalten wir:

2.755731998149713e-07 10 3.678794205188751e-01

-2.505210972003624e-08 11 3.678793907165527e-01

2.087675810003020e-09 12 3.678793907165527e-01

ex 3.678794411714423e-01

...6 gültigeZiern.

•

^Für

x = − 5

-5.000000000000000e+00 1 -4.000000000000000e+00

1.250000000000000e+01 2 8.500000000000000e+00

-2.083333396911621e+01 3 -1.233333396911621e+01

2.604166793823242e+01 4 1.370833396911621e+01

-2.333729527890682e-02 15 1.118892803788185e-03

7.292904891073704e-03 16 8.411797694861889e-03

1.221854423194557e-10 28 6.737461313605309e-03

0.000000000000000e+00 100 6.737461313605309e-03

ex 6.737946999085467e-03

nurnoh 4 gültigeZiern!

Noh kleineres Argument

•

^Für

x = − 20

^und float-Genauigkeitsind ... -2.000000000000000e+01 1 -1.900000000000000e+01

2.000000000000000e+02 2 1.810000000000000e+02

-1.333333374023438e+03 3 -1.152333374023438e+03

6.666666992187500e+03 4 5.514333496093750e+03

-2.666666796875000e+04 5 -2.115233398437500e+04

-2.611609750000000e+06 31 -1.011914250000000e+06

1.632256125000000e+06 32 6.203418750000000e+05

-9.892461250000000e+05 33 -3.689042500000000e+05

5.819095000000000e+05 34 2.130052500000000e+05

-3.325197187500000e+05 35 -1.195144687500000e+05

1.847331718750000e+05 36 6.521870312500000e+04

-4.473213550681976e-07 65 7.566840052604675e-01

1.355519287926654e-07 66 7.566841244697571e-01

-4.046326296247571e-08 67 7.566840648651123e-01

1.190095932912527e-08 68 7.566840648651123e-01

ex 2.061153622438557e-09

keine Ziern mehr gültig. Das Ergebnisist um 8 Gröÿenordnungen daneben!

Erhöhen der Genauigkeit

•

^Für

x = − 20

^und double-Genauigkeiterhält man -1.232613988175268e+07 27 -5.180694836889297e+06

8.804385629823344e+06 28 3.623690792934047e+06

1.821561256740375e-24 94 6.147561828914626e-09

-3.834865803663947e-25 95 6.147561828914626e-09

ex 2.061153622438557e-09

Immernohum einen Faktor 3daneben!

•

^{Erst mitvierfaher} Genauigkeit erhältman

7.0937168371834023136209731491427e-42 118 2.0611536224385583392700458752947e-09

ex 2.0611536224385578279659403801558e-09

15gültige Ziern(beia 30Ziern Rehengenauigkeit).

Fragen

•

^{Was bedeutet überhaupt} Rehengenauigkeit?

•

^Welhe Genauigkeitkönnen wir erwarten?

•

^{Wo kommen diese Fehler her?}

•

^{Wiewerden solhe} Kommazahlen dargestellt und verarbeitet?

ObigeBerehnungenwurdenmitdenPaketenqdundarpre(beidehttp://rd. lbl.

gov/ ~dhbailey/mpdist/ )durhgeführt.qderlaubtbiszuvierfahedoubleGenauigkeit,

arpre beliebige Genauigkeit. Die GNU multipreision library (http://gmplib.org/ )

isteine Alternative.

3 Gauÿelimination

Algorithmus 4.1 (Gauÿ-Elimination)

Input:

A ∈ R ⁿ × n

,

b ∈ R ⁿ

Output:

x ∈ R ⁿ

for (

k = 1

^;

k < n

^;

k = k + 1

^{) do}

Finde

r ∈

^{{k,...,n} so dass}

a rk 6 = 0

^{;(sonst Fehler);}

if (

r 6 = k

^{) then {taushe Zeile}

k

^mitZeile

r

^}

for (

j = k

^;

j ≤ n

^;

j = j + 1

^)do

t = a kj

^;

a kj = a rj

^;

a rj = t

^;

end for

t = b k

^;

b k = b r

^;

b r = t

^;

endif

{Update der unteren Zeilen}

for (

i = k + 1

^;

i ≤ n

^;

i = i + 1

^)do

q = a ik /a kk

^;

for (

j = k + 1

^;

j ≤ n

^;

j = j + 1

^{) do}

a ij = a ij − q · a kj

^;

end for

b i = b i − q · b k

^;

endfor

{Rükwärtseinsetzen}

for (

k = n

^;

k ≥ 1

^;

k = k − 1

^{) do}

t = 0

^;

for (

j = k + 1

^;

j ≤ n

^;

j = j + 1

^)do

t = t + a kj · x j

^;

endfor

x k = (b k − t)/a kk

end for

Beispiel 4.2

Hier sind keine Zeilenvertaushungen notwendig. Das Pivotelement ist jeweils durh

einenKasten gekennzeihnet.



 



2 4 6 8 40

16 33 50 67 330 4 15 31 44 167 10 29 63 97 350



 

 →



 



2 4 6 8 40

0 1 2 3 10

0 7 19 28 87 0 9 33 57 150



 



→



 



2 4 6 8 40 0 1 2 3 10 0 0 5 7 17 0 0 15 30 60



 

 →



 



2 4 6 8 40 0 1 2 3 10 0 0 5 7 17 0 0 0 9 9



 



Shlieÿlihliefert Rükwärtseinsetzen:

x 4 = 9/9 = 1 , x 3 = (17 − 7 · 1)/5 = 2 ,

x 2 = (10 − 2 · 2 − 3 · 1)/1 = 3 , x 1 = (40 − 4 · 3 − 6 · 2 − 8 · 1)/2 = 4 .

Algorithmus 4.7.1 zur LR-Zerlegung

Input:

A ∈ R ⁿ × n

,

b ∈ R ⁿ

(wird übershrieben) Output:

L ∈ R ^{n × n} in a _ij , j < i, l _ii = 1

^implizit

R ∈ R ^{n × n} in a ij , j ≥ i

p : { 0, . . . , n − 1 } → { 0, . . . , n − 1 }

for (

k = 1

^;

k < n

^;

k = k + 1

^{) do}

Finde

r ∈

^{{k,...,n} so dass}

a rk 6 = 0

^{;(sonst Fehler);}

if (

r 6 = k

^{) then {taushe Zeile}

k

^mitZeile

r

^}

for (

j = 1

^;

j ≤ n

^;

j = j + 1

^{) do}

t = a kj

^;

a kj = a rj

^;

a rj = t

^;

end for

endif

p k = r

^{; {merke}Permutation}

{Update der unteren Zeilen}

for (

i = k + 1

^;

i ≤ n

^;

i = i + 1

^)do

a ik = a ik /a kk

^;

for (

j = k + 1

^;

j ≤ n

^;

j = j + 1

^{) do}

a _ij = a _ij − a _ik · a _kj

^;

end for

endfor

end for

{Permutation von

b

^}

for (

k = 1

^;

k < n

^;

k = k + 1

^{) do}

if (

p k 6 = k

^)then

t = b k

^;

b k = b p k

^;

b p k = t

^;

endif

end for

{Vorwärtseinsetzen}

for (

k = 1

^;

k ≤ n

^;

k = k + 1

^)do

t = 0

^;

for (

j = 1

^;

j < k

^;

j = j + 1

^{) do}

t = t + a kj · x j

^;

endfor

x k = b k − t

^;

end for

{Rükwärtseinsetzen}

for (

k = n

^;

k ≥ 1

^;

k = k − 1

^{) do}

t = 0

^;

for (

j = k + 1

^;

j ≤ n

^;

j = j + 1

^)do

t = t + a kj · x j

^;

endfor

x k = (x k − t)/a kk

^;

end for

Algorithmus 4.7.2 zur LR-Zerlegung (Andere Permutation)

Input:

A ∈ R ^{n × n}

,

b ∈ R ⁿ

(wird übershrieben) Output:

L ∈ R ⁿ × n in a ij , j < i, l ii = 1

^implizit

R ∈ R ^{n × n} in a ij , j ≥ i

p : { 0, . . . , n − 1 } → { 0, . . . , n − 1 }

for (

k = 1

^;

k ≤ n

^;

k = k + 1

^)do

p k = k

^;

end for

for (

k = 1

^;

k < n

^;

k = k + 1

^{) do}

Finde

r ∈

^{{k,...,n} so dass}

a rk 6 = 0

^{;(sonst Fehler);}

if (

r 6 = k

^{) then {taushe Zeile}

k

^mitZeile

r

^}

for (

j = 1

^;

j ≤ n

^;

j = j + 1

^{) do}

t = a kj

^;

a kj = a rj

^;

a rj = t

^;

end for

t = p _k

^;

p _k = p _r

^;

p _r = t

^;

endif

{Update der unteren Zeilen}

for (

i = k + 1

^;

i ≤ n

^;

i = i + 1

^)do

a ik = a ik /a kk

^;

for (

j = k + 1

^;

j ≤ n

^;

j = j + 1

^{) do}

a ij = a ij − a ik · a kj

^;

end for

endfor

end for

{Permutation von

b

^}

for (

k = 1

^;

k ≤ n

^;

k = k + 1

^)do

x _k = b _{p k}

^;

end for

{Vorwärtseinsetzen}

for (

k = 1

^;

k ≤ n

^;

k = k + 1

^)do

t = 0

^;

for (

j = 1

^;

j < k

^;

j = j + 1

^{) do}

t = t + a kj · x j

^;

endfor

x k = x k − t

^;

end for

{Rükwärtseinsetzen}

for (

k = n

^;

k ≥ 1

^;

k = k − 1

^{) do}

t = 0

^;

for (

j = k + 1

^;

j ≤ n

^;

j = j + 1

^)do

t = t + a kj · x j

^;

endfor

x k = (x k − t)/a kk

^;

end for

Algorithmus 4. zu LS with QR-Zerlegung

Input:

x ∈ R ⁿ

,

j ∈ N ⁺

Output:

v ∈ R ⁿ

HouseholderVetor(

x

^,

j

⁾

begin

n = length(x)

^;

v = 0

^;

µ = 0

^;

for (

k = j

^;

k ≤ n

^;

k = k + 1

^{) do}

µ = µ + x k ∗ x k

^;

v k = x k

^;

end for

if (

µ 6 = 0

^{) then}

β = x 1 + sign(x 1 )µ

^;

for (

k = j + 1

^;

k ≤ n

^;

k = k + 1

^{) do}

v k = v k /β

^;

end for

end if

v 1 = 1

return

v

^;

end

Algorithmus 4. zu LS with QR-Zerlegung

Input:

A ∈ R ^{m × n}

,

v ∈ R ^m

Output:

A ∈ R ^m × n

(wird übershrieben)

HouseholderMultiplikation(

A

^,

v

⁾

begin

β = − 2/(v ^T v)

^;

w = βA ^T v

^;

A = A + vw ^T

^;

end

Algorithmus 4. zu LS with QR-Zerlegung

Input:

A ∈ R ^{m × n}

Output:

n

^Vektoren

v ^(k) ∈ R ^{m × n} in a ij , j < i, v ₁ ^(k) = 1

^implizit

R ∈ R ^m × n in a ij , j ≥ i

HouseholderQR(

A

⁾

begin

for

j = 1

^;

j ≤ n

^;

j = j + 1

^do

v k

w = βA ^T v

^;

A = A + vw ^T

^;

end for

end

4 Interpolation und Approximation

4.1 Motivation

Die Abbildungen 9 bis11 zeigen eine Anwendung von Polynomen bei der Kurvenkom-

pressionin der Computergraphik.

DieLage einesstarrenKörpers imRaumwird durh6Zahlenfestgelegt(3fürdiePo-

sitionund3fürdieOrientierung),diesihmitderZeitändernkönnen.Eineäquidistante

Shrittweite erforderteinenhohenSpeiheraufwand,umbeishnellenPositionsänderun-

gen eine gute Genauigkeit erreihen zu können. Bei einer adaptiven Shrittweitenwahl

werden möglihst wenig Zeitpunkte ausgewählt, aber so, dass ein vorgegebener Fehler

niht übershritten wird. Diese Anwendung haben Eri Shneider, Manuel Jerger und

BenjaminJillihimRahmeneines Software-Praktikums imSommersemester 2008erar-

beitet (Vielen Dank für dietollen Bilder!).

4.2 Polynominterpolation

Die Abbildung12zeigt dieMonome bis zum Grad 6.

Abbildung13zeigtdieLagrange-PolynomevomGrad6beiäquidistanten Stützstellen

auf

[0, 1]

^.

Zu interpolieren seidie folgendeWertetabellemit 4Einträgen:

x i y i

0 1.0000 2 0.4546 7 0.0938 10 − 0.0544

Abbildung14zeigtdaszugehörigeInterpolationspolynomsowiedieskaliertenLagrange-

Polynome

y i L ⁽³⁾ _i

^.

Abbildung16illustriertdasWahsender Lagrange-PolynomeweitwegvonderStütz-

stelle

x

^{an der}

L ^(n)(x) _i = 1

^gilt.

Beispiel 5.8 zur numerishen Dierentiation

Wir wollen die zweite Ableitung von

f(x) = sinh(x)

^für

x = 0.6

^{mit dem zweiten}

Dierenzenquotient ermitteln:

d ²

dx ² sinh(x) ≈ sinh(x + h) − 2 sinh(x) + sinh(x − h)

h ²

Abbildung 9:Kurvenkompression in der Computergraphik: Die Szene.

Abbildung10:KurvenkompressioninderComputergraphik:Stützpunktederunkompri-

mierten Kurve.

Abbildung11:Kurvenkompression in der Computergraphik: Stützpunkte der kompri-

mierten Kurve.

0 0.2 0.4 0.6 0.8 1 1.2

0 0.2 0.4 0.6 0.8 1

y

x i=0

i=1 i=2 i=3 i=4 i=5 i=6

Abbildung 12:Die Monome biszum Grad 6.

zur Erinnerung:

sinh(x) = 1

2 (e ^x − e ^{− x} ), d

dx sinh = cosh = 1

2 (e ^x + e ^{− x} ), d ²

dx ² sinh(x) = sinh(x).

Mit double Genauigkeit erhältman den Wert

sinh(0.6) = 6, 366535821482 · 10 ^{− 1} .

Dagegen liefertdie numerishe Dierentiation diefolgende Tabelle

h

Dierenzenquotient

1 · 10 ^{− 1} 6.371 · 10 ^{− 1} 1 · 10 ^{− 2} 6.3665888 · 10 ^{− 1} 1 · 10 ^{− 3} 6.366536352 · 10 ^{− 1} 1 · 10 ^{− 4} 6.3665358540 · 10 ^{− 1}

1 · 10 ^{− 5} 6.3665017 · 10 ^{− 1}

^Auslöshung

1 · 10 ^{− 6} 6.3671 · 10 ^{− 1}

.

.

.

1 · 10 ^{− 10} 1.1102 · 10 ⁴ !

-1.5 -1 -0.5 0 0.5 1 1.5

0 0.2 0.4 0.6 0.8 1

y

x

i=0 i=1 i=2 i=3 i=4 i=5 i=6

Abbildung13: DieLagrange-Polynome

L ⁽⁶⁾ _i (x)

^{vomGrad 6.}

NumerisheDierentationistsehranfälliggegenüberRundungsfehlern. Möglihe Ab-

hilfebietet dieExtrapolation.

4.3 Spline Interpolation

Wirbetrahten dieInterpolation der folgendendrei Funktionen

f 1 (x) = exp( − x ² )

ⁱⁿ

[ − 10, 10],

⁽⁸⁾

f ₂ (x) =

cos ² (x) | x | < ^π / 2

0 | x | ≥ ^π / 2

in

[ − π, π],

⁽⁹⁾

f 3 (x) =

− 1 x < ¹ / ²

+1 x ≥ ¹ / ²

ⁱⁿ

[0, 1],

⁽¹⁰⁾

mittelsPolynomen,

S _h ^1,0

^und

S _h ^3,2

^{(mit natürlihen} Randbedingungen).

0 0.2 0.4 0.6 0.8 1

y

f1(x)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

y

f2(x)

-1 -0.5 0 0.5 1

y

f3(x)

-0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0 2 4 6 8 10

y

x

Interpolationspolynom p

p y0*L0 y1*L1 y2*L2 y3*L3 Daten

Abbildung 14:Interpolationspolynom zur Wertetabelleaus Beispiel .

Die Abbildung17 zeigt dieInterpolation der Funktion

f 1 (x)

^.

Die Abbildung18 zeigt dieInterpolation der Funktion

f 2 (x)

^.

Die Abbildung19 zeigt dieInterpolation der Funktion

f 3 (x)

^.

Wir lernen:

•

InterpolationmitPolynomensteigendenGradesanäquidistantenStützstellenshlägt inallen Fällenfehl, d.h. der Interpolationsfehler steigt mitdem Grad an.

•

^{KubisheSplines}konvergieren undlieferneinenglattenVerlauf.Allerdingskommt es zu mögliherweise unphysikalishen Unter- bzw. Übershwingern. Diese sind

aberimFalle von

f 3 (x)

^{um die}Sprungstelle lokalisiert.

•

^{Stükweise lineareFunktionen haben diesen Defekt niht.}

Wir wollen nun den Interpolationsfehler noh experimentell bestimmen.

Fehler bei Interpolation der Funktion

f 1 (x) = e ^{− x 2}

^:

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

0 1 2 3 4 5 6

y

x

Polynominterpolation von f(x) = sin(2x)

sin(2*x) n=4 n=5 n=12

-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1

-4 -2 0 2 4

y

x

Polynominterpolation von f(x) = 1/(1+x*x)

1.0/(1.0+x*x) n=4 n=8 n=12

Abbildung15:Interpolation der Funktionen

sin(2x)

^{(oben) und}

_1+x ^{1 2}

^{(unten) mit äqui-}

distanten Stützstellen und vershiedenen Polynomgraden.

-3 -2 -1 0 1 2 3 4 5 6 7

0 0.2 0.4 0.6 0.8 1

y

x

Lagrange Polynome mit steigendem Grad i=2, n=4 i=3, n=6 i=4, n=8 i=5, n=10

-200 -150 -100 -50 0 50 100

0 0.2 0.4 0.6 0.8 1

y

x

Lagrange Polynome mit steigendem Grad i=5, n=10 i=6, n=12 i=7, n=14 i=8, n=16

Abbildung16: Die Lagrange-Polynome

L ⁽ⁿ⁾ _n/2

^für

n = 4, 6, 8, 10, 12, 14, 16

^.

n S _h ^1,0 S _h ^3,2 P n

4 6.045 e − 01 7.420 e − 01 8.038 e − 01 6 4.447 e − 01 5.612 e − 01 9.999 e − 01 8 3.002 _e − 01 3.918 _e − 01 2.311 _e + 00 10 1.774 e − 01 2.464 e − 01 5.949 e + 00 16 1.060 e − 01 2.753 e − 02

32 6.946 _e − 02 7.083 _e − 03 64 2.241 e − 02 3.316 e − 04 128 5.974 e − 03 1.918 e − 05 256 1.517 e − 03 1.173 e − 06 512 3.809 e − 04 7.289 e − 08 1024 9.533 e − 05 4.549 e − 09

Angegeben ist der maximale Fehler an einem Punkt. Polynome konvergieren niht.

Stükweise linear konvergiert mit

h ²

^{(d. h.}

e 2n /e n = ( ¹ / ² ) ²

^{), kubishe Splines mit}

h ⁴

(d. h.

e _2n /e _n = ( ¹ / 2 ) ⁴

^).

InbeidenFällengiltdiesnur,wenn

n

^{genügendgroÿ,mansprihtvon}asymptotisher Konvergenz.

Fehler bei Interpolation der Funktion

f 2 (x) =

cos ² (x) x < ^π / ² 0 x ≥ ^π / ²

^:

n S _h ^1,0 S _h ^3,2

4 1.052e − 01 1.649e − 01 8 1.052e − 01 4.498e − 02 16 3.518e − 02 8.434e − 03 32 9.423e − 03 1.945e − 03 64 2.396e − 03 4.764e − 04 128 6.015e − 04 1.184e − 04 256 1.505e − 04 2.958e − 05 512 3.764e − 05 7.394e − 06 1024 9.412e − 06 1.848e − 06

In diesem Fall konvergiert der maximale Fehler auh im Falle kubisher Splines nur

mit

h ²

^.

Dies liegt daran, dass

f ₂ ^′′ (x)

^{unstetigamPunkt}

x = ^π / 2

^{ist(springt von 2 auf0).}

Die dritte Ableitungexistiert niht mehr.

Für dieInterpolationmitstükweisen Polynomen merken wir uns:

Jehöherder(abshnittsweise)PolynomgradumsoshnellerkonvergiertdasVerfahren.

Imallgemeinen erhält man

O(h ^k+1 )

^{Konvergenz für Polynome vom Grad}

k

^.

Dies giltallerdingsnurdann, wenndiezu interpolierende Funktiongenügendof die-

renzierbarist. Istdies nihtder Fallso lohntalsoauhdieVerwendung vonPolynomen

4.4 Praktishes zur Diskreten Fourier Analyse

Abbildung 20 zeigt einige Beispiele für Spektren. Die Konstante im Zeitbereih hat

einen Puls alsSpektrum. Umgedreht hat ein Puls im Ortsbereih einkonstantes Spek-

trum.Shlieÿlihwird nohdas Spektrumeines Dreieks- bzw. Rehteksignals gezeigt.

Abbildung21zeigtdieInterpolationvonDreik-bzw.Rehteksignal beiVorgabevon

jeweilsaht Datenpunkten.

Abbildung 22 illustriert die Verbesserung der Annäherung an die zu interpolierende

Funktionbeisteigendem Parameter

n

^.

Abbildung22illustriertdieVerbesserungder AnnäherungbeiunstetigenFunktionen,

wenn an der Sprungstelle der Mittelwert vorgeshrieben wird. Wir verwenden einmal

n = 15

(Sprungstelleist Interpolationspunkt, Mittelwert wird vorgeshrieben) und

n = 16

(Sprungstelleist kein Interpolationspunkt).

-8000 -6000 -4000 -2000 0 2000 4000 6000

-200 -150 -100 -50 0 50 100 150 200

y

x

4.5 Gauÿ-Approximation

Beispiel 5.23: Fourierreihe

Für

N = 2m + 1, m ∈ N

ist

Ψ F = n

√ 1

2π , 1

√ π cos(x), . . . , 1

√ π cos(mx),

√ 1

π sin(x), . . . , 1

√ π sin(mx)

eine OrthonormalbasisvonFunktionenauf dem Intervall

[ − π, π]

^.

Damit giltdann:

g (x) = a ₀ 2 +

X m

k=1

{ a _k cos(kx) + b _k sin(kx) }

und

a _k = 1 π

Z π

− π

f (x) cos(kx)dx, k = 0, . . . , m

b k = 1 π

Z π

− π

f (x) sin(kx)dx, k = 1, . . . , m

FürunendlihvieleGlieder(

m = ∞

^{)nenntman dieReihe F}ourier-Reihe.Diesekonver- giert gegenein Element aus

L ² ( − π, π)

^.

Beispiel 5.24: Finite-Elemente Diskretisierung

Die numerishe Lösung der partiellenDierentialgleihung

∂ ² u

∂x ² + ∂ ² u

∂y ² = f in Ω

(Poisongleihung) mitHilfe der Methode der Finiten-Elementeführtauf dieAufgabe:

Finde

u ∈ S

^{, so dass}

Z

Ω

∇ u · ∇ ϕdx = Z

Ω

f ϕdx ∀ ϕ ∈ S

5 Quadratur

5.1 Newton-Cotes Formeln

Wirbetrahten folgendebestimmteIntegrale:

(i)

^{Eine einfahe, unendlihoft} dierenzierbare Funktion:

Z ^π / ² 0

sin(x) dx = 1.

(ii)

^{Eine glatte Funktionaber mitgroÿen höheren} Ableitungen:

Z 1

− 1

1

10 ^{− 5} + x ² dx = 9.914588332462438 · 10 ² . (iii)

^{Eine niht unendlih oft} dierenzierbareFunktion (Halbkreis):

Z 1

− 1

√ 1 − x ² dx = ^π / ² .

Summierte Trapezregelfür

(i)

^.

I

^{Fehler #Fktausw.}

9.480594489685199e-0 1 5.1941e-02 3

9.871158009727754e-0 1 1.2884e-02 5

9.967851718861696e-0 1 3.2148e-03 9

9.991966804850723e-0 1 8.0332e-04 17

9.997991943200188e-0 1 2.0081e-04 33

9.999498000921015e-0 1 5.0200e-05 65

9.999874501175253e-0 1 1.2550e-05 129

9.999968625352869e-0 1 3.1375e-06 257

9.999992156341920e-0 1 7.8437e-07 513

...

9.999999999995609e-0 1 4.3909e-13 524289

1.000000000003847e+00 3.8467e-12 1048577

Fehlervierteltsihjeweils,unddasvonAnfangan.Wenigerals

10 ^{− 13}

^{wirdmitdouble}

Genauigkeitniht erreiht.

Summierte Simpsonregelfür

(i)

^.

I

^{Fehler #Fktausw.}

1.000134584974194e+00 1.3458e-04 5

1.000008295523968e+00 8.2955e-06 9

1.000000516684706e+00 5.1668e-07 17

1.000000032265001e+00 3.2265e-08 33

1.000000002016129e+00 2.0161e-09 65

1.000000000126001e+00 1.2600e-10 129

1.000000000007874e+00 7.8739e-12 257

1.000000000000491e+00 4.9094e-13 513

1.000000000000030e+00 2.9976e-14 1025

1.000000000000006e+00 5.7732e-15 2049

1.000000000000002e+00 1.7764e-15 4097

Fehler reduziert sih jeweilsum den Faktor

16 = ( ¹ / 2 ) ⁴

^{,und das fast von Anfangan.}

Summierte Trapezregelfür

(ii)

^.

I

^{Fehler #Fktausw.}

1.000009999900001e +05 9.9010e+04 3

5.000449983500645e +04 4.9013e+04 5

2.501113751079469e +04 2.4020e+04 9

1.252430268327760e +04 1.1533e+04 17

6.300548144658167e +03 5.3091e+03 33

3.227572909110977e +03 2.2361e+03 65

1.765586982280199e +03 7.7413e+02 129

1.160976493727309e +03 1.6952e+02 257

1.003813438906513e +03 1.2355e+01 513

9.915347090712996e +02 7.5876e-02 1025

9.914588358257512e +02 2.5795e-06 2049

9.914588331667655e +02 7.9478e-08 4097

9.914588332263689e +02 1.9875e-08 8193

9.914588332412698e +02 4.9740e-09 16385

Fehlerverhalten am Anfangunklar, erst spät stellt sih

h ²

^ein.

Summierte Simpsonregelfür

(ii)

^.

I

^{Fehler #Fktausw.}

3.333899978334190e +04 3.2348e+04 5

1.668001673605744e +04 1.5689e+04 9

8.362024407438566e +03 7.3706e+03 17

4.225963298451690e +03 3.2345e+03 33

2.203247830595247e +03 1.2118e+03 65

1.278258340003273e +03 2.8680e+02 129

9.594396642096787e +02 3.2019e+01 257

9.514257539662473e +02 4.0033e+01 513

9.874417991262286e +02 4.0170e+00 1025

9.914335447439017e +02 2.5289e-02 2049

9.914588322804369e +02 9.6581e-07 4097

9.914588332462367e +02 7.0486e-12 8193

9.914588332462367e +02 7.0486e-12 16385

Bis 4096 Auswertungen ist Simpson shlehter als Trapez. Asymptotishe Konver-

genzrate stellt sih erst für genügendkleines

h

^ein.

Summierte Trapezregelfür

(iii)

^.

I

^{Fehler #Fktausw.}

1.000000000000000e+00 5.7080e-01 3

1.366025403784439e+00 2.0477e-01 5

1.497854534051220e+00 7.2942e-02 9

1.544909572178587e+00 2.5887e-02 17

1.561626518913870e+00 9.1698e-03 33

1.567551211438566e+00 3.2451e-03 65

1.569648456389842e+00 1.1479e-03 129

1.570390396198308e+00 4.0593e-04 257

1.570652791478614e+00 1.4354e-04 513

1.570745576359828e+00 5.0750e-05 1025

1.570778383269506e+00 1.7944e-05 2049

1.570789982705718e+00 6.3441e-06 4097

1.570794083803873e+00 2.2430e-06 8193

1.570795533774854e+00 7.9302e-07 16385

Die Konvergenzordnung

h ²

^{wird niht erreiht, sondern nur ein}

h ^α

^mit

α < 2

^(siehe

unten).

Summierte Simpsonregelfür

(iii)

^.

I

^{Fehler #Fktausw.}

1.488033871712585e+00 8.2762e-02 5

1.541797577473481e+00 2.8999e-02 9

1.560594584887709e+00 1.0202e-02 17

1.567198834492299e+00 3.5975e-03 33

1.569526108946797e+00 1.2702e-03 65

1.570347538040268e+00 4.4879e-04 129

1.570637709467796e+00 1.5862e-04 257

1.570740256572051e+00 5.6070e-05 513

1.570776504653564e+00 1.9822e-05 1025

1.570789318906069e+00 7.0079e-06 2049

1.570793849184461e+00 2.4776e-06 4097

1.570795450836595e+00 8.7596e-07 8193

1.570796017098507e+00 3.0970e-07 16385

Die Simpsonregelzeigt dieselbe Konvergenzordnung wie dieTrapezregel!

Konvergenzabshätzung hat dieForm

| I (f ) − I _h ⁽ⁿ⁾ (f ) | ≤ Ch ^m+1 .

Experimentelle Bestimmung der Konvergenzrate:

Mit dem Ansatz

e h = | I(f) − I _h ⁽ⁿ⁾ (f ) | = Ch ^α

^gilt

e ^h / ²

e h

= C( ^h / ² ) ^α

Ch ^α = ( ¹ / ² ) ^α

und darauserhalten wir

α = log e ^h / ²

e h

log 1

2

.

Das sobestimmte

α

^heiÿt experimental order of onvergene (EOC).

5.2 Fehlerkontrolle

WirhabeninSatz6.4gezeigt,wiederFehlermitmehrStützstellenabnimmt.Diesnennt

man eine a-priori Fehlershranke.

So erhalten wir etwa für diesummierte Trapezregel:

| I(f) − I _h ⁽¹⁾ (f ) | = | − b − a

12 f ^′′ (ξ)h ² | ≤ b − a 12 max

ξ ∈ [a,b] | f ^′′ (ξ) |

| {z }

=:C

h ² .

C

^{ist allerdingsim}allgemeinen shwer zu bestimmen.

In der Praxis würde man aber gerne wissen, bei wievielen Stützstellen (bei welhem

h

^{)der Fehler kleiner alseine} vorgegebene Toleranz ist.

Dazu wollen wir eine Methode zur a-posteriori Fehlershätzung vorstellen.

Idee: Die Simpson-Summe konvergiert shneller als die Trapezsumme (

f

^genügend

glatt),hat also fürgenügend kleines

h

^{einen kleineren Fehler.}

Wir wollen den Fehler inder Trapezsumme zum Gitter

h/2

^{abshätzen. Dazu shie-}

ben wir dieAuswertung der Simpsonsumme dazwishen:

| I (f) − I ⁽¹⁾ h 2

(f) | = | I(f) − I _h ⁽²⁾ (f ) + I _h ⁽²⁾ (f ) − I ⁽¹⁾ h 2

(f ) | .

Nunnutze die Dreieksungleihung:

| I (f ) − I ⁽¹⁾ h 2

(f) | ≤ | I(f ) − I _h ⁽²⁾ (f ) | + | I _h ⁽²⁾ (f ) − I ⁽¹⁾ h 2

(f ) | .

Nunnimmtmanan,dassdieSimpsonsummegenaueristalsdieTrapezsumme:

| I(f ) − I _h ⁽²⁾ (f) | ≤ ω | I(f) − I ⁽¹⁾ h

2

(f) |

^mit

0 < ω < 1

^:

| I(f) − I ⁽¹⁾ h 2

(f ) | ≤ ω | I(f ) − I ⁽¹⁾ h 2

(f ) | + | I _h ⁽²⁾ (f) − I ⁽¹⁾ h 2

(f ) | .

Auösennah dem Fehler inder Trapezsumme liefert:

| I(f ) − I ⁽¹⁾ h 2

(f ) | ≤ 1

1 − ω | I _h ⁽²⁾ (f ) − I ⁽¹⁾ h 2

(f) | .

Algorithmishe Realisierung der Fehlerkontrolle

Besonders ökonomish lässt sih die Fehlerkontrolle realisieren, wenn man folgende

BeziehungenzwishenSimpson-Rgel,MittelpunktregelundTrapezregelbenutzt(deshalb

haben wir oben die Trapezregelzu

h/2

^{und die} Simpsonsumme zu

h

verwendet):

DiesummierteSimpson-RegellässtsihaussummierterTrapez-undMittelpunktregel

zusammensetzen:

I _h ⁽²⁾ (f ) = 1

3 I _h ⁽¹⁾ (f ) + 2

3 I _h ⁽⁰⁾ (f ).

Die summierteTrapezregel zu demnähst feineren Gittererhält manaus summierter

Trapez- und Mittelpunktregel des gröberen Gitters:

I ⁽¹⁾ h 2

= 1

2 I _h ⁽¹⁾ (f ) + 1

2 I _h ⁽⁰⁾ (f ) h = b − a

^;

N = 1

^;

I1 = h(f(a) + f(b))/2

^;

while (

h > ε

^{) do}

I0 = 0

^;

for (

i = 0

^,

i < N

^,

i = i + 1

^{) do}

I 0 = I0 + hf (a + (i + 0.5)h)

^; {Mittelpunktsumme}

endfor

I2 = ¹ ₃ I1 + ² ₃ I0

^;{Simpson-Summe zu

h

^}

I1 = ¹ ₂ I1 + ¹ ₂ I0

^;{Trapez-Summe zu

h/2

^}

h = ¹ ₂ h

^;

N = 2N

^;

if (

1

1 − ω (I2 − I1) ≤ T OL

^{) then}

return

I1

^{; {Liefere T}rapez-Summe zu

h

^}

endif

end while

Die Fehlerkontrolle lässt sih soohne zusätzlihen Aufwand erledigen.

5.3 Gauÿ- und adaptive Quadratur

Adaptive Quadratur

Quadratur mit konstanter Shrittweite ist bei manhen Integranden inezient. Be-

trahte z. B.

f(x) = _{10 −} ^{5 1} _+x ²

^.

0 20000 40000 60000 80000 100000

-1 -0.5 0 0.5 1

1/((1e-5)+(x**2))

In soeinem Fallmöhte man die Shrittweite adaptiv, d. h. angepasst an den spe-

ziellenIntegranden wählen.

Dazu bedient man sih eines lokalen Fehlershätzers (oder Indikators), der angibt,

anwelher Stelle dieShrittweite weiter verkleinert werden muss.

Esbietetsihan,dies nohmiteinerFehlerkontrolle zukombinieren.Grob ergibtsih

das folgende Vorgehen:

(1) Wähle eine Unterteilung

G 0 = { x ⁽⁰⁾ _i | i = 0, . . . , N 0 }

^{. Berehne das Integral}

I 0

bezüglih der Unterteilung

G 0

^{. Setze}

k = 0

^.

(2) Berehne eine Shätzung für den Fehler

E k

ⁱⁿ

I k

^.Falls

E k < T OL

^{sind wir fertig.}

(3) Verfeinere dieUnterteilung

G k

^zu

G k+1

^{durh HinzufügenvonPunkten angepasst}

anden Integranden und setze

k = k + 1

^.

(4) Berehne

I _k

^zu

G _k

^{und gehe nah(2).}

Prinzip von Arhimedes

Wir betrahten das Prinzip vonArhimedes 2

:

2

ArhimedesvonSyrakus,287v.Chr.-212v.Chr.,griehisherMathematiker,PhysikerundIngenieur.

00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000

11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 0000000000000000000

0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000

1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111

00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000 00000000000000000000000000000000000000

11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111 11111111111111111111111111111111111111

0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000 0000000000000000000

1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111 1111111111111111111

0 0 1 1

a b = x l

f(x)

I 0

I ₁

I 2 , negativ!

I 3

I 2 = ¹ ₂

f (x m ) − f (x l ) + f (x r ) 2

| {z }

g



x r − x m + x m − x l

| {z }

h





x m

x r

g h

1. Berehne rekursiv

I(f ) = I ₀ + I ₁ + I ₂ + I ₃ + . . .

2. Brehe rekursiven Astab, falls

| I j |

^kleingenug.

Beispiel 6.8 zur numerishen Quadratur

Vershiedene Quadraturen für

(i)

^{aus letztem Beispiel.}

Methode

I

^{Fehler #Fktausw.}

Gauss4 9.999101667698898e-01 8.9833e-05 4

9.999944679581383e-01 5.5320e-06 8

9.999996555171785e-01 3.4448e-07 16

9.999999784895880e-01 2.1510e-08 32

Gauss6 1.000000118910998e+00 1.1891e-07 6

1.000000001828737e+00 1.8287e-09 12

1.000000000028461e+00 2.8461e-11 24

1.000000000000444e+00 4.4409e-13 48

Arhi 9.480594489685199e-01 5.1941e-02 3

9.871158009727754e-01 1.2884e-02 5

9.967851718861697e-01 3.2148e-03 9

9.990131153231707e-01 9.8688e-04 15

9.997876171856270e-01 2.1238e-04 31

Das Verfahren hoherOrdnung zahlt sih aus.

Methode

I

^{Fehler #Fktausw.}

Trapez 1.000009999900001e+05 9.9010e+04 3

3.227572909110977e+03 2.2361e+03 65

1.765586982280199e+03 7.7413e+02 129

1.160976493727309e+03 1.6952e+02 257

1.003813438906513e+03 1.2355e+01 513

9.915347090712996e+02 7.5876e-02 1025

9.914588358257512e+02 2.5795e-06 2049

Arhi 1.767335925226728e+03 7.7588e+02 25

1.004348965298925e+03 1.2890e+01 37

9.946212584262852e+02 3.1624e+00 81

9.922788393957054e+02 8.2001e-01 173

9.916266302474447e+02 1.6780e-01 361

9.914967844457766e+02 3.7951e-02 769

9.914672523966888e+02 8.4192e-03 1625

9.914606991793892e+02 1.8659e-03 3465

9.914592092819358e+02 3.7604e-04 7629

Fehlerreduktion mit Arhi ist von Anfang an quadratish, allerdings überholt die

Methode

I

^{Fehler #Fktausw.}

Gauss4 1.592226038754547e+00 2.1430e-02 4

1.570801362699711e+00 5.0359e-06 1024

1.570798107100650e+00 1.7803e-06 2048

1.570796956200537e+00 6.2941e-07 4096

1.570796549318533e+00 2.2252e-07 8192

Gauss6 1.578036347519909e+00 7.2400e-03 6

1.570801237513435e+00 4.9107e-06 768

1.570798062869299e+00 1.7361e-06 1536

1.570796940567470e+00 6.1377e-07 3072

1.570796543792309e+00 2.1700e-07 6144

Arhi 1.366025403784439e+00 2.0477e-01 5

1.570774639679624e+00 2.1687e-05 365

1.570791453003758e+00 4.8738e-06 765

1.570795219591928e+00 1.1072e-06 1605

1.570796082320714e+00 2.4447e-07 3433

Hohe Ordnung lohnt sihniht wegen mangelnder Dierenzierbarkeit.

5.4 Mehrdimensionale Quadratur

Für Rehteke (

d = 2

^{), Quader(}

d = 3

^{), ...lassen sih die Formelnfür die} eindimensio- nale Integration leihterweitern:

Z d c

Z b a

f (x, y) dxdy ≈ Z d

c

X n

i=0

w i f (x i , y) dy

= X n

i=0

w i

Z d c

f (x i , y) dx

≈ X n

i=0

w i

X n

j=0

w j f(x i , y j )

!

= X n

i=0

X n

j=0

w i w j

| {z }

w ij

f (x i , y j )

Koordinatentransformation

Allerdings sind niht alle Gebiete Rehteke (anders als in 1D!). Dann lassen sih

Integrale mittelsKoordinatentransformationberehnen:

Z

Ω

f(x, y)dx dy = Z 1

− 1

Z 1

− 1

f

ϕ(ξ, η), ψ(ξ, η)

∂(ϕ, ψ)

∂(ξ, η) dξdη

wobei dieTransformation

ϕ(ξ, η) ψ(ξ, η)

: [ − 1, 1] × [ − 1, 1] → Ω

das Gebiet

[ − 1, 1] × [ − 1, 1]

^auf

Ω

^abbildet.

1

−1

−1 1

η

ξ

φ(ξ, η) ψ(ξ, η)

x

y

Weiter ist

∂ (ϕ, ψ)

∂(ξ, η) = det

∂ϕ

∂ξ (ξ, η) ^∂ψ _∂ξ (ξ, η)

∂ϕ

∂η (ξ, η) ^∂ψ _∂η (ξ, η)

!

6

= 0

dieDeterminanteder (transponierten) Jaobimatrix 3

der Transformation.

Es lassen sih auh Simplizes(=Dreiek, Tetraeder) zur Unterteilung verwenden und

fürdiese direkte Integrationsformelnherleiten.

Komplexe Gebiete

Bei komplexen Gebieten, z. B. Gebieten mit Löhern oder einer komplizierten Geo-

metriereiht das niht.

⇒

^{Zerlegung in}Teilgebiete, diesih aufRehteke transformieren lassen,aber:

•

Gittergenerierung niht trivialund shwierig automatish zu mahen.

•

^ErfordertBeshreibung des Gebietes

Ω

^.

•

Zusätzliher Geometriefehlerdurhniht exakte Darstellungdes Gebietes.

−1 −1 1

1

Ω

Auh inmehr als einer Raumdimensionenkann man hierarhishadaptiv verfeinern:

Links:AdaptivesDreieks- und Viereksgitter.Rehts: Adaptives Tetraedergitter mit

Bisektionsverfeinerung.

Fluh der Dimension

•

^Ist

d

^{sehr groÿ,so sind diehier} behandelten Methoden nihtbrauhbar.

•

^Betrahte

Ω = [0, 1] ^d

^{. Zerlegt man}

[0, 1]

^{in zwei T}eilintervalle je Rihtung, so hat man den d-dimensionalen Würfel in

2 ^d

^{Teilwürfel zerlegt.}

⇒

^{Der Aufwand steigt}

exponentiellin

d

^{an. Dies bezeihnet man alsFlux der Dimension.}

•

^{Eine Möglihkeitist dann die}Monte-CarloIntegration

I(f) ≈ C N

X N

i=1

f (ξ i )

mitZufallszahlen

ξ i ∈ Ω

^.

6 Ein kleiner Programmierkurs

6.1 Hallo Welt

Programmierumgebung

•

^{Wirbenutzen die}Programmiersprahe C++.

•

^{Wirbehandeln nurdie} Programmierungunter LINUXmitden GNUCompilern.

•

^{Windows: On your own.}

•

^Wirsetzen Grundfertigkeitim Umgang mitLINUX-Rehnern voraus:

Shell, Kommandozeile,Starten von Programmen.

Dateien, Navigieren im Dateisystem.

Erstellen von Textdateien miteinem Editor ihrer Wahl.

•

^{Ideedes Kurses: Lernen anBeispielen,keine rigorose} Darstellung.

•

^{BlutigeAnfänger solltenzusätzlih einBuhlesen (siehe} Literaturliste).

Workow

C++ ist eine kompilierte Sprahe. Um ein Programm zur Ausführung zu bringen

sind folgendeShritte notwendig:

1. Erstelle/Ändereden Programmtext miteinemEditor.

2. Übersetze den Programmtextmitdem C++-Übersetzer(auhC++-Compiler)

ineinMashinenprogramm.

3. Führedas Programmaus. Das Programm gibt sein Ergebnis auf dem Bildshirm

oder ineine Dateiaus.

4. Interpretiere Ergebnisse. Dazu benutzen wir weitere Programme wie gnuplot.