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.
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
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).$$
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.
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 .$$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.$$
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.
Nutzen Sie pyMG um die Systemmatrix für das Helmholtz-Problem in 1D aufzustellen für gegebene Parameter $ n$ und $\sigma$.
In [1]:
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
In [2]:
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
In [3]:
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))))
Plotten Sie mithilfe von matrix_plot die Systemmatrizen für $\sigma = 0$ und $n=10$.
In [4]:
matrix_plot(to_dense(system_matrix_hh1d(10,0)))
matrix_plot(to_dense(system_matrix_hh1d_periodic(10,0)))
In [5]:
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()
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$.
In [6]:
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)
Frage: Wie unterscheiden sich die Spektren der verschiedenen Systemmatrizen?
Weitaus spannender sind die Spektren der Iterationsmatrizen eines Glätters spannender.
In [7]:
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)
In [8]:
matrix_plot(iteration_matrix_wjac(10,0))
In [110]:
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))
In [10]:
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)
In [11]:
# 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()
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.
In [12]:
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)
In [13]:
matrix_plot(iteration_matrix_gs(10,0,True))
In [14]:
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))
In [15]:
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)
In [16]:
# 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()
Frage: Was sieht man und was ist eigentlich mit der Differenz zwischen den beiden?
In [17]:
plt.semilogy(sigma_range, np.asarray(sr_gs_periodic)-np.asarray(sr_gs),'k-')
Out[17]:
In [18]:
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)
In [19]:
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))
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.
In [20]:
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)])
In [21]:
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):]])
In [111]:
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))
In [112]:
print high_eigvals_periodic
In [24]:
high_eigvals_dirichlet - high_eigvals_periodic
Out[24]:
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?
In [25]:
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"
In [26]:
from project.linear_transfer import LinearTransfer
from project.linear_transfer_periodic import LinearTransferPeriodic
Im Folgenden werden wir nun mithilfe des pymg frameworks die Zweigitter-Iterationsmatrix für ein einfaches Multigrid
aufstellen. Wir beginnen mit der Grobgitterkorrektur.
In [27]:
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)
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.
In [28]:
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))
In [113]:
plot_3_eigvalueplots(coarse_grid_correction(31,15,-1000),coarse_grid_correction(31,15,0),coarse_grid_correction(31,15,100))
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?
In [30]:
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)
In [31]:
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))
In [114]:
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))
In [33]:
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))
In [34]:
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)
In [35]:
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))
In [36]:
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))
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.
In [37]:
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'))
In [39]:
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'))
In [43]:
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)
In [86]:
sr_2grid_var_sigma = map(lambda sig : spec_rad(two_grid_it_matrix(15,7,sig)), sigma_range)
In [61]:
plt.semilogy(sigma_range, sr_2grid_var_sigma,'k-')
plt.title('$n_f = 15, n_c = 7$')
plt.xlabel('$\sigma$')
plt.ylabel("spectral radius")
Out[61]:
In [71]:
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)
In [77]:
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)
Out[77]:
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$.
In [78]:
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))
In [94]:
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'))
In [96]:
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'))
In [84]:
sr_2grid_var_sigma_periodic = map(lambda sig : spec_rad(two_grid_it_matrix_periodic(16,8,sig)), sigma_range)
In [116]:
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")
Out[116]:
In [98]:
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)
In [99]:
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)
Out[99]:
Aufgabe: Plotten sie die Differenz zwischen periodischem und nicht periodischem Fall.
In [101]:
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)
Out[101]:
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.
In [102]:
def hs_norm(A):
n = A.shape[0]
return sp.linalg.norm(A,'fro')/np.sqrt(n)
Aufgabe: Überprüfen Sie empirisch ob die
asymptotisch äquivalent sind für $\sigma = \{ -1000, 0.001, 1000 \}$.
Systemmatrizen:
In [107]:
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)
In [108]:
plt.plot(hs_sysmat_m1000)
plt.plot(hs_sysmat_0)
plt.plot(hs_sysmat_p1000)
Out[108]:
Glättung:
Jacobi
In [132]:
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)
In [133]:
plt.plot(n_range, hs_wjac_m1000)
plt.plot(n_range, hs_wjac_0)
plt.plot(n_range, hs_wjac_p1000)
Out[133]:
Gauss-Seidel
In [136]:
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)
In [137]:
plt.plot(n_range, hs_gs_m1000)
plt.plot(n_range, hs_gs_0)
plt.plot(n_range, hs_gs_p1000)
Out[137]:
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.
In [160]:
def einmal_einpacken(A):
return np.r_[[np.zeros(A.shape[0]+1)],np.c_[np.zeros(A.shape[0]),A]]
In [165]:
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)
In [171]:
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))
Out[171]:
Zweigitter
In [178]:
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)
In [183]:
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)
Out[183]: