teratail
回答率
85.30%
質問するログイン新規登録
トップに関する質問python matplotlib グラフが表示されない

編集履歴

質問編集履歴

4

書式の改善を行いました

2018/05/11 11:51

投稿

astromelt0416
astromelt0416

スコア15

title CHANGED
File without changes
body CHANGED
@@ -66,15 +66,19 @@
66
66
 
67
67
  J = 1 # J>0 to make it ferromagnetic
68
68
 
69
+ # Intitialize the XY network
69
70
  def Init():
70
71
  return np.random.rand(L, L)*2*np.pi
71
-
72
+ #return np.ones([L, L])
73
+
74
+ # periodic neighbor
72
75
  def next(x):
73
76
  if x == L-1:
74
77
  return 0
75
78
  else:
76
79
  return x+1
77
80
 
81
+ # construct the bond lattice
78
82
  def FreezeBonds(Ising,T,S):
79
83
  iBondFrozen = np.zeros([L,L],dtype=int)
80
84
  jBondFrozen = np.zeros([L,L],dtype=int)
@@ -88,11 +92,13 @@
88
92
  jBondFrozen[i][j] = 1
89
93
  return iBondFrozen, jBondFrozen
90
94
 
95
+ # H-K algorithm to identify clusters
91
96
  def properlabel(prp_label,i):
92
97
  while prp_label[i] != i:
93
98
  i = prp_label[i]
94
99
  return i
95
100
 
101
+ # Swendsen-Wang cluster
96
102
  def clusterfind(iBondFrozen,jBondFrozen):
97
103
  cluster = np.zeros([L, L],dtype=int)
98
104
  prp_label = np.zeros(L**2,dtype=int)
@@ -103,27 +109,28 @@
103
109
  ibonds = np.zeros(4,dtype=int)
104
110
  jbonds = np.zeros(4,dtype=int)
105
111
 
112
+ # check to (i-1,j)
106
113
  if (i > 0) and iBondFrozen[i-1][j]:
107
114
  ibonds[bonds] = i-1
108
115
  jbonds[bonds] = j
109
116
  bonds += 1
110
-
117
+ # (i,j) at i edge, check to (i+1,j)
111
118
  if (i == L-1) and iBondFrozen[i][j]:
112
119
  ibonds[bonds] = 0
113
120
  jbonds[bonds] = j
114
121
  bonds += 1
115
-
122
+ # check to (i,j-1)
116
123
  if (j > 0) and jBondFrozen[i][j-1]:
117
124
  ibonds[bonds] = i
118
125
  jbonds[bonds] = j-1
119
126
  bonds += 1
120
-
127
+ # (i,j) at j edge, check to (i,j+1)
121
128
  if (j == L-1) and jBondFrozen[i][j]:
122
129
  ibonds[bonds] = i
123
130
  jbonds[bonds] = 0
124
131
  bonds += 1
125
132
 
126
-
133
+ # check and label clusters
127
134
  if bonds == 0:
128
135
  cluster[i][j] = label
129
136
  prp_label[label] = label
@@ -136,19 +143,19 @@
136
143
  minlabel = plabel
137
144
 
138
145
  cluster[i][j] = minlabel
139
-
146
+ # link to the previous labels
140
147
  for b in np.arange(bonds):
141
148
  plabel_n = cluster[ibonds[b]][jbonds[b]]
142
149
  prp_label[plabel_n] = minlabel
143
-
150
+ # re-set the labels on connected sites
144
151
  cluster[ibonds[b]][jbonds[b]] = minlabel
145
152
  return cluster, prp_label
146
153
 
147
-
154
+ # flip the cluster spins
148
155
  def flipCluster(Ising,cluster,prp_label):
149
156
  for i in np.arange(L):
150
157
  for j in np.arange(L):
151
-
158
+ # relabel all the cluster labels with the right ones
152
159
  cluster[i][j] = properlabel(prp_label,cluster[i][j])
153
160
  sNewChosen = np.zeros(L**2)
154
161
  sNew = np.zeros(L**2)
@@ -157,7 +164,7 @@
157
164
  for j in np.arange(L):
158
165
  label = cluster[i][j]
159
166
  randn = np.random.rand()
160
-
167
+ # mark the flipped label, use this label to flip all the cluster elements with this label
161
168
  if (not sNewChosen[label]) and randn < 0.5:
162
169
  sNew[label] = +1
163
170
  sNewChosen[label] = True
@@ -171,14 +178,15 @@
171
178
 
172
179
  return Ising,flips
173
180
 
174
-
181
+ # Swendsen-Wang Algorithm in Ising model (with coupling constant dependency on sites)
182
+ # One-step for Ising
175
183
  def oneMCstepIsing(Ising, S):
176
184
  [iBondFrozen, jBondFrozen] = FreezeBonds(Ising, T, S)
177
185
  [SWcluster, prp_label] = clusterfind(iBondFrozen, jBondFrozen)
178
186
  [Ising, flips] = flipCluster(Ising, SWcluster, prp_label)
179
187
  return Ising
180
188
 
181
-
189
+ # Decompose XY network to two Ising networks with project direction proj
182
190
  def decompose(XY,proj):
183
191
  x = np.cos(XY)
184
192
  y = np.sin(XY)
@@ -190,7 +198,7 @@
190
198
  S_y = np.absolute(y_rot)
191
199
  return Isingx, Isingy, S_x, S_y
192
200
 
193
-
201
+ # Compose two Ising networks to XY network
194
202
  def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
195
203
  x_rot_new = np.multiply(Isingx_new,S_x)
196
204
  y_rot_new = np.multiply(Isingy_new,S_y)
@@ -203,29 +211,31 @@
203
211
  proj = np.random.rand()
204
212
  [Isingx, Isingy, S_x, S_y] = decompose(XY, proj)
205
213
  Isingx_new = oneMCstepIsing(Isingx, S_x)
206
-
214
+ #Isingy_new = oneMCstepIsing(Isingy, S_y)
215
+ # Here we only use the flopping of Ising in projection direction, without the perpendicular direction
216
+ #XY_new = compose(Isingx_new, Isingy_new, proj, S_x, S_y)
207
217
  XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
208
218
  return XY_new
209
219
 
210
-
220
+ # Calculate the energy for XY network
211
221
  def EnMag(XY):
212
222
  energy = 0
213
223
  for i in np.arange(L):
214
224
  for j in np.arange(L):
215
-
225
+ # energy
216
226
  energy = energy - (np.cos(XY[i][j]-XY[(i-1)%L][j])+np.cos(XY[i][j]-XY[(i+1)%L][j])+np.cos(XY[i][j]-XY[i][(j-1)%L])+np.cos(XY[i][j]-XY[i][(j+1)%L]))
217
227
  magx = np.sum(np.cos(XY))
218
228
  magy = np.sum(np.sin(XY))
219
229
  mag = np.array([magx,magy])
220
230
  return energy * 0.5, LA.norm(mag)/(L**2)
221
231
 
222
-
232
+ # Swendsen Wang method for XY model
223
233
  def SWang(T):
224
234
  XY = Init()
225
-
235
+ # thermal steps to get the equilibrium
226
236
  for step in np.arange(ESTEP):
227
237
  XY = oneMCstepXY(XY)
228
-
238
+ # finish with thermal equilibrium, and begin to calc observables
229
239
  E_sum = 0
230
240
  M_sum = 0
231
241
  Esq_sum = 0
@@ -258,6 +268,7 @@
258
268
  M_sus = np.append(M_sus, 1/T*(Msq_mean-M_mean**2))
259
269
  SpcH = np.append(SpcH, 1/T**2*(Esq_mean-E_mean**2))
260
270
 
271
+ # plot the figures
261
272
  T = Trange
262
273
 
263
274
  plt.figure()
@@ -290,4 +301,5 @@
290
301
 
291
302
  np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
292
303
 
304
+
293
305
  ```

3

書式を改善しました

2018/05/11 11:51

投稿

astromelt0416
astromelt0416

スコア15

title CHANGED
File without changes
body CHANGED
@@ -51,6 +51,8 @@
51
51
 
52
52
  ### 補足情報（FW/ツールのバージョンなど）
53
53
 
54
+ ```
55
+
54
56
  ```python
55
57
  import matplotlib
56
58
  matplotlib.use('Agg')
@@ -287,7 +289,5 @@
287
289
  plt.show()
288
290
 
289
291
  np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
290
- ```
291
292
 
292
-
293
293
  ```

2

コードの全文を貼りました

2018/05/11 11:47

投稿

astromelt0416
astromelt0416

スコア15

title CHANGED
File without changes
body CHANGED
@@ -51,4 +51,243 @@
51
51
 
52
52
  ### 補足情報（FW/ツールのバージョンなど）
53
53
 
54
+ ```python
55
+ import matplotlib
56
+ matplotlib.use('Agg')
57
+ import numpy as np
58
+ from numpy import linalg as LA
59
+ import matplotlib.pyplot as plt
60
+
61
+ L = 8
62
+ ESTEP = 1000
63
+ STEP = 10000
64
+
65
+ J = 1 # J>0 to make it ferromagnetic
66
+
67
+ def Init():
68
+ return np.random.rand(L, L)*2*np.pi
69
+
70
+ def next(x):
71
+ if x == L-1:
72
+ return 0
73
+ else:
74
+ return x+1
75
+
76
+ def FreezeBonds(Ising,T,S):
77
+ iBondFrozen = np.zeros([L,L],dtype=int)
78
+ jBondFrozen = np.zeros([L,L],dtype=int)
79
+ for i in np.arange(L):
80
+ for j in np.arange(L):
81
+ freezProb_nexti = 1 - np.exp(-2 * J * S[i][j] * S[next(i)][j] / T)
82
+ freezProb_nextj = 1 - np.exp(-2 * J * S[i][j] * S[i][next(j)] / T)
83
+ if (Ising[i][j] == Ising[next(i)][j]) and (np.random.rand() < freezProb_nexti):
84
+ iBondFrozen[i][j] = 1
85
+ if (Ising[i][j] == Ising[i][next(j)]) and (np.random.rand() < freezProb_nextj):
86
+ jBondFrozen[i][j] = 1
87
+ return iBondFrozen, jBondFrozen
88
+
89
+ def properlabel(prp_label,i):
90
+ while prp_label[i] != i:
91
+ i = prp_label[i]
92
+ return i
93
+
94
+ def clusterfind(iBondFrozen,jBondFrozen):
95
+ cluster = np.zeros([L, L],dtype=int)
96
+ prp_label = np.zeros(L**2,dtype=int)
97
+ label = 0
98
+ for i in np.arange(L):
99
+ for j in np.arange(L):
100
+ bonds = 0
101
+ ibonds = np.zeros(4,dtype=int)
102
+ jbonds = np.zeros(4,dtype=int)
103
+
104
+ if (i > 0) and iBondFrozen[i-1][j]:
105
+ ibonds[bonds] = i-1
106
+ jbonds[bonds] = j
107
+ bonds += 1
108
+
109
+ if (i == L-1) and iBondFrozen[i][j]:
110
+ ibonds[bonds] = 0
111
+ jbonds[bonds] = j
112
+ bonds += 1
113
+
114
+ if (j > 0) and jBondFrozen[i][j-1]:
115
+ ibonds[bonds] = i
116
+ jbonds[bonds] = j-1
117
+ bonds += 1
118
+
119
+ if (j == L-1) and jBondFrozen[i][j]:
120
+ ibonds[bonds] = i
121
+ jbonds[bonds] = 0
122
+ bonds += 1
123
+
124
+
125
+ if bonds == 0:
126
+ cluster[i][j] = label
127
+ prp_label[label] = label
128
+ label += 1
129
+ else:
130
+ minlabel = label
131
+ for b in np.arange(bonds):
132
+ plabel = properlabel(prp_label,cluster[ibonds[b]][jbonds[b]])
133
+ if minlabel > plabel:
134
+ minlabel = plabel
135
+
136
+ cluster[i][j] = minlabel
137
+
138
+ for b in np.arange(bonds):
139
+ plabel_n = cluster[ibonds[b]][jbonds[b]]
140
+ prp_label[plabel_n] = minlabel
141
+
142
+ cluster[ibonds[b]][jbonds[b]] = minlabel
143
+ return cluster, prp_label
144
+
145
+
146
+ def flipCluster(Ising,cluster,prp_label):
147
+ for i in np.arange(L):
148
+ for j in np.arange(L):
149
+
150
+ cluster[i][j] = properlabel(prp_label,cluster[i][j])
151
+ sNewChosen = np.zeros(L**2)
152
+ sNew = np.zeros(L**2)
153
+ flips = 0 # get the number of flipped spins to calculate the Endiff and Magdiff
154
+ for i in np.arange(L):
155
+ for j in np.arange(L):
156
+ label = cluster[i][j]
157
+ randn = np.random.rand()
158
+
159
+ if (not sNewChosen[label]) and randn < 0.5:
160
+ sNew[label] = +1
161
+ sNewChosen[label] = True
162
+ elif (not sNewChosen[label]) and randn >= 0.5:
163
+ sNew[label] = -1
164
+ sNewChosen[label] = True
165
+
166
+ if Ising[i][j] != sNew[label]:
167
+ Ising[i][j] = sNew[label]
168
+ flips += 1
169
+
170
+ return Ising,flips
171
+
172
+
173
+ def oneMCstepIsing(Ising, S):
174
+ [iBondFrozen, jBondFrozen] = FreezeBonds(Ising, T, S)
175
+ [SWcluster, prp_label] = clusterfind(iBondFrozen, jBondFrozen)
176
+ [Ising, flips] = flipCluster(Ising, SWcluster, prp_label)
177
+ return Ising
178
+
179
+
180
+ def decompose(XY,proj):
181
+ x = np.cos(XY)
182
+ y = np.sin(XY)
183
+ x_rot = np.multiply(x,np.cos(proj))+np.multiply(y,np.sin(proj))
184
+ y_rot = -np.multiply(x,np.sin(proj))+np.multiply(y,np.cos(proj))
185
+ Isingx = np.sign(x_rot)
186
+ Isingy = np.sign(y_rot)
187
+ S_x = np.absolute(x_rot)
188
+ S_y = np.absolute(y_rot)
189
+ return Isingx, Isingy, S_x, S_y
190
+
191
+
192
+ def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
193
+ x_rot_new = np.multiply(Isingx_new,S_x)
194
+ y_rot_new = np.multiply(Isingy_new,S_y)
195
+ x_new = np.multiply(x_rot_new,np.cos(proj))-np.multiply(y_rot_new,np.sin(proj))
196
+ y_new = np.multiply(x_rot_new,np.sin(proj))+np.multiply(y_rot_new,np.cos(proj))
197
+ XY_new = np.arctan2(y_new,x_new)
198
+ return XY_new
199
+
54
- ここにより詳細な情報を記載してください。
200
+ def oneMCstepXY(XY):
201
+ proj = np.random.rand()
202
+ [Isingx, Isingy, S_x, S_y] = decompose(XY, proj)
203
+ Isingx_new = oneMCstepIsing(Isingx, S_x)
204
+
205
+ XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
206
+ return XY_new
207
+
208
+
209
+ def EnMag(XY):
210
+ energy = 0
211
+ for i in np.arange(L):
212
+ for j in np.arange(L):
213
+
214
+ energy = energy - (np.cos(XY[i][j]-XY[(i-1)%L][j])+np.cos(XY[i][j]-XY[(i+1)%L][j])+np.cos(XY[i][j]-XY[i][(j-1)%L])+np.cos(XY[i][j]-XY[i][(j+1)%L]))
215
+ magx = np.sum(np.cos(XY))
216
+ magy = np.sum(np.sin(XY))
217
+ mag = np.array([magx,magy])
218
+ return energy * 0.5, LA.norm(mag)/(L**2)
219
+
220
+
221
+ def SWang(T):
222
+ XY = Init()
223
+
224
+ for step in np.arange(ESTEP):
225
+ XY = oneMCstepXY(XY)
226
+
227
+ E_sum = 0
228
+ M_sum = 0
229
+ Esq_sum = 0
230
+ Msq_sum = 0
231
+ for step in np.arange(STEP):
232
+ XY = oneMCstepXY(XY)
233
+ [E,M] = EnMag(XY)
234
+
235
+ E_sum += E
236
+ M_sum += M
237
+ Esq_sum += E**2
238
+ Msq_sum += M**2
239
+
240
+ E_mean = E_sum/STEP/(L**2)
241
+ M_mean = M_sum/STEP
242
+ Esq_mean = Esq_sum/STEP/(L**4)
243
+ Msq_mean = Msq_sum/STEP
244
+
245
+ return XY, E_mean, M_mean, Esq_mean, Msq_mean
246
+
247
+ M = np.array([])fig = plt.gcf()
248
+ E = np.array([])
249
+ M_sus = np.array([])
250
+ SpcH = np.array([])
251
+ Trange = np.linspace(0.1, 2.5, 10)
252
+ for T in Trange:
253
+ [Ising, E_mean, M_mean, Esq_mean, Msq_mean] = SWang(T)
254
+ M = np.append(M, np.abs(M_mean))
255
+ E = np.append(E, E_mean)
256
+ M_sus = np.append(M_sus, 1/T*(Msq_mean-M_mean**2))
257
+ SpcH = np.append(SpcH, 1/T**2*(Esq_mean-E_mean**2))
258
+
259
+ T = Trange
260
+
261
+ plt.figure()
262
+ plt.plot(T, E, 'rx-')
263
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
264
+ plt.ylabel(r'$\langle E \rangle$ per site $(J)$')
265
+ plt.savefig("E8.pdf", format='pdf', bbox_inches='tight')
266
+
267
+ plt.figure()
268
+ plt.plot(T, SpcH, 'kx-')
269
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
270
+ plt.ylabel(r'$C_V$ per site $(\frac{J^2}{k_B^2})$')
271
+ plt.savefig("Cv8.pdf", format='pdf', bbox_inches='tight')
272
+
273
+ plt.figure()
274
+ plt.plot(T, M, 'bx-')
275
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
276
+ plt.ylabel(r'$\langle|M|\rangle$ per site $(\mu)$')
277
+ plt.savefig("M8.pdf", format='pdf', bbox_inches='tight')
278
+
279
+ plt.figure()
280
+ plt.plot(T, M_sus, 'gx-')
281
+ plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
282
+ plt.ylabel(r'$\chi$ $(\frac{\mu}{k_B})$')
283
+ plt.savefig("chi8.pdf", format='pdf', bbox_inches='tight')
284
+
285
+ plt.tight_layout()
286
+ fig = plt.gcf()
287
+ plt.show()
288
+
289
+ np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])
290
+ ```
291
+
292
+
293
+ ```

1

初心者マークを入れ忘れていました。python初心者です。

2018/05/11 11:45

投稿

astromelt0416
astromelt0416

スコア15

title CHANGED
File without changes
body CHANGED
File without changes