<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Interaktives Übungsblatt
Zusammenfassung der Mathematik in der Übung
Matrixnorms
Spektrale Norm
Ist definiert als induzierte Vektornorm
$$\| A \| 
Step1:  Plotten Sie mithilfe von matrix_plot die Systemmatrizen für $\sigma = 0$ und $n=10$.
Step2:  Aufgabe
Step3:  Frage
Step4:  Frage 
Step5:  Frage
Step6:  Nur die hohen Frequenzen
Step7:  Frage
Step8:  Aufgabe
Step9:  Zweigitter-Iterationsmatrix
Step10:  Im Folgenden werden wir nun mithilfe des pymg frameworks die Zweigitter-Iterationsmatrix für ein einfaches Multigrid 
aufstellen. Wir beginnen mit der Grobgitterkorrektur.
Step11:  Buntebilderaufgabe
Step12:  Aufgabe
Step13:  Buntebilderaufgabe
Step14:  Aufgabe
Step15:  Aufgabe
Step16:  Asymptotische Äquivalenz  zwischen periodisch und nicht-periodisch
Wir sehen, dass die Spektralradiien auf den ersten Blick gut übereinstimmen. Wir wollen nun empirisch ergründen ob die Matrizenklassen der periodischen und nicht periodischen Fällen möglicherweise zueinander asymptotisch äquivalent sind.
Aufgabe
Step17:  Aufgabe
Step18:  Glättung
Step19:  Gauss-Seidel
Step20:  Grobgitterkorrektur
Hier trifft man mal wieder auf das Problem, dass die Freiheitsgrade im periodischen und nicht periodischen Fall unterschiedlich verteilt sind. Für einen Vergleich wird die Matrix um eine Nullzeile und Nullspalte erweitert.
Step21:  Zweigitter
<ASSISTANT_TASK:>
Python Code:
import sys
sys.path.append("/home/moser/MG_2016/pyMG/")
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pymg
from project.helmholtz1d import Helmholtz1D
from project.helmholtz1d_periodic import Helmholtz1D_Periodic
from project.gauss_seidel import GaussSeidel
from project.weighted_jacobi import WeightedJacobi
from project.pfasst.plot_tools import eigvalue_plot_list, matrix_plot, matrix_row_plot
from project.pfasst.transfer_tools import to_dense
from project.pfasst.matrix_method_tools import matrix_power
def system_matrix_hh1d(n,sig):
    hh1d = Helmholtz1D(n, sig)
    return hh1d.A
def system_matrix_hh1d_periodic(n,sig):
    hh1d = Helmholtz1D_Periodic(n, sig)
    return hh1d.A
def spec_rad(A):
    return np.max(np.abs(sp.linalg.eigvals(to_dense(A))))
Explanation: Interaktives Übungsblatt
Zusammenfassung der Mathematik in der Übung
Matrixnorms
Spektrale Norm
Ist definiert als induzierte Vektornorm
$$\| A \| := \max_{x \neq 0}\frac{\| Ax \|2}{\| x \|_2} = \max{\| x \|2 = 1} \| Ax \|_2.$$
Kann aber auch als maximaler Singulärwert 
$$\| A \|_2 = \sqrt{\mu\max}$$
dargestellt werden.
Hilbert-Schmidt Norm
Wir definieren die Hilbert-Schmidt Norm einer Matrx $A \in K^{n \times n}$ als
$$ |A| = \left( \frac{1}{n}\sum_{i = 0}^{n-1}\sum_{i = 0}^{n-1} |a_{i,j}|^2 \right)^{1/2}.$$
Es gilt
1. $|A| = \left( \frac{1}{n}\mbox{Spur}(A^A) \right)^{1/2}$
1. $|A| = \left( \frac{1}{n}\sum_{k=0}^{n-1}\lambda_k\right)^{1/2}$, wobei $\lambda_k$ die Eigenwerte von $A^A$ sind
1. $|A| \leq \|A\|$
Rayleigh-quotient
Sei $H$ eine hermitsche Matrix, so definiert sich der zugehörige Rayleigh quotient als
$$R_H(x) = \frac{x^Hx}{x^\cdot x}.$$
Der Rayleigh-quotient ist vor allem nützlich fü̇r Eigenwertabschätzungen.
So lässt sich der größte und kleinste Eigenwert $\lambda_M,\lambda_m$ abschätzen durch
$$ \lambda_m = \min_x R_H(x) \; \mbox{und} \; \lambda_M = \max_x R_H(x).$$
Funktion der Wiener Klasse
Sei ${t_k}{-\infty}^{\infty}$ eine absolut summierbare Folge, d.h. $\sum{k=-\infty}^{\infty}|t_k| < \infty$.
Sei desweiteren $f(\lambda) = \lim_{n \to \infty} \sum_{k=-n}^{n}t_k e^{-ik\lambda}$. So ist $f(\lambda)$ stetig und Riemann-Integrierbar.
Zyklische Matrizen
Alle zyklischen Matrizen $C \in K^{n \times n}$ 
$$C:=\begin{pmatrix} 
    c_0      &c_{1}    &c_{2}   &\ldots    &c_{n-1}\ 
    c_{n-1}  &c_0      &c_{1}   &\ldots    &c_{n-2}\ 
    c_{n-2}  &c_{n-1}  &c_0     &\ldots    &c_{n-3}\ 
             &\ddots   &\ddots  &\ddots\ 
    c_{1}    &c_{2}    &c_{3}   &\ldots    &c_0\end{pmatrix}. $$
haben die gleichen Eigenvektoren
$$v^{(m)} = \frac{1}{\sqrt{n}} \begin{pmatrix}1 \ e^{-2\pi i m/n} \ \vdots \  e^{-2\pi i m(n-1)/n} \end{pmatrix} .$$
und die gleiche Formel funktioniert für die Berechnung der Eigenwerte aller zyklischen Matrizen
$$ \lambda_m = \sum_{k=0}^{n-1}c_k e^{-2\pi i m k /n}. $$
Eine weitere Möglichkeit die Eigenwerte zu berechnen ist es die generierende Funktion 
$$f(\lambda) = \sum_{k=-s}^{s} t_k \exp(-i\lambda)^k.$$
zu verwenden. Z.B. für die zyklischen Matrizen
$$C_{n}(f)=\begin{pmatrix} 
t_0    & \ldots & t_{-s}&  \mathbf{0}m     &t{s} &\ldots & t_{1}   \
\vdots &        &       & \ddots&      &\ddots & \vdots  \
t_{s}  &        &       &       &      &       & t_{s}   \
 \mathbf{0}m^T       &\ddots  &       &   t_0 &      & \ddots&     \mathbf{0}_m^T     \ 
t{-s} &        &       &       &      &       & t_{-s}  \
\vdots & \ddots &       & \ddots&      &       & \vdots  \
t_{-1} & \ldots & t_{-s}&  \mathbf{0}_m      &  t_s & \ldots & t_0
\end{pmatrix} \in K^{n \times n}, \mbox{mit}\; n > 2s+1.$$
So sind die Eigenwerte von $C_n(f)$ $\psi_{n,j} = f(\frac{2 \pi j}{n})$.
Was einem sofort auch Abschätzungen für die Eigenwerte gibt und falls $f$ reelwertig und $F$ zusätzlich stetig ist gilt
$$ \lim_{n \to \infty} \frac{1}{n} \sum_{j=0}^{n-1} F(\psi_{n,j}) = \frac{1}{2\pi}\int_{0}^{2\pi} F(f(\lambda)) \mathrm{d}\lambda .$$
Toeplitz Matrizen
Zu jeder Toeplitz Matrix 
$$T_{n+1}= \begin{pmatrix} 
t_0     & \ldots  & t_{-n}  \
\vdots  & \ddots  & \vdots  \
t_{n} & \ldots  & t_{0}
\end{pmatrix} $$ 
ergibt sich eine generierende Funktion 
$$f(\lambda) = \sum_{k=-n}^{n} t_k \exp(-i\lambda)^k.$$ 
Für diese gilt dann
$$ T_n(f) = {\frac{1}{2 \pi} \int_0^{2\pi} f(\lambda) e^{-i(k-j)\lambda} \mathrm{d}\lambda ; k,j = 0,1,\ldots,n-1}.$$
Dieses $f$ kann genutzt werden um die Norm und Eigenwerte von $T_n$ abzuschätzen.
$$ \| T_n(f) \|\leq 2 M_{|f|}, \quad \mbox{wobei}\quad  M_{|f|} = \max_{\lambda \in K} |f(\lambda)|.$$
Diese Matrizen können künstlich vergrößert werden ohne grundlegende Eigenschaften zu ändern.
Es entstehen zusätzlich sogar asymptotische Eigenschaften. Sei z.B. 
$$T_{n}(f)=\begin{pmatrix} 
t_0    & \ldots & t_{-s}& & & & \
\vdots &        &       & \ddots & & &  \
t_{s}  &        &       &       &      &       &   \
       &\ddots  &       &   t_0 &      & \ddots&       \ 
& & &       &      &       & t_{-s}  \
& & & \ddots & &       & \vdots  \
 & & &       &  t_s & \ldots & t_0
\end{pmatrix} \in K^{n \times n}, \mbox{mit}\; n > 2s+1.$$
so gilt 
$$ \lim_{n \to \infty} \frac{1}{n} \sum_{j=0}^{n-1} F(\tau_{n,j}) = \frac{1}{2\pi}\int_{0}^{2\pi} F(f(\lambda)) \mathrm{d}\lambda.$$
Asymptotisch äquivalente Folgen von Matrizen
Seien ${A_n}$ und ${B_n}$ Folgen von $n\times n$ Matrizen, welche 
beschränkt bzgl. der starken Norm sind:
$$ \|A_n\|,\|B_n\| \leq M \le \infty, n=1,2,\ldots $$ 
und bzgl. der schwachen Norm konvergieren
$$\lim_{n \to \infty} |A_n -B_n| = 0.$$
Wir nennen diese Folgen asymptotisch äquivalent und notieren dies als $A_n \sim B_n$. 
Für ${A_n}$ , ${B_n}$ und ${C_n}$, 
welche jeweils die Eigenwerte ${\alpha_{n,i}}$,${\beta_{n,i}}$ und ${\zeta_{n,i}}$ haben gelten
folgende Zusammenhänge.
Wenn $A_n \sim B_n$, dann $\lim_{n \to \infty} |A_n| = \lim_{n \to \infty} |B_n| $
Wenn $A_n \sim B_n$ und $B_n \sim C_n$, dann  $A_n \sim C_n$
Wenn $A_nB_n \sim C_n$ und $\|A_n^{-1}\|\leq K \le \infty$, dann gilt $B_n \sim A_n^{-1}C_n$
Wenn $A_n \sim B_n$, dann $\exists -\infty \le m,M\le \infty$, s.d. $m\leq \alpha_{n,i}, \beta_{n,i}\leq M \; \forall n\geq 1 \mbox{und}\; k\geq 0$  
Wenn $A_n \sim B_n$, dann gilt $ \lim_{n \to \infty} \frac{1}{n} \sum_{k=0}^{n-1} (\alpha_{n,k}^s - \beta_{n,k}^s) = 0$
Theory in action
Systemmatrizen
Nutzen Sie pyMG um die Systemmatrix für das Helmholtz-Problem in 1D aufzustellen für gegebene Parameter $ n$ und $\sigma$.
End of explanation
matrix_plot(to_dense(system_matrix_hh1d(10,0)))
matrix_plot(to_dense(system_matrix_hh1d_periodic(10,0)))
def plot_3_eigvalueplots(A_p,A_z,A_m):
    eig_p.append(sp.linalg.eigvals(to_dense(A_p)))
    eig_z.append(sp.linalg.eigvals(to_dense(A_z)))
    eig_m.append(sp.linalg.eigvals(to_dense(A_m)))
    
    real_part_p = np.real(eig_p[-1])
    img_part_p = np.imag(eig_p[-1])
    
    real_part_z = np.real(eig_z[-1])
    img_part_z = np.imag(eig_z[-1])
    
    real_part_m = np.real(eig_m[-1])
    img_part_m = np.imag(eig_m[-1])
    fig1, (ax1, ax2, ax3) = plt.subplots(ncols=3,figsize=(15,3))
    
    ax1.plot(real_part_p,img_part_p,'ro')
    ax1.set_xlabel("real part")
    ax1.set_ylabel("img part")
    ax1.set_title('eigenvalues')
    
    ax2.plot(real_part_z,img_part_z,'bo')
    ax2.set_xlabel("real part")
    ax2.set_ylabel("img part")
    ax2.set_title('eigenvalues')
    
    ax3.plot(real_part_m,img_part_m,'go')
    ax3.set_xlabel("real part")
    ax3.set_ylabel("img part")
    ax3.set_title('eigenvalues')
    
    fig1.tight_layout()
    plt.show()
def plot_2_eigvalueplots(A_p,A_z):
    eig_p.append(sp.linalg.eigvals(to_dense(A_p)))
    eig_z.append(sp.linalg.eigvals(to_dense(A_z)))
    
    real_part_p = np.real(eig_p[-1])
    img_part_p = np.imag(eig_p[-1])
    
    real_part_z = np.real(eig_z[-1])
    img_part_z = np.imag(eig_z[-1])
    
    fig1, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,3))
    
    ax1.plot(real_part_p,img_part_p,'ro')
    ax1.set_xlabel("real part")
    ax1.set_ylabel("img part")
    ax1.set_title('eigenvalues')
    
    ax2.plot(real_part_z,img_part_z,'bo')
    ax2.set_xlabel("real part")
    ax2.set_ylabel("img part")
    ax2.set_title('eigenvalues')
    
    fig1.tight_layout()
    plt.show()
Explanation: Plotten Sie mithilfe von matrix_plot die Systemmatrizen für $\sigma = 0$ und $n=10$.
End of explanation
eig_p=[]
eig_m=[]
eig_z=[]
for n in [5,10,20]:
    A_p = system_matrix_hh1d(n,100.0)
    A_z = system_matrix_hh1d(n,0.0)
    A_m = system_matrix_hh1d(n,-100.0)
    
    plot_3_eigvalueplots(A_p, A_z, A_m)
Explanation: Aufgabe: Plotten Sie mithilfe von plot_3_eigvalueplots die Eigenwerte der Systemmatrix für $n \in [5,10,20]$ und $\sigma = 100$,$\sigma = -100$ und $\sigma = 0$.
End of explanation
def iteration_matrix_wjac(n, sigma, periodic=True):
    
    if periodic:
        A = system_matrix_hh1d_periodic(n,sigma)
    else:
        A = system_matrix_hh1d(n,sigma)
        
    wjac = WeightedJacobi(A, 2.0/3.0)
    P_inv = wjac.Pinv
    return np.eye(n) - P_inv.dot(A)
matrix_plot(iteration_matrix_wjac(10,0))
n=10
for sigma in [100,0,-100]:
    plot_2_eigvalueplots(iteration_matrix_wjac(n, sigma,periodic=True),iteration_matrix_wjac(n, sigma,periodic=False))
sigma_range = np.linspace(-100,100,100)
sr_wjac_periodic = map(lambda sig : spec_rad(iteration_matrix_wjac(n, sig,periodic=True)), sigma_range)
sr_wjac = map(lambda sig : spec_rad(iteration_matrix_wjac(n, sig,periodic=False)), sigma_range)
# Achsen festhalten
fig1, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,4))
ax1.plot(sigma_range, sr_wjac_periodic,'k-')
ax1.set_xlabel('$\sigma$')
ax1.set_ylabel("spectral radius")
ax1.set_title('periodic')
ax2.plot(sigma_range, sr_wjac,'k-')
ax2.set_xlabel('$\sigma$')
ax2.set_ylabel("spectral radius")
ax2.set_title('non-periodic')
fig1.tight_layout()
plt.show()
Explanation: Frage: Wie unterscheiden sich die Spektren der verschiedenen Systemmatrizen?
Iterationsmatrizen des Glätters
Weitaus spannender sind die Spektren der Iterationsmatrizen eines Glätters spannender.
End of explanation
def iteration_matrix_gs(n, sigma, periodic=True):
    
    if periodic:
        A = system_matrix_hh1d_periodic(n,sigma)
    else:
        A = system_matrix_hh1d(n,sigma)
        
    gs = GaussSeidel(A)
    P_inv = gs.Pinv
    return np.eye(n) - P_inv.dot(A)
matrix_plot(iteration_matrix_gs(10,0,True))
n=10
for sigma in [100,0,-100]:
    plot_2_eigvalueplots(iteration_matrix_gs(n, sigma,periodic=True),iteration_matrix_gs(n, sigma,periodic=False))
sr_gs_periodic = map(lambda sig : spec_rad(iteration_matrix_gs(n, sig,periodic=True)), sigma_range)
sr_gs = map(lambda sig : spec_rad(iteration_matrix_gs(n, sig,periodic=False)), sigma_range)
# Achsen festhalten
fig1, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,4))
ax1.plot(sigma_range, sr_gs_periodic,'k-')
ax1.set_xlabel('$\sigma$')
ax1.set_ylabel("spectral radius")
ax1.set_title('periodic')
ax2.plot(sigma_range, sr_gs,'k-')
ax2.set_xlabel('$\sigma$')
ax2.set_ylabel("spectral radius")
ax2.set_title('non-periodic')
fig1.tight_layout()
plt.show()
Explanation: Frage : Wie verhalten sich die Spektren für das periodische Problem zu den Problemen mit Dirichletrandbedingungen?
Schreiben Sie eine Funktion um die Iterationsmatrix für Gauß-Seidel abhängig von $\sigma$ und $n$ zu berechnen und finden Sie heraus wie sich der Spektralradius für verschiedene $\sigma$ und den periodischen, sowie nicht periodischen Fall verhält.
End of explanation
plt.semilogy(sigma_range, np.asarray(sr_gs_periodic)-np.asarray(sr_gs),'k-')
Explanation: Frage: Was sieht man und was ist eigentlich mit der Differenz zwischen den beiden?
End of explanation
def transformation_matrix_fourier_basis(N):
    psi = np.zeros((N,N),dtype=np.complex128)
    for i in range(N):
        for j in range(N):
            psi[i,j] = np.exp(2*np.pi*1.0j*j*i/N)
    return psi/np.sqrt(N)
PSI_trafo = transformation_matrix_fourier_basis(n)
PSI_trafo_inv = sp.linalg.inv(PSI_trafo)
wjac_trafo = np.dot(PSI_trafo_inv, np.dot(iteration_matrix_wjac(n,0),PSI_trafo))
gs_trafo = np.dot(PSI_trafo_inv, np.dot(iteration_matrix_gs(n,0),PSI_trafo))
matrix_plot(np.real(wjac_trafo))
matrix_plot(np.real(gs_trafo))
Explanation: Nur die hohen Frequenzen
End of explanation
def plot_fourier_transformed(A):
    A = to_dense(A)
    n = A.shape[0]
    PSI_trafo = transformation_matrix_fourier_basis(n)
    PSI_trafo_inv =  sp.linalg.inv(PSI_trafo)
    A_traf = np.dot(PSI_trafo_inv, np.dot(A,PSI_trafo))
    matrix_row_plot([A,np.abs(A_traf)])
def get_high_theta_eigvals(A, plot=False):
    A = to_dense(A)
    n = A.shape[0]
    PSI_trafo = transformation_matrix_fourier_basis(n)
    PSI_trafo_inv =  sp.linalg.inv(PSI_trafo)
    A_traf = np.dot(PSI_trafo_inv, np.dot(A,PSI_trafo))
    if plot:
        matrix_plot(np.abs(A_traf))
    eigvals =  np.asarray(map(lambda k : A_traf[k,k],range(n)))
    return eigvals[np.floor(n/4):np.ceil(3.0*n/4)]
def get_low_theta_eigvals(A, plot=False):
    A = to_dense(A)
    n = A.shape[0]
    PSI_trafo = transformation_matrix_fourier_basis(n)
    PSI_trafo_inv =  sp.linalg.inv(PSI_trafo)
    A_traf = np.dot(PSI_trafo_inv, np.dot(A,PSI_trafo))
    if plot:
        matrix_plot(np.abs(A_traf))
    eigvals =  np.asarray(map(lambda k : A_traf[k,k],range(n)))
    return np.hstack([eigvals[:np.floor(n/4)],eigvals[np.ceil(3.0*n/4):]])
high_eigvals_periodic = np.abs(get_high_theta_eigvals(iteration_matrix_gs(20,0,True),True))
high_eigvals_dirichlet = np.abs(get_high_theta_eigvals(iteration_matrix_gs(20,0,False),True))
print high_eigvals_periodic
high_eigvals_dirichlet - high_eigvals_periodic
Explanation: Frage: Was ist hier passiert?
Die hohen Eigenwerte extrahiert man nun durch auspicken der richtigen Diagonalwerte nach der Transformation, falls die Matrix zyklisch ist. Wir giessen
das ganze in eine Funktion.
End of explanation
do_plot = False
for n in [5,10,50]:
    for sig in [-100,0,100]:
        h_eigs_periodic = np.abs(get_high_theta_eigvals(iteration_matrix_gs(n,sig,True),do_plot))
        h_eigs_dirichlet = np.abs(get_high_theta_eigvals(iteration_matrix_gs(n,sig,False),do_plot))
        print "Sigma :\t\t",sig
        print "n:\t\t",n
#         print "\n",h_eigs_periodic-h_eigs_dirichlet,"\n"
        print "sum:\t\t", np.sum(np.abs(h_eigs_periodic-h_eigs_dirichlet))
        print "max:\t\t", np.max(np.abs(h_eigs_periodic-h_eigs_dirichlet)),"\n"
Explanation: Aufgabe: Schauen Sie sich die Differenz (periodisch und nicht-periodisch) der Eigenwerte der Gauss-Seidel Iterationsmatrix, die mit den hohen Frequenzen assoziiert werden können , für $n\in {5,10,50}$ und $\sigma \in {-100,0,100}$ an. Was fällt auf?
End of explanation
from project.linear_transfer import LinearTransfer
from project.linear_transfer_periodic import LinearTransferPeriodic
Explanation: Zweigitter-Iterationsmatrix
End of explanation
def coarse_grid_correction(n,nc, sigma):
    A_fine = to_dense(system_matrix_hh1d(n,sigma))
    A_coarse = to_dense(system_matrix_hh1d(nc,sigma))
    A_coarse_inv = sp.linalg.inv(A_coarse)
    lin_trans = LinearTransfer(n, nc)
    prolong = to_dense(lin_trans.I_2htoh)
    restrict = to_dense(lin_trans.I_hto2h)
    return np.eye(n)- np.dot(prolong.dot(A_coarse_inv.dot(restrict)), A_fine)
Explanation: Im Folgenden werden wir nun mithilfe des pymg frameworks die Zweigitter-Iterationsmatrix für ein einfaches Multigrid 
aufstellen. Wir beginnen mit der Grobgitterkorrektur.
End of explanation
plot_fourier_transformed(coarse_grid_correction(31,15,-1000))
plot_fourier_transformed(coarse_grid_correction(31,15,0))
plot_fourier_transformed(coarse_grid_correction(31,15,1000))
plot_3_eigvalueplots(coarse_grid_correction(31,15,-1000),coarse_grid_correction(31,15,0),coarse_grid_correction(31,15,100))
Explanation: Buntebilderaufgabe: Nutze matrix_plot um für $n=31$, $n_c=15$ und verschiedene $\sigma\in[-1000,1000]$ um die Grobgitterkorrekturiterationsmatrizen und deren Fourier-transformierten zu plotten.
End of explanation
def coarse_grid_correction_periodic(n,nc, sigma):
    A_fine = to_dense(system_matrix_hh1d_periodic(n,sigma))
    A_coarse = to_dense(system_matrix_hh1d_periodic(nc,sigma))
    A_coarse_inv = sp.linalg.inv(A_coarse)
    lin_trans = LinearTransferPeriodic(n, nc)
    prolong = to_dense(lin_trans.I_2htoh)
    restrict = to_dense(lin_trans.I_hto2h)
    return np.eye(n)- np.dot(prolong.dot(A_coarse_inv.dot(restrict)), A_fine)
matrix_plot(coarse_grid_correction_periodic(31,15,-1000))
matrix_plot(coarse_grid_correction_periodic(31,15,-0.00))
matrix_plot(coarse_grid_correction_periodic(31,15,1000))
plot_fourier_transformed(coarse_grid_correction_periodic(32,16,-1000))
plot_fourier_transformed(coarse_grid_correction_periodic(32,16,-0.00))
plot_fourier_transformed(coarse_grid_correction_periodic(32,16,1000))
prolong = to_dense(LinearTransferPeriodic(16, 8).I_2htoh)
restrict = to_dense(LinearTransferPeriodic(16, 8).I_hto2h)
matrix_plot(to_dense(LinearTransferPeriodic(16, 8).I_2htoh))
matrix_plot(to_dense(LinearTransferPeriodic(16, 8).I_hto2h))
matrix_plot(prolong.dot(restrict))
def mat_of_interest(n):
    prolong = to_dense(LinearTransferPeriodic(n, n/2).I_2htoh)
    restrict = to_dense(LinearTransferPeriodic(n, n/2).I_hto2h)
    return prolong.dot(restrict)
prolong = to_dense(LinearTransferPeriodic(15, 7).I_2htoh)
restrict = to_dense(LinearTransferPeriodic(15, 7).I_hto2h)
matrix_plot(to_dense(LinearTransferPeriodic(15, 7).I_2htoh))
matrix_plot(to_dense(LinearTransferPeriodic(15, 7).I_hto2h))
matrix_plot(prolong.dot(restrict))
def two_grid_it_matrix(n,nc, sigma, nu1=3,nu2=3,typ='wjac'):
    cg = coarse_grid_correction(n,nc,sigma)
    if typ is 'wjac':
        smoother = iteration_matrix_wjac(n,sigma, periodic=False)
    if typ is 'gs':
        smoother = iteration_matrix_gs(n,sigma, periodic=False)
    
    pre_sm = matrix_power(smoother, nu1)
    post_sm = matrix_power(smoother, nu2)
    
    return pre_sm.dot(cg.dot(post_sm))    
Explanation: Aufgabe: Schreiben Sie die Grobgitterkorrektur für das periodische Problem und plotten Sie nochmal für verschiedene $\sigma$.
Frage: Was genau passiert bei $\sigma = 0$ und in der Nähe davon?
End of explanation
plot_fourier_transformed(two_grid_it_matrix(15,7,-100,typ='wjac'))
plot_fourier_transformed(two_grid_it_matrix(15,7,0,typ='wjac'))
plot_fourier_transformed(two_grid_it_matrix(15,7,100,typ='wjac'))
plot_fourier_transformed(two_grid_it_matrix(15,7,-100,typ='gs'))
plot_fourier_transformed(two_grid_it_matrix(15,7,0,typ='gs'))
plot_fourier_transformed(two_grid_it_matrix(15,7,100,typ='gs'))
eig_p=[]
eig_m=[]
eig_z=[]
for n,nc in zip([7,15,31],[3,7,15]):
    A_p = two_grid_it_matrix(n,nc,100.0)
    A_z = two_grid_it_matrix(n,nc,0.0)
    A_m = two_grid_it_matrix(n,nc,-100.0)
    
    plot_3_eigvalueplots(A_p, A_z, A_m)
sr_2grid_var_sigma = map(lambda sig : spec_rad(two_grid_it_matrix(15,7,sig)), sigma_range)
plt.semilogy(sigma_range, sr_2grid_var_sigma,'k-')
plt.title('$n_f = 15, n_c = 7$')
plt.xlabel('$\sigma$')
plt.ylabel("spectral radius")
nf_range = map(lambda k: 2**k-1,range(3,10))
nc_range = map(lambda k: 2**k-1,range(2,9))
sr_2grid_m1000 = map(lambda nf,nc : spec_rad(two_grid_it_matrix(nf,nc,-1000)), nf_range, nc_range)
sr_2grid_0 = map(lambda nf,nc : spec_rad(two_grid_it_matrix(nf,nc,0)), nf_range, nc_range)
sr_2grid_p1000 = map(lambda nf,nc : spec_rad(two_grid_it_matrix(nf,nc,1000)), nf_range, nc_range)
plt.semilogy(nf_range, sr_2grid_m1000,'k-',nf_range, sr_2grid_0,'k--',nf_range, sr_2grid_p1000,'k:')
plt.xlabel('$n_f$')
plt.ylabel("spectral radius")
plt.legend(("$\sigma = -1000$","$\sigma = 0$","$\sigma = 1000$"),'upper right',shadow = True)
Explanation: Buntebilderaufgabe:Nutze matrix_plot um für $n=31$, $n_c=15$ und verschiedene $\sigma\in[-1000,1000]$ um die Zweigittermatrizen und deren Fourier-transformierten zu plotten.
End of explanation
def two_grid_it_matrix_periodic(n,nc, sigma, nu1=3,nu2=3,typ='wjac'):
    cg = coarse_grid_correction_periodic(n,nc,sigma)
    if typ is 'wjac':
        smoother = iteration_matrix_wjac(n,sigma, periodic=True)
    if typ is 'gs':
        smoother = iteration_matrix_gs(n,sigma, periodic=True)
    
    pre_sm = matrix_power(smoother, nu1)
    post_sm = matrix_power(smoother, nu2)
    return pre_sm.dot(cg.dot(post_sm))    
plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,-100,typ='wjac'))
plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,0.01,typ='wjac'))
plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,100,typ='wjac'))
plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,-100,typ='gs'))
plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,-0.01,typ='gs'))
plot_fourier_transformed(two_grid_it_matrix_periodic(16,8,100,typ='gs'))
sr_2grid_var_sigma_periodic = map(lambda sig : spec_rad(two_grid_it_matrix_periodic(16,8,sig)), sigma_range)
plt.plot(sigma_range, np.asarray(sr_2grid_var_sigma)-np.asarray(sr_2grid_var_sigma_periodic),'k-')
plt.title('Differenz periodisch und nicht periodisch')
plt.xlabel('$\sigma$')
plt.ylabel("spectral radius")
nf_range = map(lambda k: 2**k,range(3,10))
nc_range = map(lambda k: 2**k,range(2,9))
sr_2grid_m1000_p = map(lambda nf,nc : spec_rad(two_grid_it_matrix_periodic(nf,nc,-1000)), nf_range, nc_range)
sr_2grid_0_p = map(lambda nf,nc : spec_rad(two_grid_it_matrix_periodic(nf,nc,0.01)), nf_range, nc_range)
sr_2grid_p1000_p = map(lambda nf,nc : spec_rad(two_grid_it_matrix_periodic(nf,nc,1000)), nf_range, nc_range)
plt.semilogy(nf_range, sr_2grid_m1000_p,'k-',nf_range, sr_2grid_0_p,'k--',nf_range, sr_2grid_p1000_p,'k:')
plt.xlabel('$n_f$')
plt.ylabel("spectral radius")
plt.legend(("$\sigma = -1000$","$\sigma = 0$","$\sigma = 1000$"),'upper right',shadow = True)
Explanation: Aufgabe: Schreiben Sie eine Funktion two_grid_it_matrix_periodic für den periodischen Fall und führen plotten Sie den Spektralradius über $\sigma$ und den Spektralradius über $n$ für 3 verschiedene $\sigma$.
End of explanation
plt.semilogy(nf_range, np.abs(np.asarray(sr_2grid_m1000_p) - np.asarray(sr_2grid_m1000)),'k-',
             nf_range, np.abs(np.asarray(sr_2grid_0_p) - np.asarray(sr_2grid_0)),'k--',
             nf_range, np.abs(np.asarray(sr_2grid_p1000_p) - np.asarray(sr_2grid_p1000)),'k:')
plt.xlabel('$n_f$')
plt.ylabel("spectral radius")
plt.legend(("$\sigma = -1000$","$\sigma = 0$","$\sigma = 1000$"),'upper right',shadow = True)
Explanation: Aufgabe: Plotten sie die Differenz zwischen periodischem und nicht periodischem Fall.
End of explanation
def hs_norm(A):
    n = A.shape[0]
    return sp.linalg.norm(A,'fro')/np.sqrt(n)
Explanation: Asymptotische Äquivalenz  zwischen periodisch und nicht-periodisch
Wir sehen, dass die Spektralradiien auf den ersten Blick gut übereinstimmen. Wir wollen nun empirisch ergründen ob die Matrizenklassen der periodischen und nicht periodischen Fällen möglicherweise zueinander asymptotisch äquivalent sind.
Aufgabe: 
Schreiben sie eine Funktion hs_norm, welche die Hilbert-Schmidt Norm berechnet.
End of explanation
n_range = np.arange(10,100)
hs_sysmat_m1000 = map(lambda n: hs_norm(to_dense(system_matrix_hh1d(n,-1000))-to_dense(system_matrix_hh1d_periodic(n,-1000))),n_range)
hs_sysmat_0 = map(lambda n: hs_norm(to_dense(system_matrix_hh1d(n,0.001))-to_dense(system_matrix_hh1d_periodic(n,0.001))),n_range)
hs_sysmat_p1000 = map(lambda n: hs_norm(to_dense(system_matrix_hh1d(n,1000))-to_dense(system_matrix_hh1d_periodic(n,1000))),n_range)
plt.plot(hs_sysmat_m1000)
plt.plot(hs_sysmat_0)
plt.plot(hs_sysmat_p1000)
Explanation: Aufgabe: 
Überprüfen Sie empirisch ob die
Systemmatrizenklassen 
Glättungsiterationsmatrizenklassen 
Grobgitterkorrekturmatrizenklassen 
Zweigitteriterationsmatrizenklassen 
asymptotisch äquivalent sind für $\sigma = { -1000, 0.001, 1000 }$.
Systemmatrizen:
End of explanation
n_range = 2**np.arange(1,11)
hs_wjac_m1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_wjac(n,-1000))-to_dense(iteration_matrix_wjac(n,-1000,False))),n_range)
hs_wjac_0 = map(lambda n: hs_norm(to_dense(iteration_matrix_wjac(n,0))-to_dense(iteration_matrix_wjac(n,0,False))),n_range)
hs_wjac_p1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_wjac(n,1000))-to_dense(iteration_matrix_wjac(n,1000,False))),n_range)
plt.plot(n_range, hs_wjac_m1000)
plt.plot(n_range, hs_wjac_0)
plt.plot(n_range, hs_wjac_p1000)
Explanation: Glättung:
Jacobi
End of explanation
n_range = 2**np.arange(1,11)
hs_gs_m1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_gs(n,-1000))-to_dense(iteration_matrix_gs(n,-1000,False))),n_range)
hs_gs_0 = map(lambda n: hs_norm(to_dense(iteration_matrix_gs(n,0))-to_dense(iteration_matrix_gs(n,0,False))),n_range)
hs_gs_p1000 = map(lambda n: hs_norm(to_dense(iteration_matrix_gs(n,1000))-to_dense(iteration_matrix_gs(n,1000,False))),n_range)
plt.plot(n_range, hs_gs_m1000)
plt.plot(n_range, hs_gs_0)
plt.plot(n_range, hs_gs_p1000)
Explanation: Gauss-Seidel
End of explanation
def einmal_einpacken(A):
    return np.r_[[np.zeros(A.shape[0]+1)],np.c_[np.zeros(A.shape[0]),A]]
n_f_range = 2**np.arange(3,10)
n_c_range = 2**np.arange(2,9)
hs_cgc_m1000 =  map(lambda nf,nc: hs_norm(einmal_einpacken(coarse_grid_correction(nf-1,nc-1,-1000))-coarse_grid_correction_periodic(nf,nc,-1000)),n_f_range ,n_c_range)
hs_cgc_0 =  map(lambda nf,nc: hs_norm(einmal_einpacken(coarse_grid_correction(nf-1,nc-1,0))-coarse_grid_correction_periodic(nf,nc,0.001)),n_f_range ,n_c_range)
hs_cgc_p1000 =  map(lambda nf,nc: hs_norm(einmal_einpacken(coarse_grid_correction(nf-1,nc-1,1000))-coarse_grid_correction_periodic(nf,nc,1000)),n_f_range ,n_c_range)
plt.semilogy(n_f_range, hs_cgc_m1000)
plt.semilogy(n_f_range, hs_cgc_0)
plt.semilogy(n_f_range, hs_cgc_p1000)
# plt.semilogy(n_f_range, 1/np.sqrt(n_f_range))
Explanation: Grobgitterkorrektur
Hier trifft man mal wieder auf das Problem, dass die Freiheitsgrade im periodischen und nicht periodischen Fall unterschiedlich verteilt sind. Für einen Vergleich wird die Matrix um eine Nullzeile und Nullspalte erweitert.
End of explanation
n_f_range = 2**np.arange(3,12)
n_c_range = 2**np.arange(2,11)
hs_2grid_m1000 =  map(lambda nf,nc: hs_norm(
        einmal_einpacken(two_grid_it_matrix(nf-1,nc-1,-1000))-two_grid_it_matrix_periodic(nf,nc,-1000))
                    ,n_f_range ,n_c_range)
hs_2grid_0 =  map(lambda nf,nc: hs_norm(
        einmal_einpacken(two_grid_it_matrix(nf-1,nc-1,0.001))-two_grid_it_matrix_periodic(nf,nc,0.001))
                ,n_f_range ,n_c_range)
hs_2grid_p1000 =  map(lambda nf,nc: hs_norm(
        einmal_einpacken(two_grid_it_matrix(nf-1,nc-1,1000))-two_grid_it_matrix_periodic(nf,nc,1000))
                    ,n_f_range ,n_c_range)
plt.semilogy(n_f_range, hs_2grid_m1000)
plt.semilogy(n_f_range, hs_2grid_0)
plt.semilogy(n_f_range, hs_2grid_p1000)
plt.semilogy(n_f_range, 1/np.sqrt(n_f_range)*30)
Explanation: Zweigitter
End of explanation
<END_TASK>