In [1]:
import numpy as np
import numpy.linalg as linalg
import matplotlib.pyplot as pl
%matplotlib inline
In [2]:
def npower(matrix, k): # método da potência n-ésima para um
x = np.array([1, 0, 0])
print(0, x, x.transpose() @ matrix @ x)
for i in range(k):
x = matrix @ x
x_norm = linalg.norm(x)
xk = x / x_norm
eig = xk.transpose() @ A @ xk
print(i+1, xk, eig)
return xk, eig
In [3]:
def npower_eps(matrix, lambda_diff):
x0 = np.array([1, 0, 0])
xk = matrix @ x0
xk = xk / linalg.norm(xk)
k = 1
def eig(vector): return vector.transpose() @ matrix @ vector
print(k, eig(xk))
while np.abs(eig(xk) - eig(x0)) > lambda_diff:
x0 = xk
xk = matrix @ xk
xk /= linalg.norm(xk)
k += 1
if k%2 != 0: print(k, eig(xk))
return k, eig(xk), xk, eig(x0), x0, np.abs(eig(x0) - eig(xk))
In [4]:
A = np.array([[1., 1, 1/2], [1, 1, 1/4], [1/2, 1/4, 2]])
In [5]:
npower_eps(A, 1e-5)
Out[5]:
Cálculo utilizando o método sem critério de parada, mas em vez disso, em função do k Retornando:
In [6]:
npower(A, 6)
Out[6]:
In [7]:
npower(A, 12)
Out[7]:
In [8]:
linalg.eig(A)
Out[8]:
Exercício 2 - Método de Jacobi Cíclico
In [9]:
def jacobi(A_i, eps, nitmax):
A = np.copy(A_i) # para cortar dependências
m = len(A)
iteration = 0
Q = np.identity(m)
def off(mat):
off_sum = 0
for i in range(m):
for j in range(m):
if j != i: off_sum += mat[i, j]**2
return np.sqrt(off_sum)
def frobenius_norm(mat):
norm = 0
for i in range(m):
for j in range(m):
norm += mat[i, j]**2
return np.sqrt(norm)
while (off(A) > eps and iteration < nitmax):
j, k, ajk = 0, 0, 0.
for ji in range(m-1):
for ki in range(ji+1, m):
absjk = abs(A[ji, ki])
if absjk >= ajk:
j, k, ajk = ji, ki, absjk
def CSjk(mati, j, k):
mat = np.copy(mati)
if mat[j, j] == mat[k, k]:
C, S = np.cos(np.pi / 4), np.sin(np.pi / 4)
else:
tau = 2*mat[j, k] / (mat[k, k] - mat[j, j])
chi = 1 / np.sqrt(1 + tau**2)
C = np.sqrt((1 + chi) / 2)
S = np.sign(tau) * np.sqrt((1 - chi) / 2)
return C, S
C, S = CSjk(A, j, k)
A_l = np.zeros_like(A)
for r in range(m):
if r != j and r != k:
A_l[r, j] = C * A[r, j] - S * A[r, k]
A_l[j, r] = C * A[r, j] - S * A[r, k]
A_l[r, k] = S * A[r, j] + C * A[r, k]
A_l[k, r] = S * A[r, j] + C * A[r, k]
for s in range(m):
if s != j and s != k:
A_l[r, s] = np.copy(A[r, s])
A_l[j, j] = np.copy((C**2 * A[j, j]) + (S**2 * A[k, k]) - (2 * S * C * A[j, k]))
A_l[j, k] = S * C * (A[j, j] - A[k, k]) + ((C**2 - S**2) * A[j, k])
# A_l[j, k] = 0
A_l[k, j] = np.copy(A_l[j, k])
A_l[k, k] = np.copy((S**2 * A[j, j]) + (C**2 * A[k, k]) + (2 * S * C * A[j, k]))
A = A_l
Q_l = np.zeros_like(Q)
for r in range(m):
for s in range(m):
if s != j and s!= k:
Q_l[r, s] = np.copy(Q[r, s])
Q_l[r, j] = C * Q[r, j] - S * Q[r, k]
Q_l[r, k] = S * Q[r, j] + C * Q[r, k]
Q = Q_l
iteration += 1
D = Q.transpose().dot(A_i).dot(Q)
return A, off(A), D, Q, iteration, off(A_i), frobenius_norm(A_i)
In [10]:
a, oa, d, q, it, oi, fi = jacobi(A, 1e-4, 100)
print('D:')
print( d)
print('Q:')
print(q)
print('Iterações:', it)
In [11]:
oi**2 / fi**2
Out[11]:
In [12]:
oa**2 / fi**2
Out[12]:
In [13]:
linalg.eig(A)
Out[13]:
In [14]:
oi**2 / 8.625
Out[14]:
In [15]:
oa**2 / 8.625
Out[15]:
In [16]:
np.cos(np.pi/4)
Out[16]:
In [17]:
df_5pts = lambda f, x, h: (-3*f(x+4*h) + 16*f(x+3*h) - 36*f(x+2*h) + 48*f(x+h) - 25*f(x)) / (12*h)
df2_5pts = lambda f, x, h: df_5pts(df_5pts, x, h)
In [18]:
def ui(x):
if
In [ ]:
def schrodingerpvp(xmin, xmax, subint, eps):
h = (xmax - xmin) / subint
xi = [xmin + i for i in range(subint)]
# ui = [lambda x: u(x) for x in xi]
D = np.zeros((subint-1, subint-1))
for i in range(subint-2):
D[i, i] = (2 / h**2) + xi[i]**2
if i==0: D[i, 1] = - 1 / h**2
elif i== subint-2: D[i, i-1] = - 1 / h**2
else:
D[i, i+1] = -1 / h**2
D[i, i-1] = -1 / h**2
Dt, odt, lam, vec, it, od, frobd = jacobi(D, eps, 100000)
return lam, vec, D
In [ ]:
l, v, D = schrodingerpvp(-10, 10, 50, 1e-4)
In [ ]:
linalg.eigvals(D)