質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.35%

  • Python

    9138questions

    Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

  • Python 3.x

    7345questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

  • Matplotlib

    370questions

    MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

python matplotlib 複数のコードのグラフを一つにまとめたい

受付中

回答 2

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 1,501

閲覧ありがとうございます。python3でmatplotlibを用いてグラフを描いたのですが、このグラフを一つにまとめたいです。
具体的には二つのコードがありまして
①L=8のコード(E8、Cv8、M8、chi8という4つのグラフが出てきます)

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
@author: zhshang
"""
import matplotlib
matplotlib.use('TKAgg')
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

L = 8
ESTEP = 1000
STEP = 10000

J = 1 # J>0 to make it ferromagnetic

# Intitialize the XY network
def Init():
    return np.random.rand(L, L)*2*np.pi
    #return np.ones([L, L])

# periodic neighbor
def next(x):
    if x == L-1:
        return 0
    else:
        return x+1

# construct the bond lattice
def FreezeBonds(Ising,T,S):
    iBondFrozen = np.zeros([L,L],dtype=int)
    jBondFrozen = np.zeros([L,L],dtype=int)
    for i in np.arange(L):
        for j in np.arange(L):
            freezProb_nexti = 1 - np.exp(-2 * J * S[i][j] * S[next(i)][j] / T)
            freezProb_nextj = 1 - np.exp(-2 * J * S[i][j] * S[i][next(j)] / T)
            if (Ising[i][j] == Ising[next(i)][j]) and (np.random.rand() < freezProb_nexti):
                iBondFrozen[i][j] = 1
            if (Ising[i][j] == Ising[i][next(j)]) and (np.random.rand() < freezProb_nextj):
                jBondFrozen[i][j] = 1
    return iBondFrozen, jBondFrozen

# H-K algorithm to identify clusters
def properlabel(prp_label,i):
    while prp_label[i] != i:
        i = prp_label[i]
    return i

# Swendsen-Wang cluster
def clusterfind(iBondFrozen,jBondFrozen):
    cluster = np.zeros([L, L],dtype=int)
    prp_label = np.zeros(L**2,dtype=int)
    label = 0
    for i in np.arange(L):
        for j in np.arange(L):
            bonds = 0
            ibonds = np.zeros(4,dtype=int)
            jbonds = np.zeros(4,dtype=int)

            # check to (i-1,j)
            if (i > 0) and iBondFrozen[i-1][j]:
                ibonds[bonds] = i-1
                jbonds[bonds] = j
                bonds += 1
            # (i,j) at i edge, check to (i+1,j)
            if (i == L-1) and iBondFrozen[i][j]:
                ibonds[bonds] = 0
                jbonds[bonds] = j
                bonds += 1
            # check to (i,j-1)
            if (j > 0) and jBondFrozen[i][j-1]:
                ibonds[bonds] = i
                jbonds[bonds] = j-1
                bonds += 1
            # (i,j) at j edge, check to (i,j+1)
            if (j == L-1) and jBondFrozen[i][j]:
                ibonds[bonds] = i
                jbonds[bonds] = 0
                bonds += 1

            # check and label clusters
            if bonds == 0:
                cluster[i][j] = label
                prp_label[label] = label
                label += 1
            else:
                minlabel = label
                for b in np.arange(bonds):
                    plabel = properlabel(prp_label,cluster[ibonds[b]][jbonds[b]])
                    if minlabel > plabel:
                        minlabel = plabel

                cluster[i][j] = minlabel
                # link to the previous labels
                for b in np.arange(bonds):
                    plabel_n = cluster[ibonds[b]][jbonds[b]]
                    prp_label[plabel_n] = minlabel
                    # re-set the labels on connected sites
                    cluster[ibonds[b]][jbonds[b]] = minlabel
    return cluster, prp_label

# flip the cluster spins
def flipCluster(Ising,cluster,prp_label):
    for i in np.arange(L):
        for j in np.arange(L):
            # relabel all the cluster labels with the right ones
            cluster[i][j] = properlabel(prp_label,cluster[i][j])
    sNewChosen = np.zeros(L**2)
    sNew = np.zeros(L**2)
    flips = 0 # get the number of flipped spins to calculate the Endiff and Magdiff
    for i in np.arange(L):
        for j in np.arange(L):
            label = cluster[i][j]
            randn = np.random.rand()
            # mark the flipped label, use this label to flip all the cluster elements with this label
            if (not sNewChosen[label]) and randn < 0.5:
                sNew[label] = +1
                sNewChosen[label] = True
            elif (not sNewChosen[label]) and randn >= 0.5:
                sNew[label] = -1
                sNewChosen[label] = True

            if Ising[i][j] != sNew[label]:
                Ising[i][j] = sNew[label]
                flips += 1

    return Ising,flips

# Swendsen-Wang Algorithm in Ising model (with coupling constant dependency on sites)
# One-step for Ising
def oneMCstepIsing(Ising, S):
    [iBondFrozen, jBondFrozen] = FreezeBonds(Ising, T, S)
    [SWcluster, prp_label] = clusterfind(iBondFrozen, jBondFrozen)
    [Ising, flips] = flipCluster(Ising, SWcluster, prp_label)
    return Ising

# Decompose XY network to two Ising networks with project direction proj
def decompose(XY,proj):
    x = np.cos(XY)
    y = np.sin(XY)
    x_rot = np.multiply(x,np.cos(proj))+np.multiply(y,np.sin(proj))
    y_rot = -np.multiply(x,np.sin(proj))+np.multiply(y,np.cos(proj))
    Isingx = np.sign(x_rot)
    Isingy = np.sign(y_rot)
    S_x = np.absolute(x_rot)
    S_y = np.absolute(y_rot)
    return Isingx, Isingy, S_x, S_y

# Compose two Ising networks to XY network
def compose(Isingx_new,Isingy_new,proj,S_x, S_y):
    x_rot_new = np.multiply(Isingx_new,S_x)
    y_rot_new = np.multiply(Isingy_new,S_y)
    x_new = np.multiply(x_rot_new,np.cos(proj))-np.multiply(y_rot_new,np.sin(proj))
    y_new = np.multiply(x_rot_new,np.sin(proj))+np.multiply(y_rot_new,np.cos(proj))
    XY_new = np.arctan2(y_new,x_new)
    return XY_new

def oneMCstepXY(XY):
    proj = np.random.rand()
    [Isingx, Isingy, S_x, S_y] = decompose(XY, proj)
    Isingx_new = oneMCstepIsing(Isingx, S_x)
    #Isingy_new = oneMCstepIsing(Isingy, S_y)
    # Here we only use the flopping of Ising in projection direction, without the perpendicular direction
    #XY_new = compose(Isingx_new, Isingy_new, proj, S_x, S_y)
    XY_new = compose(Isingx_new, Isingy, proj, S_x, S_y)
    return XY_new

# Calculate the energy for XY network
def EnMag(XY):
    energy = 0
    for i in np.arange(L):
        for j in np.arange(L):
            # energy
            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]))
    magx = np.sum(np.cos(XY))
    magy = np.sum(np.sin(XY))
    mag = np.array([magx,magy])
    return energy * 0.5, LA.norm(mag)/(L**2)

# Swendsen Wang method for XY model
def SWang(T):
    XY = Init()
    # thermal steps to get the equilibrium
    for step in np.arange(ESTEP):
        XY = oneMCstepXY(XY)
    # finish with thermal equilibrium, and begin to calc observables
    E_sum = 0
    M_sum = 0
    Esq_sum = 0
    Msq_sum = 0
    for step in np.arange(STEP):
        XY = oneMCstepXY(XY)
        [E,M] = EnMag(XY)

        E_sum += E
        M_sum += M
        Esq_sum += E**2
        Msq_sum += M**2

    E_mean = E_sum/STEP/(L**2)
    M_mean = M_sum/STEP
    Esq_mean = Esq_sum/STEP/(L**4)
    Msq_mean = Msq_sum/STEP

    return XY, E_mean, M_mean, Esq_mean, Msq_mean

M = np.array([])
E = np.array([])
M_sus = np.array([])
SpcH = np.array([])
Trange = np.linspace(0.1, 2.5, 10)
for T in Trange:
    [Ising, E_mean, M_mean, Esq_mean, Msq_mean] = SWang(T)
    M = np.append(M, np.abs(M_mean))
    E = np.append(E, E_mean)
    M_sus = np.append(M_sus, 1/T*(Msq_mean-M_mean**2))
    SpcH = np.append(SpcH, 1/T**2*(Esq_mean-E_mean**2))

# plot the figures
T = Trange

plt.figure()
plt.plot(T, E, 'rx-')
plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
plt.ylabel(r'$\langle E \rangle$ per site $(J)$')
plt.savefig("E8.pdf", format='pdf', bbox_inches='tight')

plt.figure()
plt.plot(T, SpcH, 'kx-')
plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
plt.ylabel(r'$C_V$ per site $(\frac{J^2}{k_B^2})$')
plt.savefig("Cv8.pdf", format='pdf', bbox_inches='tight')

plt.figure()
plt.plot(T, M, 'bx-')
plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
plt.ylabel(r'$\langle|M|\rangle$ per site $(\mu)$')
plt.savefig("M8.pdf", format='pdf', bbox_inches='tight')

plt.figure()
plt.plot(T, M_sus, 'gx-')
plt.xlabel(r'Temperature $(\frac{J}{k_B})$')
plt.ylabel(r'$\chi$ $(\frac{\mu}{k_B})$')
plt.savefig("chi8.pdf", format='pdf', bbox_inches='tight')

plt.tight_layout()
fig = plt.gcf()
plt.show()

np.savetxt('output.data',np.c_[T,E,SpcH,M,M_sus])

②L=16のコード(コード省略、①のコードにおいてはじめのL=8をL=16に変更するだけです

で出てくる計8つのグラフを各種類(E,M,Cv,chi)ごとに4つのグラフにまとめたいです
例えば、以下の写真のようにしたいと考えています。
Cv(L=8,16,32,64の4つのグラフを合体させたパターン)
詳しい方、よろしくお願い致します。

追記

コード①とコード②からはそれぞれ以下のようなグラフが出てきます(Cv8とCv16のみ貼っています)
イメージ説明イメージ説明
このように別々のコードから出てきたグラフを1つのグラフにまとめたいと考えています。(縦軸の値を合わせて・・・)

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 2

0

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/12 10:42

    回答ありがとうございます。サイトの方を確認させていただきましたが、少し私の考えていた内容とは異なっていました。コード①とコード②で出る別々の結果(グラフ)を一つのグラフにまとめて表示させたいです。後ほど補足情報を質問欄に足しておきます。

    キャンセル

0

plt.subplotsでaxisを4つ用意して、そのplotメソッドを使って上書きしていけばよさそうに思います。

# あらかじめ2x2の4つの軸オブジェクトを用意する。
fig, axs = plt.subplots(2, 2)
symbols = ["co", "rx", "ys", "m^"]
labels = [8, 16, 32, 64]

# L = 8の結果
step = 0
xx = [T] * 4
yy = [E, SpcH, M, M_sus]

symbol, label = symbols[step], labels[step]
for x, y, ax in zip(xx, yy, axs.ravel()):  # axsは2x2になっているので、ravel()で1次元にする
    ax.plot(x, y, symbol, label=label)

# L = 16の結果
step = 1
xx = [T] * 4
yy = [E/2, SpcH/2, M/2, M_sus/2]    

symbol, label = symbols[step], labels[step]
for x, y, ax in zip(xx, yy, axs.ravel()):  # axsは2x2になっているので、ravel()で1次元にする
    ax.plot(x, y, symbol, label=label)

# legendを付け足す
for ax in axs.ravel():
    ax.legend()

plt.show()

結果

こんな感じ

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/05/22 11:23

    tachikome様 回答ありがとうございます。コードの方を確認させていただいたのですが、yy = [E/2, SpcH/2, M/2, M_sus/2] となっているところが気になりました。L=8とL=16は系のサイズ(L)が異なるだけで同じアルゴリズムなのですが、決して÷2のような関係にあるわけではありません。分布の形自体は同じなのですが、数値がL=8とL=16とで一対一に対応しないのです・・・。。欲しいグラフは確かに回答者様のように”4つのパラメーターにそれぞれ2つずつ折れ線がある”もので間違いないです。

    キャンセル

  • 2018/05/22 11:24 編集

    それは異なるプロットを作るためのタミーデータです。適当に結果の数値に置き換えて使ってください。

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.35%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

同じタグがついた質問を見る

  • Python

    9138questions

    Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

  • Python 3.x

    7345questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

  • Matplotlib

    370questions

    MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。