Commit 6a0be7a6 authored by Jonas David Grams's avatar Jonas David Grams
Browse files

Lösungshinweise Blatt 13 hinzugefügt

parent fcc15937
import numpy as np
import numpy.linalg as LA
n = 4
m = 100
B = np.random.randn(n, n)
A = B.T @ B
eigA, v = LA.eig(A)
eigA = np.sort(eigA)
eigA = eigA[::-1]
print(np.diag(eigA))
for i in range(m):
C = LA.cholesky(A)
A = C.T @ C
print(A)
print(np.diag(eigA))
import numpy as np
import numpy.linalg as LA
from scipy.linalg import lu
n = 4
m = 100
B = np.random.randn(n, n)
A = B.T @ B
eigA, v = LA.eig(A)
eigA = np.sort(eigA)
eigA = eigA[::-1]
print(np.diag(eigA))
for i in range(m):
PL, U = lu(A, permute_l=True)
A = U @ PL
print(A)
print(np.diag(eigA))
import numpy as np
from scipy.linalg import eigvals
import matplotlib.pyplot as plt
def qd_max (q, e, m, eps):
n = q.size
Q = np.zeros((n,m))
E = np.zeros((n,m))
Q[:,0] = q
E[:,0] = e
iter = 0
for i in range(1,m):
Q[0,i] = Q[0,i-1] + E[0,i-1]
for k in range(n-1):
E[k,i] = E[k,i-1]/Q[k,i] * Q[k+1,i-1]
Q[k+1,i] = Q[k+1,i-1] - E[k,i] + E[k+1,i-1]
if np.max(E[:-1,i]) < eps:
return iter, Q, E
iter += 1
return iter, Q, E
def qd_min (q, e, m, eps):
n = q.size
Q = np.zeros((n,m))
E = np.zeros((n,m))
Q[:,0] = q
E[:,0] = e
iter = 0
for i in range(1,m):
Q[0,i] = Q[0,i-1] + E[0,i-1]
for k in range(n-1):
E[k,i] = E[k,i-1]/Q[k,i] * Q[k+1,i-1]
Q[k+1,i] = Q[k+1,i-1] - E[k,i] + E[k+1,i-1]
if np.min(E[:-1,i]) < eps:
return iter, Q, E
iter += 1
def treiber():
n = 10
eps = 1e-16
m = 2000
q = (np.random.rand(n))
e = np.zeros(n)
e[:-1] = np.random.rand(n-1)
L = np.diag(q, 0) + np.diag(np.ones(n-1), 1)
R = np.diag(np.ones(n), 0) + np.diag(e[:-1], -1)
T = L @ R
eigs = eigvals(T)
eigs = np.sort(eigs)
m_min, Q_min, E_min = qd_min(q, e, m, eps)
m_max, Q_max, E_max = qd_max(q, e, m, eps)
idx_min = np.argsort(Q_min[:,m_min])
idx_max = np.argsort(Q_max[:,m_max])
brel_min = np.zeros((n,m))
babs_min = np.zeros((n,m))
brel_max = np.zeros((n,m))
babs_max = np.zeros((n,m))
print(Q_min[idx_min,m_min])
print(eigs)
for k in range(m):
babs_min[:,k] = np.abs(Q_min[idx_min,k] - eigs)
brel_min[:,k] = babs_min[:,k]/np.abs(eigs)
babs_max[:,k] = np.abs(Q_max[idx_max,k] - eigs)
brel_max[:,k] = babs_max[:,k]/np.abs(eigs)
print(babs_min[:,m_min])
plt.subplot(7, 1, (1,3))
plt.xlabel('Realteil')
plt.ylabel('Imaginärteil')
plt.title('min Abbruchkriterium')
plt.plot(np.real(eigs), np.imag(eigs), 'r+')
for j in range(n):
c = max(0, min(-np.log10(brel_min[j,m_min])/16,1))
plt.plot(np.real(Q_min[idx_min[j],m_min]), np.imag(Q_min[idx_min[j],m_min]), 'o', c=(0, 1-c, c, 1), fillstyle='none')
plt.legend(['tatsächliche Eigenwerte', 'QD-Werte'], loc='best')
plt.subplot(7, 1, (5,7))
plt.xlabel('Realteil')
plt.ylabel('Imaginärteil')
plt.title('max Abbruchkriterium')
plt.plot(np.real(eigs), np.imag(eigs), 'r+')
for j in range(n):
c = max(0, min(-np.log10(brel_max[j,m_max])/16,1))
plt.plot(np.real(Q_max[j,m_max]), np.imag(Q_max[j,m_max]), 'o', c=(0, 1-c, c, 1), fillstyle='none')
plt.legend(['tatsächliche Eigenwerte', 'QD-Werte'], loc='best')
plt.savefig('Eigenwerte.pdf', bbox_inches='tight')
plt.close()
plt.subplot(7,1,(1,3))
plt.yscale('log')
plt.title('min, absoluter Fehler')
for j in range(n):
plt.plot(babs_min[j,:m_min+1], label='${:.4f}$'.format(np.real(eigs[j])))
plt.ylim(1e-16, 1e3)
plt.xlabel('Schritte')
plt.ylabel('Fehler')
plt.legend(loc='best', ncol=4)
plt.subplot(7,1,(5,7))
plt.yscale('log')
plt.title('min, relativer Fehler')
for j in range(n):
plt.plot(brel_min[j,:m_min+1], label='${:.4f}$'.format(np.real(eigs[j])))
plt.ylim(1e-16, 1e3)
plt.xlabel('Schritte')
plt.ylabel('Fehler')
plt.legend(loc='best', ncol=4)
plt.savefig('Fehler_min.pdf', bbox_inches='tight')
plt.close()
plt.subplot(7,1,(1,3))
plt.yscale('log')
plt.title('max, absoluter Fehler')
for j in range(n):
plt.plot(babs_max[j,:m_max+1], label='${:.4f}$'.format(np.real(eigs[j])))
plt.ylim(1e-16, 1e3)
plt.xlabel('Schritte')
plt.ylabel('Fehler')
plt.legend(loc='best', ncol=4)
plt.subplot(7,1,(5,7))
plt.yscale('log')
plt.title('max, relativer Fehler')
for j in range(n):
plt.plot(brel_max[j,:m_max+1], label='${:.4f}$'.format(np.real(eigs[j])))
plt.ylim(1e-16, 1e3)
plt.xlabel('Schritte')
plt.ylabel('Fehler')
plt.legend(loc='best', ncol=4)
plt.savefig('Fehler_max.pdf', bbox_inches='tight')
plt.close()
if __name__ == '__main__':
treiber()
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment