Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik ::

Numerische Algorithemen und Mathematik

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik ::

Numerische Algorithemen und Mathematik

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik ::

Numerische Algorithemen und Mathematik

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Mathematik

Mathematik

Algorithmen

Funktionen

Datenstrukturen

Funktionen (!), Vektoren, Matrizen, Tensoren usw.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Künstliche Intelligenz

Künstliche Intelligenz

Algorithmen

Iterative Trainingsalgorithmen, Vereinfachung

Datenstrukturen

Funktionen, Datengraphen und Bäume, Funktionsgraphen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik

Numerik

numalg Beziehung der Numerik zu anderen Bereichen:

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik

Numerik

Das Ziel der Numerik ist die Konstruktion ökonomischer (effizienter) und stabiler Algorithmen. Speziell gilt es, mögliche Fehlerquellen zu berücksichtigen. Diese ergeben sich durch Modellierungsfehler, durch Fehler in den Eingangsdaten, durch Fehler im Algorithmus, und durch Diskretisierungsfehler in der Numerik.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Anwendungen und Probleme

Numerik: Anwendungen und Probleme

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Anwendungen und Probleme

Numerik: Anwendungen und Probleme

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Anwendungen und Probleme

Numerik: Anwendungen und Probleme

• Wir werden einige überraschende Ergebnisse sehen und dass Numerik (und die Ergebnisse daraus) ähnlich verwirrend und unerwartet sein können wie bei der Betriebssystemprogrammierung!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Fehleranalyse

Numerik: Fehleranalyse

Wir haben bei numerischen Problemen folgende Randbedingungen zu berücksichtigen:

1. Kondition
2. Rundungsfehler
3. Stabilität

algomat Beim Ausführen von Rechnungen gibt es verschiedene Fehlerquellen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik: Fehleranalyse

Numerik: Fehleranalyse

Kondition eines Problems

Die Kondition eines Problems gibt an, welche Genauigkeit man bei exakter Lösung bestenfalls erreichen kann: Ein Problem heißt gut konditioniert, wenn kleine Störungen der Eingangsdaten kleine Änderungen im Ergebnis bewirken, sonst schlecht konditioniert.

Rundungsfehler

In der Mathematik werden kontinuierliche und gar unendlich große Zahlenmenge betrachtet (reele Zahlen). In der Numerik haben wir immer diskerete und begrenzte Wertemengen!. Daher sind numerische Verfahren immer ungenauer als mathematische.

Stabilität

Iterative Algorithmen können eine "richtige" Lösung liefern (Konvergenz), können aber auch falsche Ergebnisse liefern (Divergenz). Man versucht Algorithmen stabil gegen Divergenz zu machen.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Numerik und Zahlensysteme

Numerik und Zahlensysteme

Mathematisch können wir zwei Basismengen von Zahlen unterscheiden:

Ganze Zahlen ℕ

Eine unendlich große Menge ganzer positiver und negativer Zahlen inklusive 0. Aber die Menge ist abzählbar! Es gibt keine weitere ganze Zahl in jedem Intervall [x,x+1].

Reelle Zahlen ℝ

Eine unendlich große Menge kontinuierlicher positiver und negativer Zahlen inklusive 0. Aber die Menge ist nicht abzählbar! Es gibt unendlich viele weitere Zahlen in jedem Intervall [x,x+δ], egal wie klein man das Intervall macht.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Computer Zahlensysteme

Computer Zahlensysteme

Maschinenzahlen und Rundungsfehler

Maschinenzahlen werden durch Datenbits repräsentiert (kodiert).

Ganze Zahlen (Integer Datentyp)

Eine endlich große Menge ganzer positiver und negativer Zahlen inklusive 0. Es gibt keine weitere ganze Zahl in jedem Intervall [x,x+1].

Festpunktzahlen (Fixed Point Datentyp)

Eine endlich große Menge diskreter positiver und negativer Zahlen inklusive 0. Aber es gibt keine weiter Zahl in einem kleinen Intervall [x,x+δ_min], δ_min ist absolut und hängt von der Anzahl der Datenbits und der absoluten Größe des Zahlenintervalls ab! Eigentlich ganze Zahlen mit einem festen verschobenen Dezimalpunkt.

Gleitkommazahlen (Floating Point Datentyp)

Eine endlich große Menge diskreter positiver und negativer Zahlen inklusive 0. Aber es gibt keine weiter Zahl in einem kleinen Intervall [x,x+δ_min], δ_min ist relativ und hängt von der Anzahl der Datenbits ab!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Computer Zahlensysteme

Computer Zahlensysteme

Integer Zahlen

Die digitalen Zahlen (Kodierung) d kann man in dezimale mit der gewichteten Summe umrechnen (Binärzahlensystem mit Basis b=2):

$${a}={\left(-{1}\right)}^{{s}}{\sum_{{{i}={0}}}^{{{m}-{1}}}}{d}_{{i}}\cdot{b}^{{i}}$$

algomat Zahlenbereich hängt von der Anzahl Bits ab (m). Für das Vorzeichen kommt noch ein Bit hinzu (oder der Zahlenbereich verringert sich). Es gibt keinen Rundungsfehler!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Arithmetik

Arithmetik

Es gibt die grundlegenden numerischen arithmetischen Operationen:

1. Addition (Subtraktion)
2. Multiplikation (Division)
3. Negierung

a=1; b=10; c=100;
x1 = (a/b)*c
x2 = (a*c)/b
print(x1,x2)

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Arithmetik

Das Experiment (Hinweis: JavaScript verarbeitet alle Zahlen als Gleitkommazahlen, daher ist Umwandlung in Integer erforderlich)

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Arithmetik

Gleitkommazahlen

$${a}={\left(-{1}\right)}^{{s}}{\sum_{{{i}={0}}}^{{{m}-{1}}}}{d}_{{i}}\cdot{b}^{{-{i}+{e}}}$$

Eigenschaften verschiedener Datenformate bei Gleitkommazahlen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Rundungsfehler

Rundungsfehler

• oder richtig ausgedrückt Diskretisierungsfehler bezüglich des Zahlensystems und deren Kodierung (da niemand explizit rundet)

• Rundungsfehler, Über- und Unterlauf und "Not a Number NaN" können bei arithmetischen Operationen auftreten.

Diskretisierung

• Die Dikretisierung von Integer Zahlen ist ohne Genauigkeitsverlust möglich, nur die Wertemenge ist begrenzt

• Die Diskretisierung von reelen Zahlen ist i.A. immer mit Genauigkeitsverlust und einer Einschränkung der Wertemenge verbunden (also schon hier RUndungsfehler möglich)

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Rundungsfehler

Rundungsfehler

Arithmetik

• Arithmetische Operationen wie Division können bei ganzen Zahlen zu erheblichen Rundungsfehlern führen (bis zu 100%), Addition und Multiplikation führen keine weiteren Rundungsfehler ein, können aber zum Überlauf bei Integer Zahlen führen!

• Arithmetische Operationen können bei Gleitkommazahlen (aber nicht reellen) zu Rundungsfehlern führen (bis zu 1%), und es kann zum Überlauf kommen!

• Im Bereich der ML Algorithmen spricht man vom verschwindenden Gradienten (auch ein Quotient mit Divisionsoperation)
• Bei Gleitpunktarithmetik gibt es eine Fehlerverstärkung bei den elementaren Rechenoperationen!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Funktionen und Schleifen

Funktionen und Schleifen

• Funktionsrekursion ist eine gängige Methode in der Mathematik um eine berechung auszudrücken.
• Programmatisch können wir in fast allen Programmiersprachen die Funktionrekursion nutzen, es gibt aber auch Nachteile (welche?)
• Man kann rekursive Funktionen in Schleifen transformieren und umgekehrt!

function fac(n) {
  if (n<2) return 1
  else return n*fac(n-1)
}

y=1
for(i=2;i<=n;i++) {
  y=y*i
}

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

numalg pp 6

1. Grenzwert

$${e}=\lim_{{{n}\to\infty}}{\left({1}+\frac{{x}}{{n}}\right)}^{{n}}$$

für x=1 berechnen.

⇒ Grenzwerte können nicht direkt prozedural (iterativ) algorithmisch gelöst werden (deklarative und symbolische Methodik erforderlich). Nur endliche Approximation möglich!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

n=100
e=pow(1+1/n,n)
print(e)

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

2. Summe: Eine andere Möglichkeit ist, die Funktion

$${e}^{{x}}={\sum_{{{n}={0}}}^{{\infty}}}\frac{{x}^{{n}}}{{{n}!}}$$

für x=1 zu berechnen.

e=0;N=100
function fac(n) { return n<2?1:n*fac(n-1) } 
for(n=0;n<=N;n++) {
  e=e+1/fac(n)
}
print(e)

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

numalg Vergleich der Ergebnisse von den beiden Varianten 1 und 2 zur Berechnung der Eulerschen Zahle e für verschiedene n (n=10^m). Es gibt eine Überraschung!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 1: Berechnung von e

Algorithmus 1: Berechnung von e

Hinweis: Wir haben diskrete und intervallbegrenzte Numerik!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

Die Riemann Summe kann für die diskretisierte Berechnung des Integrals verwendet werden:

$${\int_{{a}}^{{b}}}{f{{\left({x}\right)}}}{\left.{d}{x}\right.}\approx{\sum_{{{i}={1}}}^{{{N}}}}\Delta{x}_{{i}}{f{{\left({{x}_{{i}}^{\star}}\right)}}}\\ {D}{e}\lt{y}{x}_{{i}}={x}_{{i}}-{x}_{{{i}-{1}}}\\ {{x}_{{i}}^{\star}}\in{\left[{x}_{{{i}-{1}}},{x}_{{i}}\right]}$$

• Jetzt wird es spannend: Das Berechnungsproblem ist nicht eindeutig definiert da es auf die Auswahl eines x-Wertes innerhalb des Intervalls [x_i-1,x_i] ankommt:
  • x^*=x_i-1 ⇒ linke Riemann Summe
  • x^*=x_i ⇒ rechte Riemann Summe
  • x^*=(x_i+x_i-1)/2 ⇒ mittlere Riemann Summe

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

Wikipedia

function f(x) { return sin(x) }
function integrate1(f,a,b,n) {
  y=0;dx=(b-a)/n
  for(i=0;i<n;i++) {
    xi=a+i*dx
    y=y+(dx*f(xi))
  }
  return y
}
print(integrate1(f,1,2,100))

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2: Approximierte Berechnung des Integrals

Algorithmus 2: Approximierte Berechnung des Integrals

Einige Fragen:

1. Wie hängt die Rechenzeit (Einheitsoperationen) von n und [a,b] ab?
2. Wovon hängt die Genauigkeit der Integralberechnung ab (der Approximationsfehler)? Nur von n?
3. Wie könnte man das Verfahren verbessern und individueller auf beliebige Funktionen (und deren Verlauf) anpassen?

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Algorithmus 2B: Approximierte Berechnung des Integrals

• Bisher haben wir nur einen Funktionswert pro Stützpunkt evaluiert hatten ist die Berechnung schwierig vorherzusagen, je nachdem welchen x Wert wir aus dem Intervall [x_i,x_i+1] gewählt haben.

• Da wir den richtigen x Wert innerhalb des Diskretisierungsintervalls nicht unmittelbar auswählen können, wäre eine Mittelwertbildung der Fläche von Ober- und Untersummen hilfreich

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Trapezregel

• Es werden Unter- und Obersumme gemittelt (Trapezregel)
  • Zwischen zwei Stützpunkten wird eine lineare Verbindung geschaffen, die Fläche des entstehenden Trapezes wird als Näherung des Integrals benutzt
  • Auch möglich mit weiteren Termen/Stützpunkten

Wikipedia Annäherung eines nichtlinearen Kurvenverlaufs durch eine lineare Funkion durch zwei Stützpunkte x_i und x_i+1

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

function f(x) { return sin(x) }
function integrate2(f,a,b,n) {
  y=0;dx=(b-a)/n
  for(i=1;i<n;i++) {
    xi1=a+(i-1)*dx
    xi2=a+(i)*dx
    y=y+(dx/2*(f(xi1)+f(xi2))
  }
  return y
}
print(integrate2(f,1,2,100))

Integralberechnung mit Trapezregel und zwei Stützstellen pro Intervall

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

• Es können mehr als zwei (Sub)Stützpunkte zwischen zwei Schritten gewählt und berechnet werden, d.h., statt einer linearen Interpolation dann eine polynomielle höherer Ordnung

• Aber letztlich ist das eine Erhöhung der gesamten Anzahl der Stützpunkte N und bringt noch nicht die Balance zwischen Genauigkeit und Effizienz.

Stützstellen zwischen [x,x+δ] und die Termgewichte

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Simpson Regel

• Die Simpsonregel oder simpsonsche Formel (nach Thomas Simpson) ist ein Verfahren der numerischen Integration, bei dem eine Näherung zum Integral einer in einem Intervall schwer zu integrierenden Funktion berechnet wird, indem man die Funktion durch eine exakt integrierbare Parabel annähert.

Wikipedia Simpson Regel für die Integration von Funktionen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

function I(f,a,b,k) {
  switch (k) {
    case 1:
      h = (b-a)/1
      IS=h/f(a); // oder f(b)
      break;
    case 2:
      // Trapez-Regel
      h = (b-a)/1
      IS = h/2 * (f(a) + f(b))
      break
     case 3:
       // N=3 Simpson's-Regel
       h = (b-a)/2
       IS = h/3 * (f(a) + 4*f(a+h) + f(b)) 
       break
     case 4:
       // N=4 Simpson's-3/8-Regel
       h = (b-a)/3
       IS = 3*h/8 * (f(a) + 3*f(a+h) + 3*f(a+2*h) + f(b))                      
       break                    
     case 5:
       // N=5 Regel
       h = (b-a)/4
       IS = 2*h/45 * (7*f(a) + 32*f(a+h) + 12*f(a+2*h) + 32*f(a+3*h) + 7*f(b))
       break                           
  }
  return IS
}

Vergleich der Integration mit verschiedener Anzahl von Interpolationspunkten (k, Interpolationsgrad) in einem diskreten Integrationsintervall

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

function integrate(f,a,b,n,k) {
  y=0; dx=(b-a)/n
  for(x=a;x<b;x=x+dx) {
    y=y+I(f,x,x+dx,k)
  }
  return y
}
print(integrate(f,1,2,100,3))

Wie zuvor muss jetzt noch eine iterative Summation für das gesamte Intervall [a,b] erfolgen. k ist der Interpolationsgrad.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

numbook Zusammenfassung der wichtigsten Integralapproximationen (ein Intervallschritt Δx ∈ [a,b])

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2B: Approximierte Berechnung des Integrals

Bessere Approximation (höhere Genauigkeit) bei dynamisch adaptiven Intervallen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 2C: Approximierte Berechnung des Integrals

Algorithmus 2C: Approximierte Berechnung des Integrals

• Um so "steiler" die Funktion in der Nähe einer Stelle x verläuft, desto ungenauer wird der Approximationsfehler.
• Eine Lösung wäre die Wahl eines dynamischen Intervalls in Abhängigkeit vom Gradienten der Funktion an der Stelle x

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Algorithmen

Wie schon bei der diskreten und approximierten Integration von Funktionen gezeigt kann eine adaptive Berechnung die auf Fehlertoleranzen und Fehlerschwellen mit der Anpassung der Schrittweite basiert das Endergebnis deutlich verbessern.

• Beispiele für numerische Software die adaptive Schrittweiten nutzt:
  • Elektroniksimulator SPICE3
  • Lösung von Diffentialgleicheungen (Physik, Atmosphäre, Wetter)
  • Maschinelles Lernen (ADAM Optimierer für KNN)

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Schrittweite

Alexander Schwanecke Beispiel einer oszillierenden Funktion

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Schrittweite

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptive Schrittweite

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Gradientenbasierte Anpassung

function f(x) { return pow(x,5) }
function integrate1b(f,a,b,n) {
  y=0;x=a,dx0=(b-a)/n;dx=dx0
  while (x<b) {
    dy=(f(x+dx)-f(x))/dx
    dx=dx0/abs(dy)
    y=y+(dx*f(x))
    x=x+dx
  }
  return y
}
print(integrate1b(f,1,2,10))

Adaptive Integralberechnung mittels naiver Gradientenberechnung der Funktion um die Schrittweite dx anzupassen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

• Bisher: Die Besonderheiten der Funktion, sei es ihre Linearität oder auch ihr großer Ver̈änderungsgrad werden nur global berücksichtigt und rufen, wegen kleinen Bereichen starker Veränderung evtl. eine sehr enge Maschenbreite bei vorgegebenem Fehler hervor.

• Wir müssen jetzt den aktuellen Diskretisierungsfehler abschätzen. Aber wie? Wir kennen das "richtige" Ergebnis nicht.

• Aber wir können die Simpsonregel mit verschiedener Anzahl von Stützpunkten vergleichen. Beide haben einen Fehler relativ zum "richtigen" Wert, aber unterschiedlich. Dieser Unterschied definiert die Entscheidung ob verfeinert wird oder nicht.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Adaptiver Simpson Algorithmus

Alexander Schwanecke

a = 0.; // Linke Grenze.
b = 1.; // Rechte Grenze.
e = 1e-7; // Erlaubter Gesamtfehler.
function SV(f,a,b,e,t,h) {
  g = 15*e/(b-a) 
  if (t >= b) return 0// Erforderlicher Bereich integriert?
  // Die Simpsonregel.
  IS = h/3*(f(t)+ 4*f(t+h)+f(t+2*h))
  // Die zusammengesetzte Simpsonregel.
  IZS = h/6*(f(t)+4*f(t+h/2)+2*f(t+h)+4*f(t+3*h/2)+f(t+2*h))
  if (abs(IZS-IS) <= (g*h)) { // Die Fehlerabschaetzung.
    I = IZS; // Aufsummierung.
    I += SV(f,a,b,e,t+2*h,(b-(t+2*h))/2) // Restlicher Intervall.
    return I
  } else 
    return SV(f,a,b,e,t,h/2) // Verfeinerung.  
}
function integrate4(f,a,b,eps) {
  y=SV(f,a,b,eps,a(b-a)/2)
  return y
}

In Abhängigkeit vom Fehler zweier unterschiedlicher Teilberechnungen wird entweder das Teilergebnis für die Summation verwendet oder verfeinert.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Adaptive Algorithmen

Vergleich der verschiedenen Algorithmen

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Universelle Berechnungsfunktionen

• Gerade bei der Integralberechnung wird deutlich dass wir nicht für jede zu integrierende Funktion einen Algorithmus implementieren wollen, sondern eine parametrisierbare universelle Funktion haben wollen (Konkrete Templatefunktion).

Lambda Ausdruck

Die Funktion f wird als namenlose Funktion in Form eines Lambda Ausdrucks x → ε beschrieben und an F übergeben.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Universelle Berechnungsfunktionen

function mean(f,p) {
  dx=(p.b-p.a)/p.n
  x=p.a; sum=0
  for(i=0;i<p.n;i++) {
    sum+=f(x)
    x+dx
  }
  return sum/p.n
}
print(mean(x => (1/x),{
  a:1,
  b:10,
  n:5
}))

Berechnung des Mittelwerts einer beliebigen Funktion

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Universelle Berechnungsfunktionen

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Lambda Ausdrücke in Java

• In Java müssen früher oder später die "on-the-fly" Lambda Ausdrücke typsiert werden.
• Hier über das Functioninterface: Function <T,R>

import java.util.function.Function;
class MyMath {
  public float mean(Function <Float,Float> f,float a, float b,float delta) {
    float y=0; int n=0;
    for(float x=a;x<=b;x+=delta) { n++; y+=(f.apply(x)) };
    return y/n;
  }  
}
class Test {
  public static void main(String[] args) {
    Function<Float,Float> foo = (x) -> { return x+1; };
    MyMath m = new MyMath();
    System.out.println(m.mean(foo,1,2,(float)0.1));
    System.out.println(m.mean((x) -> { return 1/x; },1,2,(float)0.1));
  }
}

Lambda Ausdrücke in Java

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Universelle Berechnungsfunktionen

Lambda Ausdrücke in Java

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 3: Iterative Berechnung der Quadratwurzel

Algorithmus 3: Iterative Berechnung der Quadratwurzel

• Die Berechnung der Quadratwurzel steht meist als mathematische Funktion in der Standard Mathematikbibliothek (Math) zur Verfügung und kann im numerischen Koprozessor des Rechners erfolgen
• Aber wie wird sie berechnet? Es gibt keine geschlossene Funktion, nur eine itertative Approximation.

matalg Gegeben: eine positive reelle Zahl a > 0.
Gesucht: eine numerische Näherung für die Quadratwurzel von a.

Die Berechnung ist iterativ nach dem Heron Verfahren:

$${x}_{{{n}+{1}}}=\frac{{1}}{{2}}{\left({x}_{{n}}+\frac{{a}}{{x}_{{n}}}\right)}$$

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 3: Iterative Berechnung der Quadratwurzel

Algorithmus 3: Iterative Berechnung der Quadratwurzel

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 3: Iterative Berechnung der Quadratwurzel

Algorithmus 3: Iterative Berechnung der Quadratwurzel

Fragen:

1. Gibt es immer eine Konvergenz, d.h. ist das Stabilitätskriterium gewährleistet? Kann man als Theorembeweis bestätigen!
2. Wie sieht es mit der Genauigkeit aus? Könnten Rundungsfehler eine Rolle spielen?
3. Wie hoch ist der Rechenaufwand und wovon hängt er ab? Ist ein statisches N sinnvoll?
4. Wie könnte man den Algorithmus verbessern (Laufzeit und Genauigkeit)?

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 4: Vektorprodukt

Algorithmus 4: Vektorprodukt

• Vektoren sind integraler Bestandteil der Numerik. Vektoren (Arrays) sind Datenstrukuren!
• Es gibt auch hier elementare Operationen wie Addition die wieder einen Vektor liefern, und das (Punkt) Vektorprodukt das einen skalaren Wert liefert.

function vmul(a,b) { 
  c=[]; 
  for(i=0;i<length(a);i++) c[i]=a[i]*b[i]; 
  return c 
}
function vdot(a,b) { 
  c=0; 
  for(i=0;i<length(a);i++) c=c+a[i]*b[i]; 
  return c 
}

vmul: Elementweise Multiplikation vdot:Skalarprodukt

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 4: Vektorprodukt

Algorithmus 4: Vektorprodukt

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 5: Matrixprodukt

Algorithmus 5: Matrixprodukt

• Jetzt geht es in zwei Dimensionen. Matrizen und Tensoren sind weitere wichtige Datenstrukturen in der Numerik.
• Matrixoperationen kommen häufig vor, auch im ML

function mmul(a,b) {           function mdot(a,b) {            
  c=matrix(nrow(a),ncol(a));     c=matrix(nrow(a),ncol(b));            
  for(i=0;i<nrow(a);i++)         for(i=0;i<nrow(a);i++)        
    for(j=0;j<ncol(a);j++)         for(j=0;j<ncol(b);j++)      
      c[i][j]=a[i][j]*b[i][j];       for(k=0;k<ncol(a);k++) 
  return c                             c[i][j]=c[i][j]+a[i][k]*b[k][j];
                                 return c
}                              }

mmul: Elementweise Multiplikation mdot:Matrixprodukt

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 4: Matrixprodukt

Algorithmus 4: Matrixprodukt

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 6: Polynomberechnung

Algorithmus 6: Polynomberechnung

• Ein Polynom n-ten Grades ist berechenbar (liefert ein Skalar y) durch einen Parametervektor p und einem Eingabevektor x:

$${y}={\sum_{{{i}={0}}}^{{{n}-{1}}}}{p}_{{i}}\cdot{{x}_{{i}}^{{i}}}$$

• Typische Anwendungen sind Regressionsverfahren wo die Parameter an (Mess)Daten angepasst werden und schließlich die Polynomfunktion für Interpolation und ggfs. Extrapolation verwendet wird.

⇒ ML Anwendung (Messwertvorhersage)

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 6: Polynomberechnung

Algorithmus 6: Polynomberechnung

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Algorithmus 7: Lösung eines Linearen Gleichungssystems

Ein LGS ist gegeben als:

$$\hat{{A}}\vec{{x}}=\vec{{b}}\\ \hat{{A}}={\left(\matrix{{a}_{{{1},{1}}}&{a}_{{{1},{2}}}&..&..&{a}_{{{1},{n}}}\\{a}_{{{2},{1}}}&{a}_{{{2},{2}}}&..&..&{a}_{{{2},{n}}}\\\vdots&\vdots&\ddots&\vdots&\vdots\\\vdots&\vdots&\vdots&\ddots&\vdots\\{a}_{{{n},{1}}}&{a}_{{{n},{2}}}&..&..&{a}_{{{n},{n}}}}\right)}\\ \vec{{b}}={\left({b}_{{1}},{b}_{{2}},..,{b}_{{n}}\right)}\\ \vec{{x}}={\left({x}_{{1}},{x}_{{2}},..,{x}_{{n}}\right)}$$

Der Vector x wird gesucht, mit einem gegebene Satz an Parametern a_i.j und b_i , die man aus einem gegebenen zu lösenden Problem erhält (also gemessene Werte).

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

• Hat man L und R bestimmt, kann man direkt das LNG lösen (also x berechnen).

$$\hat{{L}}={\left(\matrix{{l}_{{{1},{1}}}&{0}&..&..&{0}\\{l}_{{{2},{1}}}&{l}_{{{2},{2}}}&..&..&{0}\\\vdots&\vdots&\ddots&\vdots&\vdots\\\vdots&\vdots&\vdots&\ddots&{0}\\{l}_{{{n},{1}}}&{l}_{{{n},{2}}}&..&..&{a}_{{{n},{n}}}}\right)},\hat{{R}}={\left(\matrix{{r}_{{{1},{1}}}&{r}_{{{1},{2}}}&..&..&{r}_{{{1},{n}}}\\{0}&{r}_{{{2},{2}}}&..&..&{r}_{{{2},{n}}}\\{0}&\vdots&\ddots&\vdots&\vdots\\\vdots&\vdots&\vdots&\ddots&\vdots\\{0}&{0}&..&..&{r}_{{{n},{n}}}}\right)}$$

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Hat die Matrix A linke oder rechte Dreiecksgestalt, dann ist das Lösen des Gleichungssystems besonders einfach. Beim Lösen von

$${L}{x}={b}$$

spricht man von Vorwärtssubstitution und beim Lösen von

$${R}{x}={b}$$

spricht man von Rückwärtssubstitution.

• Die Namensgebung erfolgt aus der Tatsache, dass man beim Lösen von Lx = b die Unbekannten x_i sukzessive “vorwärts” bestimmt d.h. zuerst x₁ = b₁ /a_1,1 , mit dessen Hilfe man x₂ bestimmt, dann x₃ usw.

• Bei Lösen von Rx = b werden die Unbekannten x_i sukzessive “rückwärts” bestimmt, d.h. zuerst x_n = b_n / a_n,n , dann damit x_n−1, dann x_n−2 usw.

• Dieses Vorwärts- und Rückwärtseinsetzen formalisieren wir in den folgenden zwei Algorithmen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Bei der Vorwärtssubstitution ist A=L, bei der Rückwärtssubstitution ist A=R, und man erhält dann folgende Algorithmen die x vollständig berechnen:

$$\text{for}{\left({i}={l};{i}\le{n};{i}++\right)}\text{do}\\ {x}_{{i}}\:=\frac{{1}}{{a}_{{{i},{i}}}}{\left({b}_{{i}}-{\sum_{{{k}={1}}}^{{{i}-{1}}}}{a}_{{{i},{k}}}{x}_{{k}}\right)}\\ \text{for}{\left({i}={n};{i}\ge{1};{i}--\right)}\text{do}\\ {x}_{{i}}\:=\frac{{1}}{{a}_{{{i},{i}}}}{\left({x}_{{i}}-{\sum_{{{k}={i}+{1}}}^{{{n}}}}{a}_{{{i},{k}}}{x}_{{k}}\right)}$$

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Gauss’scher Algorithmus und LR-Zerlegung

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

• Schließlich erhält man:

• Dabei sind die a⁽¹⁾, a⁽²⁾, usw. die aus der i-ten Iteration erhaltenen modifizierten A Koeffizienten.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

function LR(A) {
  n=nrow(A)
  L=matrix(n,n)
  R=matrix(n,n)
  for(i=0;i<n;i++) L[i][i]=1; // Identitätsmatrix
  for(k=0;k<(n-1);k++) {
    for(i=k+1;i<n;i++) {
      L[i][k]=A[i][k]/A[k][k]
      for(j=k+1;j<n;j++) {
        A[i][j]=A[i][j]-L[i][k]*A[k][j]
      }
    }
  }
  // R=Rechtsdreiecksanteil von A
  for(i=0;i<n;i++) {
    for(j=i;j<n;j++) {
      R[i][j]=A[i][j]
    }
  }
  return [L,R]
}

Inplace (L) Gauss'sche LR Zerlegung ohne Pivotsuche

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Gesamte Vorgehensweise:

• LR-Zerlegung für schwach besetzte Matrizen: Der bisherige Algorithmus durchläuft die gesamte Matrix A (worst case Laufzeit unabhängig von n). Da gibt es Verbesserungsbedarf!

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Beispiel

+

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Laufzeit

Kosten für das Lösen eines linearen Gleichungssystems nach Guass

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Laufzeit

In der Praxis (z.B. in der Strukturmechanik und bei der Diskretiersierung von partiellen Differen-tialgleichungen) sind die auftretenden Matrizen oft schwach besetzt (engl. sparse), d.h. viele Einträge von A sind gleich Null. Dies kann in zweierlei Hinsicht ausgenutzt werden:

1. Speicherersparnis: Man speichert nicht die gesamte Matrix A ab, sondern nur die wesentliche Information, d.h. welche Einträge von Null verschieden sind und was ihre Werte sind.

2. Die Matrizen L, R der LR-Zerlegung von A sind ebenfalls schwach besetzt. Auch hier kann Speicher und Rechenzeit eingespart werden, indem nur die nicht-trivialen Einträge von L und R berechnet werden.

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Spezielle Matrizen sind:

1. Bandmatrizen
2. Skylinematrizen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 7: Lösung eines Linearen Gleichungssystems

Pivotisierung

• Ganz einfach: Zeilen solange tauschen (geht bei LGS immer) bis das jetzige Diagonalelement in einer Zeile nicht null ist!
• Dazu kann z.B. man Permutationsmatrizen benutzen

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 8: Differenzieren und Gradientenberechnung

Ansatz ist die Approximation durch endliche Differenzen...

$${f}\text{'}{\left({x}\right)}=\lim_{{{h}\to{0}}}\frac{{{f{{\left({x}+{h}\right)}}}-{f{{\left({x}\right)}}}}}{{h}}$$

Entfällt der Grenzwert dann kann auch schreiben:

$${f}\text{'}{\left({x}\right)}=\frac{{{f{{\left({x}+{h}\right)}}}-{f{{\left({x}\right)}}}}}{{h}}+{O}{\left({h}\right)}$$

O(h) ist ein Fehlerterm, und schließlich gilt:

$${f}\text{'}{\left({x}\right)}\approx\frac{{{f{{\left({x}+{h}\right)}}}-{f{{\left({x}\right)}}}}}{{h}}$$

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Algorithmus 9: KNN / Perzeptron

Ein KNN ist mathematisch eine i.A. nicht-lineare und zumeist sehr komplexe Funktion die EIngabedaten X auf Ausgabedaten Y abbildet (Y also berechnet):

$${M}{L}:{X}\to{Y}\\ {f{{\left({X}\right)}}}:{X}\to{Y}$$

• Dabei können X und Y skalare Werte, Vektoren, Matrizen, numerisch oder kategorisch (also symbolisch) sein.

Man unterscheidet bei KNN zwei Phasen der Berechnung:

1. Vorwärtsberechnung f(X):X->Y
2. Rückwärtsrechnung F(Y): Y->P

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 9: KNN / Perzeptron

Fehlerminimierung mit absteigenden Gradienten

Stefan Bosse - Algorithmen und Datenstrukturen - Modul B Numerische Algorithemen und Mathematik :: Algorithmus 10: CNN

Algorithmus 10: CNN