Update PS7.py

PS7_Jung/PS7.py

@@ -12,42 +12,132 @@ import matplotlib.pyplot as plt
 T = np.matrix([[0,1/4,0,0,0], [1,0,1/2,1/2,0], [0,1/4,0,0,0], [0,1/4,0,0,0], [0,1/4,1/2,1/2,1]])
 eigenvalues, eigenvectors = np.linalg.eig(T)
 
-# plot spectrum of T
-real = [eigenvalue.real for eigenvalue in eigenvalues]
-imaginary = [eigenvalue.imag for eigenvalue in eigenvalues]
-plt.figure()
-plt.scatter(real, imaginary)
-plt.xlabel('Real')
-plt.ylabel('Imaginary')
-plt.tight_layout()
-plt.savefig('Spectrum_cat_mouse.png')
+# # plot spectrum of T
+# real = [eigenvalue.real for eigenvalue in eigenvalues]
+# imaginary = [eigenvalue.imag for eigenvalue in eigenvalues]
+# plt.figure()
+# plt.scatter(real, imaginary)
+# plt.xlabel('Real')
+# plt.ylabel('Imaginary')
+# plt.tight_layout()
+# plt.savefig('Spectrum_cat_mouse.png')
 
-# stationary distribution = last eigenvector
+# # stationary distribution = last eigenvector
 
-# c)
+# # c)
 av_columns = np.zeros(100)
 columns_indep = np.zeros((5,100))
 E_k = np.zeros(100)
+q_1 = np.array((1,0,0,0,0))
 
-
+print(f'{eigenvectors[:,0]=}')
 for k in range(1,101):
     T_k = np.linalg.matrix_power(T, k)
     # L2-Norm
     for i in range(5):
-        columns_indep[i,k-1] = np.linalg.norm(eigenvectors[0]-T_k[:,i])
+        columns_indep[i,k-1] = np.linalg.norm(eigenvectors[:,0]-T_k[:,i])
     av_columns[k-1] = np.average(columns_indep[:,k-1])
+    print(f'{(T_k.dot(q_1)-eigenvectors[:,0]).shape=}')
+    E_k[k-1] = np.linalg.norm((T_k.dot(q_1)).T-eigenvectors[:,0])
 
-    q_1 = np.array((1,0,0,0,0))
-    E_k[k-1] = np.linalg.norm(np.dot(T_k,q_1)-eigenvectors[0])
+# plt.figure()
+# for i in range(5):
+#     plt.plot(columns_indep[i,:], label=f'Column {i+1}')
+# plt.plot(av_columns, label='All')
+# plt.legend()
+# plt.tight_layout()
+# plt.savefig('Column_convergence.png')
 
-plt.figure()
-for i in range(5):
-    plt.plot(columns_indep[i,:], label=f'Column {i+1}')
-plt.plot(av_columns, label='All')
-plt.legend()
-plt.show()
+# # k =1,..,9
 
 plt.figure()
-plt.plot(np.log(E_k))
-plt.show()
+plt.plot(list(range(1,101)),np.log(E_k),'o')
+plt.xlim(1,101)
+plt.xlabel('k')
+plt.ylabel('$log(E_k)$')
+plt.tight_layout()
+plt.savefig('For_slope.png')
+
+## Problem 7.2
+# a)
+# T = np.matrix([[0,0,0,1], [1,0,0,0], [0,1,0,0], [0,0,1,0]])
+# eigenvalues, eigenvectors = np.linalg.eig(T)
+# real = [eigenvalue.real for eigenvalue in eigenvalues]
+# imaginary = [eigenvalue.imag for eigenvalue in eigenvalues]
+# plt.figure()
+# plt.scatter(real, imaginary)
+# plt.xlabel('Real')
+# plt.ylabel('Imaginary')
+# plt.tight_layout()
+# plt.savefig('Spectrum_stoch_mat_A.png')
+
+# b)
+# T = np.matrix([[0,0,0,1], [1/3,0,0,0], [1/3,1,0,0], [1/3,0,1,0]])
+# eigenvalues, eigenvectors = np.linalg.eig(T)
+# real = [eigenvalue.real for eigenvalue in eigenvalues]
+# imaginary = [eigenvalue.imag for eigenvalue in eigenvalues]
+# plt.figure()
+# plt.scatter(real, imaginary)
+# plt.xlabel('Real')
+# plt.ylabel('Imaginary')
+# plt.tight_layout()
+# plt.savefig('Spectrum_stoch_mat_B.png')
+
+# c)
+# T = np.matrix([[0,0,0,1], [1/2,0,0,0], [0,1,0,0], [1/2,0,1,0]])
+# # convergence
+# p_0_1 = np.array([1,0,1,0])
+# p_0_2 = np.array([1,0,0,0])
+# p_0_3 = np.array([0,0,1,0])
+# p_0_4 = np.array([0,1,0,0])
+# p_0_5 = np.array([0,0,0,1])
+# k = 100
+# T_k = np.linalg.matrix_power(T, k)
+# T_k_1 = np.linalg.matrix_power(T, k+1)
+
+# norm_0_1 = np.zeros(100)
+# norm_0_2 = np.zeros(100)
+# norm_0_3 = np.zeros(100)
+# norm_0_4 = np.zeros(100)
+# norm_0_5 = np.zeros(100)
+
+# for k in range(1,101):
+#     T_k = np.linalg.matrix_power(T, k)
+#     T_k_1 = np.linalg.matrix_power(T, k+1)
+#     # norm_0_1[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_1)-np.dot(T_k,p_0_1))
+#     norm_0_2[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_2)-np.dot(T_k,p_0_2))
+#     norm_0_3[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_3)-np.dot(T_k,p_0_3))
+#     norm_0_4[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_4)-np.dot(T_k,p_0_4))
+#     norm_0_5[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_5)-np.dot(T_k,p_0_5))
+
+# plt.figure()
+# # plt.plot(norm_0_1, label='1')
+# plt.plot(norm_0_2, label='2')
+# plt.plot(norm_0_3, label='3')
+# plt.plot(norm_0_4, label='4')
+# plt.plot(norm_0_5, label='5')
+# plt.legend()
+# plt.show()
+
+# # divergence
+# p_0_1 = np.array([1,1,0,0])
+# p_0_2 = np.array([1,1,1,1])
+# k = 100
+# T_k = np.linalg.matrix_power(T, k)
+# T_k_1 = np.linalg.matrix_power(T, k+1)
+
+# norm_0_1 = np.zeros(100)
+# norm_0_2 = np.zeros(100)
+# norm_0_3 = np.zeros(100)
+
+# for k in range(1,101):
+#     T_k = np.linalg.matrix_power(T, k)
+#     T_k_1 = np.linalg.matrix_power(T, k+1)
+#     norm_0_1[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_1)-np.dot(T_k,p_0_1))
+#     norm_0_2[k-1] = np.linalg.norm(np.dot(T_k_1,p_0_2)-np.dot(T_k,p_0_2))
 
+# plt.figure()
+# plt.plot(norm_0_1, label='1')
+# plt.plot(norm_0_2, label='2')
+# plt.legend()
+# plt.show()
\ No newline at end of file