%% Bonusblatt 11 - Programmieraufgaben:
%% Dünngitterquadratur (8 Punkte)
% Vorname Nachname Mat.-Nr.
% In der Vorlesung haben wir die Trapezregel als einfachste Quadraturformel kennengelernt. Wir betrachten dazu das Intervall I = [0, 1] und den Fall N = 2^j für ein j \in \mathbb{N}. In diesem Falle ist h = 1/N = 2^{-j}, und wir definieren x_i = ih, i = 0, 1, \ldots, 2^j. Wir können dann für eine stetige Funktion f das Integral
% \int_0^1 f(x) \ {\rm d} x
% mithilfe der zusammengesetzen Trapezregel annähern, welche auf dem Level j durch
% T_j[f] : \! = \frac{h}{2} f(0) + h \sum_{i = 1}^{2^j-1} f(x_i) + \frac{h}{2} f(1)
% definiert ist.
% Dieses Vorgehen lässt sich leicht auf zwei, drei und im Allgemeinen beliebig viele Dimensionen übertragen. Ist etwa f \colon \mathbb{R}^2 \to \mathbb{R}, möchten wir eine Annäherung des Integrals
% \int_0^1 \int_0^1 f(x, y) \ {\rm d}x \ {\rm d}y
% berechnen. Der Einfachkeit halber wollen wir uns dabei auf Funktionen beschränken, welche am Rand verschwinden, das heisst, wir nehmen f(0, y) = f(1, y) = f(x, 0) = f(x, 1) = 0 für alle x, y \in [0, 1] an. Die zweidimensionale Trapezregel können wir dann mittels
% (T_j \otimes T_j)[f] : \! = h^2 \sum_{i, k = 1}^{2^j-1} f(ih, kh)
% definieren.
% Offensichtlich benötigt man für die zweidimensionale Trapezregel \mathcal{O}(N^2) = \mathcal{O}(4^j) Funktionsauswertungen. Man kann zudem zeigen, dass die zweidimensionale Trapezregel eine hinreichend glatte Funktion f bis auf \mathcal{O}(h^2) = \mathcal{O}(4^{-j}) genau zu integrieren vermag. Während der Fehler also im Vergleich zum eindimensionalen Fall asymptotisch gleich bleibt, wächst die Komplexität von linear in einer Dimension auf quadratisch in zwei Dimensionen und allgemein auf \mathcal{O}(2^{dj}) in d Dimensionen. Dieses Wachstum ist in höheren Dimensionen sehr schnell, womit diese Routine im Allgemeinen nicht gut gebraucht werden kann.
% Eine Alternative bietet die Dünngitterquadratur, welche durch eine geschickte Wahl der Quadraturpunkte und -gewichte die Komplexität auf \mathcal{O}(j^{d-1}2^j) in d Dimensionen reduziert. Dabei nehmen wir eine leichte Einbusse bei der Genauigkeit in Kauf, welche nun \mathcal{O}(4^{-j}j^{d-1}) beträgt. Da alle Potenzen von j gegenüber jedem exponentiellen Wachstum und exponentiellen Abklingen asymptotisch vernachlässigbar sind, ist das kein grosses Hindernis. Zusammengefasst haben wir somit also eine Quadraturroutine, welche einen essenziell linearen Aufwand beim einem essentiell bestmöglichen Fehler hat.
% Um die Dünngitterquadratur zu berechnen, gehen wir wie folgt vor:
% Input: Funktion f \in C([0, 1]^2), Level j
% Output: Approximation Q \approx \int \int f
%
% setze Q : \! = 0
% for \ell_x = 0, 1, \ldots, j-1
%   setze h_x : \! = 2^{-(\ell_x + 1)}
%     for \ell_y = 0, 1, \ldots, j - \ell_x - 1
%       setze h_y : \! = 2^{-(\ell_y + 1)
%       for k_x = 0, 1, \ldots 2^{\ell_x} - 1
%         for k_y = 0, 1, \ldots, 2^{\ell_y} - 1 datiere auf
%    Q = Q + h_x h_y \mathcal{S}_{h_x, h_y} \big( h_x(2 k_x + 1), \ h_y (2 k_y + 1) \big) [f]
% Dabei bezeichne
% \mathcal{S}_{h_x, h_y} (x, y)[f] : \! = f(x, y) - \frac{1}{2} \big( f(x + h_x, \, y) + f(x - h_x, \, y) + f(x, \, y + h_y) + f(x, \, y - h_y) \big) \\ \qquad \qquad \qquad \qquad + \frac{1}{4} \big( f(x + h_x, \, y + h_y) + f(x + h_x, y - h_y) + f(x - h_x, y + h_y) + f(x - h_x, y - h_y) \big)
%
% den hierarchischen Überschuss an der Stelle (x, y).
%
%
% Schreiben Sie hier am Ende dieses Live Scripts MATLAB-Funktionen
% Q = trapezoidalRule2D(f, j),
% Q = sparseGridQuadrature(f, j),
% welche die Trapezregel in zwei Dimensionen sowie die Dünngitterquadratur implementieren. Testen Sie Ihre Implementierung an der Funktion f(x, y) = \sin(7 \pi x) \sin(13 \pi y) mit der exakten Lösung
% \int_0^1 \int_0^1 f(x, y) \ {\rm d}x \ {\rm d}y = \frac{4}{91 \pi^2}.
% Plotten Sie für j = 5, 6, \ldots, 15 die Fehler der Trapezregel sowie der Dünngitterquadratur in einem semilogy-Plot gegen j. Plotten Sie zusätzlich Vergleichsgeraden C_14^{-j} und C_2 j 4^{-j}, wobei C_1 und C_2 so gewählt sein sollen, dass klar ersichtlich ist, dass die Fehlerraten parallel zu den Vergleichsgeraden verlaufen.
% Vergleichen Sie in einem zweiten semilogy-Plot die Laufzeiten der Trapezregel und der Dünngitterquadratur gegen j. Plotten Sie hier auch die Vergleichsgeraden C_34^{j} und C_4 j 2^j für sinnvoll gewählte C_3 und C_4.
% Versehen Sie Ihre Plots mit einer Legende und achten Sie darauf, ihre Linien farblich zu kennzeichnen, das heisst, es muss klar ersichtlich sein, welche Linie wozu gehört. Falls Sie Ihr Livescript schwarzweiss ausdrucken, dürfen Sie gerne die Linien nachträglich farbig einzeichnen.
%

% TODO: implementieren
% f = ...

figure()
subplot(1, 2, 1)
% TODO: Fehler plotten
% semilogy(...
% legend(...
title('Error')

subplot(1, 2, 2)
% TODO: Laufzeiten plotten
% semilogy(...
% legend(...
title('Runtime')

%
% Funktionen:

function Q = trapezoidalRule2D(f, j)
    % TODO: implementieren
end

function Q = sparseGridQuadrature(f, j)
    % TODO: implementieren
end