Wissenschaftliches Programmieren: Numerische Rechnungen

Wissenschaftliches Programmieren:

Numerische Rechnungen

Jan Piclum

Wintersemester 2020/21

Inhalt

Nullstellenbestimmung

Intervallhalbierung und Regula falsi Newton-Raphson-Methode

Numerische Differentiation Numerische Integration

Newton-Cotes-Regeln Uneigentliche Integrale Gauß-Quadratur Monte-Carlo-Integration

Differentialgleichungen

Euler-Verfahren Runge-Kutta-Verfahren

Nützliche Bibliotheken

Bestimmung von Nullstellen

Lösung von nichtlinearen Gleichungen iterative Verfahren:

Intervallhalbierung Regula falsi

Newton-Raphson-Methode . . .

Newton-Raphson-Methode kann auch auf Gleichungssysteme angewandt werden

Voraussetzung: Intervall mit mindestens einer Nullstelle ist bekannt

→problematisch im mehrdimensionalen Fall Ergebnis: Nullstelle mit vorgegebener Genauigkeitǫ

Vorsicht: Sprungstellen und Singularitäten können nicht immer von Nullstellen unterschieden werden

Iterative Verfahren

bei iterativen Verfahren wird eine Funktionsvorschrift wiederholt angewandt xn+1=f(xn), xn = (f◦f◦ · · · ◦f

| {z }

n mal

)(x0)

gesucht ist ein Fixpunktx^∗=f(x^∗)

die Eigenschaften folgen aus dem Banachschen Fixpunktsatz:

Für eine AbbildungF eines Banachraumes mit

||F(x)−F(y)||< L||x−y||, L <1gilt:

1 F hat genau einen Fixpunktx^∗.

2 Die Folgexn+1=F(xn) konvergiert für jeden Startwert gegenx^∗.

3 Fehlerabschätzung:||xn−x^∗|| ≤ Lⁿ

1−L||x1−x0||

in der Praxis hängt die Konvergenz häufig vom Startwert ab

Intervallhalbierung

Algorithmus:

gegeben:[x0, y0]mitf(x0)f(y0)<0 berechne:zn+1=xn+yn

2

istf(zn+1)f(xn)<0, setzexn+1=xn,yn+1=zn+1, sonstxn+1=zn+1,yn+1=yn

wiederhole bis Genauigkeityn+1−xn+1< ǫI oder|f(zn+1)|< ǫf erreicht ist

Eigenschaften:

findet immergenaueine Nullstelle (oder Singularität) weitere Nullstellen im

Startintervall werden übersehen konvergiert linear

-6 -4 -2 0 2

-20 -10 0 10 20

x0 x1 y0

y1

Regula falsi

verbessere Wahl der neuen Intervallgrenze durch lineare Interpolation:

f˜(x) = x−xn

yn−xnf(yn) + x−yn

xn−ynf(xn)

wähle Nullstelle vonf˜:zn+1= xnf(yn)−ynf(xn) f(yn)−f(xn)

Eigenschaften:

Konvergenz ist typischerweise besser als linear

Konvergenz ist schlecht, wenn sich die Funktion nahe der Nullstelle stark ändert

0 2 4 6 8 10

-3 -2 -1 0 1 2

x0

x1

y0

y1

verbesserte Konvergenz:

ändert sichxn [yn] mehrmals hintereinander, benutze beim nächsten Schritt f(yn+1) =f(yn)/2 [f(xn+1) =f(xn)/2]

Newton-Raphson-Methode 1

Algorithmus:

gegeben Startwertx0, sowief(x)undf^′(x)

Iteration: berechne Nullstelle der Tangente anxn und setze xn+1=xn− f(xn)

f^′(xn)

wiederhole bis Genauigkeit

f(xn+1) f^′(xn+1)

< ǫerreicht ist

→ suche eines Fixpunktesx^∗ vonF(x) =x− f(x) f^′(x) es giltF^′(x^∗) = f(x^∗)f^′′(x^∗)

[f^′(x^∗)]² = 0 ⇒superstabiler Fixpunkt

Newton-Raphson-Methode 2

Eigenschaften:

in der Näheder Nullstelle ist die Konvergenz quadratisch

problematisch sind zu weit entfernte Startwerte oder Stellen mitf^′(xn)≈0

f^′ sollte analytisch bekannt sein _-1

0 1 2 3 4 5

-20 -10 0 10 20 30 40

x0

x1

x2

wird häufig benutzt, um Nullstelle zu verbessern, die mit anderen Verfahren gefunden wurden

weitere Ableitungen können auch berücksichtigt werden

→Halley-Methode

Newton-Raphson in höheren Dimensionen

F~(~x) =~0 Iterationsschritt:

~xn+1=~xn−J⁻¹(~xn)·F~(~xn), J_ij(x~n) = ∂Fi

∂xj

~ x=~xn

Konvergenz kann verbessert werden, indem der zweite Term mit einem geeigneten Faktorλ∈(0,1]multipliziert wird

(sieheNumerical Recipes, Abschnitt 9.7) Eigenschaften:

erfordert Matrixinversion

Nullstellen können nicht a priori eingeklammert werden

genaue Analyse des Systems vor Anwendung der Methode ist erforderlich, um geeignete Startwerte zu finden

Numerische Differentiation

Differenzenquotient ist ein kontinuierlicher Grenzwert:

f^′(x) = lim

h→0

f(x+h)−f(x) h numerische Auswertung erfordert Diskretisierung

benutze Taylorentwicklung für kleine h:

f(x±h) =f(x)±f^′(x)h+1

2f^′′(x)h²±1

6f⁽³⁾(x)h³+O(h⁴)

→ erste Ableitung:f^′(x)≈f(x+h)−f(x)

h +O(h)

Erste Ableitung: Vorwärtsableitung

f^′(x)≈ f(x+h)−f(x)

h +O(h)

numerische Eigenschaften:

Abbruchfehler ist von Größenordnungh

Rundungsfehler durch Auslöschung treten bei zu kleinemhauf typischerweise existiert ein Bereich, in dem das Ergebnis stabil ist Wahl vonh:

(x+h)−xmussals Gleitkommazahl darstellbar sein:

1 tmp = x + h 2 h = tmp - x

optimale Wahl:h∼√ǫfxK,

mit Genauigkeit der Funktionǫf und typischer KrümmungxK =p f /f^′′

Erste Ableitung: symmetrische Zweipunkt-Formel

Abbruchfehler kann durch symmetrische Ableitung verbessert werden:

f^′(x)≈ f(x+h)−f(x−h)

2h +O(h²)

Eigenschaften:

erfordert eventuell eine weitere Auswertung der Funktion optimale Wahl:h∼ǫ^1/3_f xK

mitdouble precisiontypischerweise zwei Größenordnungen genauer als Vorwärtsableitung

weiter Verbesserungen sind möglich (siehe Übungen), aber erfordern mehr Funktionsauswertungen

sind anfälliger für Auslöschung

Zweite Ableitung

Herleitung:

benutze wieder Taylorentwicklung

f(x+h) +f(x−h) = 2f(x) +f^′′(x)h²+O(h⁴) oder leite Formel für erste Ableitung ab

f^′′(x) = f(x+h)−2f(x) +f(x−h)

h² +O(h²)

höhere Ableitungen können analog hergeleitet werden

Numerische Integration

Integrale sind Grenzwerte von Riemannsummen:

Z b a

dx f(x) = lim

n→∞

h→0

Xn

k=1

h f(xk), xk =a+k h , h= b−a n

Diskretisierung durch Zerlegung innTeilintervalle:

Z b a

dx f(x) =

n−1

X

k=0

Z xk+h xk

dx f(x)

und Taylorentwicklung der Integranden mit diskretisierten Ableitungen

→Newton-Cotes-Regeln

Stützstellen müssen nicht äquidistant gewählt werden→Gauß-Quadratur mehrdimensionale Integrale:

Abfolge von eindimensionalen Integrationen Monte-Carlo Integration

Trapezregel 1

Z b a

dx f(x) =

n−1

X

k=0

Z xk+h xk

dx f(x)

entwicklef(x)bis zur ersten Ordnung um den Randxk des Intervalls:

f(x) = f(xk) +f^′(xk) (x−xk) +O[(x−xk)²]

= f(xk) +f(xk+h)−f(xk)

h (x−xk) +O[(x−xk)²,(x−xk)h]

integriere über das Intervall:

Z xk+h xk

dx f(x) = h f(xk) +f(xk+h)−f(xk) h

h²

2 +O(h³)

= h

2[f(xk+h) +f(xk)] +O(h³)

Trapezregel 2

für das gesamte Integral ergibt sich Z b

a

dx f(x) =

n−1

X

k=0

Z xk+h xk

dx f(x) =

n−1

X

k=0

h

2[f(xk+h) +f(xk)] +O(nh³)

= h

2[f(a) +f(b)] +h

n−1

X

i=1

f(a+ih) +O(1/n²)

Eigenschaften:

Integrand wird durch Trapeze angenähert

abgeschlossene Regel:

Auswertung auch an den Rändern

Gewichte:

1/2an den Rändern

1im Inneren ^{0 1 2 3 4 5 6}

0 10 20 30

x0 x1 x2 x3 x4

Trapezregel 3

Eine genauere Fehlerabschätzung ist mit der Euler-McLaurin-Formel möglich:

n−1X

k=1

g(k) = Z n

0

dk g(k)−1

2[g(0) +g(n)] +X

i≥1

B2i

(2i)!

g⁽²ⁱ⁻¹⁾(n)−g⁽²ⁱ⁻¹⁾(0)

mit den Bernoullizahlen:B2= 1/6,B4=−1/30, . . . Wir haben also (mitg(k) =f(xk)unddx=hdk):

Z b a

dx f(x) =h

n−1X

k=1

f(xk) +h

2[f(a) +f(b)]−h² 12

f^′(b)−f^′(a) +. . .

Es gilt:

Alle Fehlerterme sind gerade.

Der Fehler beträgt betragsmäßig höchstens das Doppelte des ersten vernachässigten Terms.

Mittelpunktsregel

entwicklef(x)bis zur ersten Ordnung um den Mittelpunktyk=xk+h/2:

Z xk+h xk

dx f(x)

=

Z xk+h xk

dx

f(yk) +f(yk+h)−f(yk)

h (x−yk) +. . .

= h f(yk) +O(h³) =h f xk+h

2

+O(h³) insgesamt ergibt sich:

Z b a

dx f(x) =

n−1X

k=0

h f xk+h

2

+O(nh³) =h

n−1X

k=0

f

a+(2k+ 1)h 2

+O(1/n²) Eigenschaften:

Integrand wird durch Rechtecke angenähert offene Regel: keine Auswertung an den Rändern Gewichte: alle1

auch hier sind alle Fehlerterme gerade

Simpsonregel 1

entwickle Integranden bis zurzweitenOrdnung um die Intervallmitte Z xk+h

xk

dx f(x)

=

Z xk+h xk

dx

f(yk) +· · ·+ 2f(xk+h)−2f(yk) +f(xk)

h² (x−yk)²+. . .

= h f(yk) + 2f(xk+h)−2f(yk) +f(xk) h²

2 3

_h 2

³

+O(h⁵)

= h

6 h

f(xk+h) + 4f xk+h

2

+f(xk)i

+O(h⁵)

benutzen^′=n/2 Intervalle undh= 2h^′ (yk=xk+h^′,xk =a+ 2kh^′):

Z xk+h xk

dx f(x) =h^′ h₁

3f(yk+h^′) + 4

3f(yk) +1

3f(yk−h^′)i

+O(h^′5)

Simpsonregel 2

insgesamt ergibt sich (yk=xk+h^′,xk=a+ 2kh^′):

Z b a

dx f(x)

=

n^′−1

X

k=0

h^′h₁

3f(yk+h^′) +4

3f(yk) +1

3f(yk−h^′)i

+O(n^′h^′5)

= h^′

3 [f(a) +f(b)] +4h^′ 3

n^′−1

X

k=0

f(a+ (2k+ 1)h^′) +2h^′ 3

n^′−1

X

k=1

f(a+ 2kh^′) +O(1/n^′4)

Eigenschaften:

Integrand wird quadratisch angenähert

Regel ist exakt für quadratische und kubische Funktionen abgeschlossene Regel

Gewichte:

1/3an den Ränderen

abwechselnd4/3und2/3im Inneren

Simpsonregel 3

alternative Herleitung der Simpsonregel:

Es seiTN das Ergebnis der Trapezregel mitN Teilintervallen.

Dann gilt für den Fehler vonTN:δN ∼1/N² Bei Halbierung der Intervalle folgt somit:

δN =1

4δN/2+O(1/N⁴). Also hat die Kombination

SN = 4

3TN −1 3TN/2

einen Fehler der OrdnungO(1/N⁴).

→Beispiel der Romberg-Integration

Iterierte Trapezregel 1

Die Trapezregel kann auch als iteratives Verfahren implementiert werden:

1. Schritt:h1=b−a Z b

a

dx f(x)≈(b−a)h₁

2f(a) +1 2f(b)i

≡T1

2. Schritt:h2= (b−a)/2→erfordert eine zusätzliche Stützstelle Z b

a

dx f(x) ≈ b−a 2

h₁

2f(a) +f_a+b 2

+1

2f(b)i

≡T2

T2 = 1 2 h

T1+ (b−a)f_a+b 2

i

(n+ 1). Schritt: hn+1=hn/2 = (b−a)/2ⁿ

→2ⁿ⁻¹ zusätzliche Stützstellen in der Mitte aller Intervalle aus Schrittn

Tn+1=1 2



Tn+hn 2ⁿ⁻¹−1

X

k=0

f(a+khn+hn+1)





Iterierte Trapezregel 2

Eigenschaften:

in jedem Schritt wird das Ergebnis des letzten Schritts wieder verwendet es müssen nur die neuen Stützstellen berechnet werden

Konvergenz und Fehler können kontrolliert werden:

→Iteration bis eine gewünschte Genauigkeit erreicht ist kann mit Romberg-Integration kombiniert werden:

interpoliere Ergebnisse mit Polynom(n−1)-ten Grades inh² Wert des Polynoms beih= 0 approximiert das Integral

Die Mittelpunktsregel kann auch iteriert werden, wenn die Breite der Teilintervalle in jedem Schritt durch 3 geteilt wird.→Übungen

Unendliches Integrationsgebiet

Integrale mit unbeschränktem Integrationsgebiet:

numerische Integration wird bei endlichemxmax abgebrochen Z ∞

a

dx f(x) = Z xmax

a

dx f(x) + Z ∞

xmax

dx f(x) und der Fehler durch den letzten Term abgeschätzt

Substitution auf endliches Gebiet zum Beispiel mitx= tan(u)oder x= 1/u

→Jacobideterminante kann singulär sein Z ∞

a

dx f(x) = Z 1/a

0

du 1 u² f₁

u

Singuläre Integranden

Betrachte Integral mit Singularität am Rand:

vermeide Auswertung am Rand zum Beispiel durch Mittelpunktsregel Z b

0

dxsin(x)

x → f(0)kann nicht einfach berechnet werden

spalte Singularität ab:

Z b

0

dxf(x) x^α =

Z b

0

dxf(x)−f(0)

x^α +

Z b

0

dxf(0) x^α mit0< α <1 undf(x)regulär beix= 0

→der letzte Term kann jetzt analytisch integriert werden

Gauß-Quadratur

Wir betrachten Integrale der Form I=

Z b a

dx W(x)f(x)≈

N−1X

j=0

wjf(xj) mit einerGewichtsfunktionW(x).

Wähle dieN Stützstellen so, dassI exakt ist, wennf ein Polynom vom Grad 2N−1 oder weniger ist.

Die Gewichtewj sind dann die Nullstellen von bestimmten orthogonalen Polynomen.

Eigenschaften:

konvergiert exponentiell mit wachsendemN, wennf glatt ist Wahl der Gewichtsfunktion kann integrable Singularitäten entfernen

Beispiele

Gauß-Legendre:

W(x) = 1, −1< x <1

→Legendre Polynome erster Art

→Basis der Gauß-Konrod-Quadratur Gauß-Tschebyschow:

W(x) = 1

√1−x², −1< x <1

→Tschebyschow Polynome Gauß-Laguerre:

W(x) =e^−x, 0< x <∞

→Laguerre Polynome

Monte-Carlo-Integration

statistisches Verfahren:

Integrand wird an zufälligen Stellen ausgewertet und darüber gemittelt

→sampling

Fehler ist nicht von der Dimension abhängig

→mehrdimensionale Integrale können effizient berechnet werden

Monte-Carlo-Methode wird auch zur Simulation von physikalischen Prozessen verwendet:

simulated annealing Metropolis-Algorithmus . . .

Eindimensionales Beispiel

Bestimme die FlächeF unter der Kurve vong zwischenaund b.

Wähle einebekannteFläche A, dieF vollständig einschliesst.

ErzeugeNA Paare von gleichverteilten Zufallszahlen (xi, yi)∈A.

Bestimme die AnzahlNF aller Paare für die gilt

xi ∈[a, b]und0≤yi≤g(xi).

Dann ist:F ≈ANF/NA. Der statistische Fehler ist dabei proportional zu1/√

NA.

0 1 2 3 4 5 6

-1 0 1 2 3

a b

Eindimensionale Integrale

Wir wählen zufällige Stützstellen, die in[a, b]gleichverteilt sind:

p(x) = ₁

b−a x∈[a, b]

0 sonst

Dann haben wir Z b

a

dx g(x) = Z

dx p(x) (b−a)g(x)

| {z }

=f(x)

=hfi^p±σI

Der Mittelwerthfip wird durch N Stichproben approximiert, wobei der statistische Fehler durch die Varianz gegeben ist:

hfi^p = 1 N

XN

i=1

f(xi) = b−a N

XN

i=1

g(xi),

σI² = 1

N hf²i^p− hfi²p

∼ 1 N

Mehrdimensionale Integrale

Der multidimensionale Fall wird analog hergeleitet:

Z

V

dⁿ~x g(~x) = Z

dⁿ~x p(~x)f(~x) =hfip±σI

p(~x)f(~x) =

g(~x) ~x∈V

0 sonst

Der statistische Fehler ist weiterhinσI = 1

√N q

hf²i^p− hfi²p∼ 1

√N

→unabhängig von der Dimension Vorsicht:

Der Fehler stammt aus der Varianz desMittelwerts.

Es handelt sichnichtum eine strenge Fehlerabschätzung des Integrals.

Der Fehler mussnichtGauß-verteilt sein.

Varianzreduktion

Importance Sampling:

es können auch andere Verteilungen benutzt werden reduziere statistischen Fehler durch Anpassen vonpang

σ²_I= 1 N

g² p²

− g

p 2!

optimale Wahl:

p(~x) =

c|g(~x)| ~x∈V

0 sonst , 1

c = Z

V

dⁿ~x|g(~x)|

Stratified Sampling:

teile Integrationsvolumen in Teilvolumen auf

die Anzahl der Punkte in den Teilvolumen kann unterschiedlich sein

→wird zum Beispiel inVEGASbenutzt

Gewöhnliche Differentialgleichungen

Wir betrachten gewöhnliche Differentialgleichungenk. Ordnung:

y^(k)(x) =g

x, y(x), y^′(x), . . . , y^(k−1)(x)

(1) Wir können die Ableitungen vony als Komponenten eines Vektors schreiben:

~y(x) =







y1 ≡ y(x) y2 ≡ y^′(x)

...

yk ≡ y^(k−1)(x)







Dann kann (1) als System vonk Differentialgleichungen 1. Ordnung geschrieben werden:

~y^′=







 y^′₁ y^′₂ ... y^′_k−1

y^′_k









=









y2

y3

... yk

g(x, y1, y2, . . . , yk)









≡f~(x, ~y)

Euler-Verfahren 1: Rekursionsgleichung

Wir betrachten die zeitabhängige DGl. 1. Ordnung:

~y˙=f~(t, ~y), ~y(0) =~y0.

Wir suchen die Lösung fürt∈[0, T].

Diskretisiere das Zeitintervall mit Schrittweiteh:

tn=nh , n= 0,1, . . . , N , N =T /h Gesucht sind dann die Werte~yn =~y(tn).

Benutze diskretisierte Ableitung:

~˙

y= ~yn+1−~yn

h +O(h) =f~(tn, ~yn) +O(h)

→Rekursionsgleichung:

~yn+1=~yn+h ~f(tn, ~yn) +O(h²)

Euler-Verfahren 2: Eigenschaften

~yn+1=~yn+h ~f(tn, ~yn) +O(h²)

Eigenschaften:

Fehler füreinenSchritt: O(h²)

akkumulierter Fehler fürN Schritte:O(N h²) =O(h)

→Verfahren 1. Ordnung

Einschrittverfahren:~yn+1 wird nur aus~yn bestimmt Genauigkeit setzt Stetigkeit der Funktion voraus

→Kann die Genauigkeit erhöht werden?

Euler-Verfahren 3: Verbesserung der Genauigkeit

Wir integrieren die DGl. über einen Zeitschritt:

~yn+1=~yn+ Z tn+1

tn

dt ~f(t, ~y(t))

Das Euler-Verfahren entspricht somit der Näherung Z tn+1

tn

dt ~f(t, ~y(t)) =h ~f(tn, ~yn) +O(h²)

Können wir stattdessen die Trapezregel benutzen?

Z tn+1

tn

dt ~f(t, ~y(t)) = h 2

hf~(tn, ~yn) +f~(tn+1, ~yn+1)i

+O(h³)

Problem:~yn+1steht in Schritt nauf beiden Seiten der Rekursionsgleichung.

Euler-Verfahren 4: Prädiktor-Korrektor-Verfahren

Lösung:

Vorhersage von~yn+1 mit einfachem Euler-Verfahren:

~k1 = h f(tn, ~yn)

~yn+1,P = ~yn+~k1

Korrektur der Vorhersage mit

~k2 = h f(tn+1, ~yn+1,P)

~yn+1 = ~yn+1 2

h~k1+~k2

i+O(h³)

→Heun-Verfahren 2. Ordnung:

Prädiktor der OrdnungO(h²)reicht für Korrektor der OrdnungO(h³) akkumulierter Fehler nachN Schritten:O(h²)

Verfahren ist numerisch stabiler

f~muss in jedem Schritt zweimal ausgewertet werden

Runge-Kutta-Verfahren

benutze andere Newton-Coates-Regeln:

Mittelpunktsregel→Runge-Kutta-Verfahren 2. Ordnung Simpsonregel→Runge-Kutta-Verfahren 4. Ordnung

Eigenschaften des Verfahrens 4. Ordnung:

gute Genauigkeit

geeignet für Bewegungsgleichungen mit wenigen Freiheitsgraden kann durch Schrittweitenanpassung verbessert werden

Runge-Kutta-Verfahren 2. Ordnung

Wir benutzen Z tn+1

tn

dt ~f(t, ~y(t)) =h ~f(tn+1/2, ~yn+1/2) +O(h³), tn+1/2=1

2(tn+tn+1) Der Zwischenschritt wird mit dem Euler-Verfahren berechnet.

Wir erhalten:

~k1 = h ~f(tn, ~yn)

~k2 = h ~f

tn+1/2, ~yn+1 2~k1

~yn+1 = ~yn+~k2+O(h³)

Der Prädiktor wird nur für einen halben Schritt vorhergesagt.

Runge-Kutta-Verfahren 4. Ordnung

Wir benutzen Z tn+1

tn

dt ~f(t, ~y(t)) = h 6

hf~(tn, ~yn) + 4f~(tn+1/2, ~yn+1/2) +f~(tn+1, ~yn+1)i

Die Rekursionsgleichung benötigt Vorhersagen für zwei Stellen:

~yn+1 = ~yn+h 6

hf~(tn, ~yn) + 2f~(tn+1/2, ~yn+1/2)

+ 2f~(tn+1/2, ~yn+1/2) +f~(tn+1, ~yn+1)i Diese erfolgen in vier Schritten:

~k1 = h ~f(tn, ~yn)

~k2 = h ~f

tn+1/2, ~yn+1 2~k1

~k3 = h ~f

tn+1/2, ~yn+1 2~k2

~k4 = h ~f(tn+1, ~yn+~k3)

~yn+1 = ~yn+1 6

h~k1+ 2~k2+ 2~k3+~k4

i+O(h⁵)

Runge-Kutta-Verfahren 4. Ordnung

0 1 2 3 4

0.0 0.5 1.0 1.5 2.0 2.5 3.0

0 1 2 3 4

0.0 0.5 1.0 1.5 2.0 2.5 3.0

~ki/happroximiert die Steigungen an den verschiedenen Stellen:

zur Berechnung von~k2 wird~yn+1/2mit Hilfe von~k1approximiert zur Berechnung von~k3 wird~yn+1/2mit Hilfe von~k2approximiert zur Berechnung von~k4 wird~yn+1 mit Hilfe von~k3approximiert

Nützliche Bibliotheken

GSL – GNU Scientific Library

→CundC++Bibliothek für numerische Rechnungen Boost

→Sammlung vonC++Bibliotheken Cuba

→Bibliothek für numerische Integration GiNaC

→Bibliothek für symbolische Rechnungen