Computational Mathematics

Allgemeine Hinweise

• Bitte modifizieren Sie nicht die vorgegebenen Methodenköpfe. Beachten Sie auch, ob Eingabevektoren zeilen- oder spaltenweise angegeben werden müssen (Sofern in der Aufgabenstellung vermerkt).
• Fügen Sie die Datei gruppe.txt bitte in die Ordner serie[Nummer]. Das erleichtert uns die Zuordnung und Ihnen ggf. den Wechsel der Gruppe während des Semesters.
• Sie erschweren die Kontrolle, durch ein Umbenennen der Verzeichnisse (SerieX, serie_X, Serie X,...). Bitte halten Sie sich an die Syntax P-NUM/serie[Nummer]/[Methodenname].m
• Falls Sie zu zweit arbeiten, splitten Sie die Abgabe bitte nicht auf beide Accounts.

Fehler in den Abgaben von Aufgabe 2.4

Die Einträge des Neville-Schema müssen so berechnet werden, dass auf die vorigen Werte zurück gegriffen werden kann. Es wurden Punkte abgezogen, falls die Lösung mittels Rekursion berechnet wurde, die zur mehrfachen Berechnung der Zwischenpunkte führte.
Aus Stabilitätsgründen wurde in der Vorlesung auf eine kleine Umformulierung des Neville-Schemas hingewiesen. Es wurden aber keine Punkte abgezogen, falls diese nicht verwendet wurde.

Fehler in den Abgaben von Aufgabe 3.3

• Der Fall, dass x einer Stützstelle entspricht, wurde nicht abgefangen oder zumindest im Funktionsheader als nicht zulässig erwähnt (-0.5P)
• Berechnung von ω'(x) mittels Symbolischem Differenzieren. Ziel der Aufgabe war es, einige Rechenoperationen durch Umformulierung das Problems zu sparen. Der vergleichsweise hohe Aufwand des Symbolischen Differenzieren hat genau den gegenteiligen Effekt (-1P).
• Der Aufruf der Methode ist nicht erfolgreich, wenn x∈ℝ^n vektorwertig ist (-1P). Die Berechnung der A_j erzeugt einen Overhead, der nur dann gerechtfertigt ist, wenn man die Koeffizienten für mehrere Auswertungen von p(x) verwenden kann.
• Falsche Indizierung von Schleifen und Vektoren. Einige Methoden setzten length(x)=length(stuetz) voraus.

Fehler in den Abgaben von Aufgabe 3.4

• Berechnung der dividierten Differenzen mittels Rekursion und ohne Speicherung der Zwischenwerte (vgl. Aufgabe 2.4, -1P).
• Mit Hilfe der Operatoren .* und ./ ist die Erweiterung auf vektorwertige Operationen meist sehr einfach. Achtung, falls ein Vektor mit einer reellen (komplexen) Zahl addiert wird, erfolgt die Addition elementweise.
  Beispiel:
  

  %Berechnung des Polynoms mit Hilfe des Horner-Schemas
  p=f(n);
  for i=1:(n-1)
    p= p.*(x-stuetz(n-i))+f(n-i);
  end

• Die Polynome wurden nicht bezüglich mit der Newtowbasis dargestellt (-2P).
• Berechnung des Polynoms ohne Horner-Schema (-2P).
• (Aus dem Code einer Abgabe)

  % Bei jedem Fkt.aufruf muessen wir leider STUETZ und FWERTE
  % uebergeben (globale Vars bei MatLab...?).

  Matlabhilfe zu globalen Variablen: help global. In den meisten Fällen ist aber von einer Verwendung globaler Variablen abzuraten. Insbesondere beim Schreiben von Bibliotheksfunktionen.
• (Aus dem Code einer Abgabe)

  % Außerdem sollte die Fkt. B natuerlich nicht von außerhalb
  % aufrufbar sein (private functions bei MatLab...?).

  Matlabhilfe: help function.
  Beispiel:
  

  % Datei A.m
  function A()
  
      function priv()
      % Verschachtelte Methoden-Deklarationen sollten vermieden werden.
         disp('Private Funktion von A.'); 
      end
  
  priv(); % Nur innerhalb von A aufrufbar.
  B();  % Nur innerhalb von Methoden aus dieser Datei aufrufbar.
  end
  
  function B()
  C();
  end
  
  function C()
  disp('Nur innerhalb von A.m sichtbar.');
  end
  

Fehler in den Abgaben zu Aufgabe 4.4

• Der zentrale Differenzenquotient besitzt die Ordnung q=2. Für die Genauigkeit der Extrapolation ist es wichtig dies auch zu nutzen und das Interpolationspolynom „bezüglich h^q“ zu bestimmen (siehe Satz 1.12 und den folgenden Abschnitt der Vorlesung). Viele haben das Neville-Schema der zweiten Serie nicht modifiziert (-1P).

Fehler in den Abgaben zu Aufgabe 5.4

• Keine Skalierung der Čebyšev-Knoten in das gegebene Intervall.
• Mögliche Antwort bei der Frage, welche Stützstellen besser geeignet sind: Der Teilausdruck π_j=0,...,n(x-x_j) aus der Fehlerabschätzung der Lagrange-Interpolation auf [-1,1] wird bei Wahl der Čebyšev-Knoten bezüglich der ∞-Norm minimiert.

Kommentare zu den Abgaben von Aufgabe 6.4

• Einige Studenten haben nur die Matlabroutine spline(...), etc., aufgerufen. Davon abgesehen, dass es dafür keine Punkte gab, berechnet die Methode auch nicht den natürlichen Spline.
• Die Methode sollte die berechneten Werte (yy) auch zurückgeben. Haben Sie immer im Hinterkopf, dass Ihre Methode in anderen Programmen als „Bibliotheksfunktion“ verwendet wird. Daher sollten Daten auch nicht innerhalb der Methode geplottet werden.
• Für schwach besetzte Matrizen gibt es in Matlab Datenstrukturen (doc sparse, doc spdiags), die nur die Nichtnulleinträge einer Matrix speichern und bei Rechenoperationen, wie z.B. der Matrix-Vektor-Multiplikation, keine quadratische Laufzeit erfordern. Sparse-Matrizen können in Matlab wie normale Matrizen behandelt werden.
• Versuchen Sie bei for-Schleifen (und if-Statements) eine hohe Verschachtelung zu vermeiden. Zum einen erschwert die Verschachtelung die Lesbarkeit des Codes. Zum anderen werden for-Schleifen von Matlab wesentlich schlechter parallelisiert, als äquivalente Befehle in „Matrix-Vektor-Schreibweise“.
  (Natürlich müssen Sie nicht alle Operationen vektorisieren.)

Fehler in den Abgaben von Aufgabe 7.4

• In vielen Abgaben wurden die (konstanten) Koeffizenten der L2-Projektion für jede Stützstelle (Schleife über xx) neu berechnet.
• Häufig wurde die Normierung der Polynome vergessen.
• Die Berechnung der Basisfunktionen mittels verschachtelter Functionhandles ist schlecht! (Warum?) Vermeiden Sie Konstruktionen der Form

  % Berechne Legendre-Polynome
      function y = legendrep(x,l)
          if (l==0)
              y=1;
          else if (l==1)
                  y=x;
              else
                  y=x.*legendrep(x,l-1)-((l-1)^2/(4*(l-1)^2-1)).*legendrep(x,l-2);
              end
          end      
      end
  

  Besser wäre

  % Berechne Legendre-Polynome
  function y=Legendre(x,k)
      if( k==0)
          y=ones(size(x));
          return
      elseif( k==1 )
          y=[ones(size(x,2));x];    
      else
          y = zeros(k+1,size(x,2));
          y(1:2,:)=[ones(1,size(x,2));x];    
              
           for i=1:k-1
               y(i+2,:) = x.*y(i+1,:) - (i)*(i)/(4*(i)*(i)-1) *y(i,:);            
           end
      end
      % Nur das k-te Polynom zurückgeben
      y=y(end,:);  
  end
  

• In einer Abgabe wurden die Legendre-Polynome per symbolischer Differentiation berechnet. Das ist extrem langsam und unbedingt zu vermeiden. Man kann diesen Ansatz aber nutzen, um die Koeffizienten der Legende-Polynome einmalig(!) zu berechnen und sie bei Programmstart zu laden. Die Speicherung der Koeffizenten hat allerdings keinen nennenswerten Vorteil gegenüber der rekursiven Berechnung der Basispolynome (Die rekursive Berechnungs-Vorschrift der Basispolynome ist dem Horner-Schema sehr ähnlich).

Kommentare zu den Abgaben von Aufgabe 8.4

• Im Konvergenzplot sollten die Quadraturfehler, d.h. |Quad(f,a,b,N) - ∫_a^b f(x) dx|, berechnet werden, um zu prüfen, ob Ihre Methoden die, theoretisch zu erwartende, Konvergenzordnung besitzen. Es sollten nicht die theoretisch erwarteten Fehlerschranken geplottet werden.
• Kommentar aus einer Abgabe

  % versucht man gemaess aufgabenstellung das integral ueber [0,2pi] zu 
  % berechnen, welches - loesbar und - gleich null ist, so erhaelt man
  % das bild fehler2. Dieses ergibt sich aus rundungsfehlern beim berechnen
  % der stuetzstellen und der entsprechenden auswertung dort, erkennbar an 
  % der Groessenordnung der Fehler, ~10^(-14). Daher kann man
  % so keine Aussagen ueber das konvergenzverhalten der jeweiligen methoden
  % machen.
  % Wir haben uns deshalb entschieden selbiges Integral ueber [0,pi]
  % (mit dem Wert 2) zu approximieren. Dabei ergab sich das bild fehler,
  % auswelchem Konvergenzraten von ~N^(-2) fuer das Trapezverfahren und ~N^(-4)
  % fuer das Simpsonverfahren. Auch hier erkennt man, dass bei kleinen
  % fehlern ind der Region von 10^(-14) die Genauigkeit wieder abnimmt wie
  % oben beschrieben.
  

  Ja, im Falle der Integration bis 2⋅π löscht sich der Quadraturfehler aufgrund der Symmetrie aus. Allerdings hat der Abbruch der Konvergenz (bei Integration bis π einen anderen Grund. Die Genauigkeit nimmt ab, weil die Addition zweier Zahlen mit (im mathematischen Sinne) relativ großen Abstand zu Rundungsfehlern führt. In Ihrem Programm ist z. B. die Stelle

  Summe1=Summe1+f(stuetz(i));
  

  kritisch für N>>1. Man könnte dieses Problem durch Umordnen, Speicherung von Zwischenergebnissen oder die Verwendung eines anderen Datentypes umgehen. Allerdings ist es wesentlich effektiver, eine bessere Quadraturregel (siehe Serie 9) zu verwenden.
• Konstrukte der Form

  x(1)=0;
  for i=1:N-1
  	x(i+1)=x(i)+h;	
  end

  lassen sich leichter (und schneller) durch

   x=xstart:(xend-xstart)/N:xend  % oder
  linspace(xstart,xend,N);
  

  definieren.
• Auch die For-Schleifen beim Aufsummieren kann man so ersetzten. Beispielsweise (aus einem abgegebenen Code) ist

  h=(b-a)/N;
  y=0;
  for i=0:N
  	y=y+2*feval(f,a+h*i);
  end
  

  äquivalent zu

  
  y = 2* sum( f( a:h:b ) )
  

Kommentare zu den Abgaben von Aufgabe 9.4

• Ein wichtiger Faktor für die Effizienz der Romberg-Integration ist die Wiederverwertung alter Werte, um die Anzahl der Funktionsauswertungen zu minimieren.
• In einigen Abgabe wurde statt dem Abbruchkriterium vor der Romberg-Integration ein N=N(tol) definiert und N Schritte des Integrationsverfahrens gestartet. Damit kann nicht garantiert werden, dass die gewünschte Toleranz (effizient) erreicht wird.
• Bei korrekter Implementierung werden vier Schritte der Romberg-Integration durchgeführt.

  0.785399	0	0	0
  0.948059	1.002280	0	0
  0.987116	1.000135	0.999991565	0
  0.996785	1.000008	0.999999876	1.000000008144021

Hinweise zu Aufgabe 9.4

Wählen Sie ε hinreichend nahe 1, damit die Anzahl der Iterationsschritte N nicht n übersteigt. In Ihrer Auswertungsollte erkennbar sein, dass N für wachsendes n (und die gewählten Matrizen/ Startwerte) eine starke Überschätzung der notwendigen Anzahl von Schritten darstellt.

Hinweise zu Aufgabe 10.4

Für die Erstellung der Graphen bietet sich die Erweiterung des Interfaces auf function [x,Nend] = mynewton(f,Df,x0,N,tol) an, um außerhalb der Methode die Anzahl der benötigten Iterationschritte zugreifen zu können.

Kommentare zu den Abgaben von Aufgabe 11.4

• Da MATLAB komplexe Zahlen unterstützt, konnte auf eine Umformulierung in ℝ^2 verzichtet werden.
• Eine gute Idee war die Vektorisierung des Startpunktes x0, um das Newtonverfahren parallel für alle Werte eines Gitters zu berechnen.
• Die Methode zum Plotten eines Bildes sollte in der Regel erst nach Berechnung aller Werte (einmalig!) aufgerufen werden, da der Aufruf zum Zeichnen in Relation zum restlichen Programm sehr lange dauert.

Einige Bilder, die mit Hilfe Ihrer Skripte erzeugt wurden: Newton-Fraktal

Hinweise zu Aufgabe 12.4

Falls Sie eine GUI zum Einstellen des Anzeigebereiches schreiben möchten, können Sie mit guide ein Matlab-Tool starten, dass die Erstellung von Graphischen Interfaces erleichtert.

Einige Bilder, die mit Hilfe Ihrer Skripte erzeugt wurden: Mandelbrotmenge

Musterlösungen

Serie 3

Serie 6

Serie 7

Serie 8

Serie 11