Help us understand the problem. What is going on with this article?

matplotlibで三角グラフを作成する 2

この記事を読むとできるようになること
matplotlibで三角グラフに3変数関数の解が描ける

matplotlibで三角グラフを作成する の続編です.
仕事柄,三角グラフに点をプロットするよりも3変数関数の縮退している解(三角グラフ上で直線になる)を描くことが多いので,その描き方のまとめです.

  • 環境
    • macOS Mojave 10.14.4
    • Python 3.7.2

例として,三角グラフ上に次の3変数関数の解を描くとします.

ax+by+cz = E \quad\cdots (1)\\
x+y+z = 1 \qquad\cdots (2)\\
0 \leq x,y,z \leq 1

3変数で式が2本なので,解は縮退します.
点$(x,y,z)$を,2次元直交座標$(X,Y)$に変換して三角グラフ上に描画することにします.
変換の式は以下です.

X = z +\frac{x}{2}\\
Y = \frac{\sqrt{3}}{2}x

まず,$x$をx = np.arange(0,1.0,0.01)などで用意しておきます.
(1)(2)より$y$を消去すると,

z =  ((E-b) -(a-b)x)/(c-b)

となります.
これを使って$(X,Y)$座標に直します.

triangle2.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import math

x = np.arange(0,1.0,0.01)

fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.set_aspect('equal', 'datalim')
plt.tick_params(labelbottom=False, labelleft=False, labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

h = np.sqrt(3.0)*0.5

#ax+by+cz = E の係数
a = 2
b = 3
c = 6
E = 4 

#####formule for triangle plot######
z = ((E-b) -(a-b)*x)/(c-b)
X = z +x*0.5
Y = h*x

ax1.plot(X,Y, color='red')

################# triangle plot settings #######################
#内側目盛
for i in range(1,10):
    ax1.plot([i/20.0, 1.0-i/20.0],[h*i/10.0, h*i/10.0], color='gray', lw=0.5)
    ax1.plot([i/20.0, i/10.0],[h*i/10.0, 0.0], color='gray', lw=0.5)
    ax1.plot([0.5+i/20.0, i/10.0],[h*(1.0-i/10.0), 0.0], color='gray', lw=0.5)

#外周
ax1.plot([0.0, 1.0],[0.0, 0.0], 'k-', lw=2)
ax1.plot([0.0, 0.5],[0.0, h], 'k-', lw=2)
ax1.plot([1.0, 0.5],[0.0, h], 'k-', lw=2)

#頂点のラベル
ax1.text(0.52, h, 'x', fontsize=16)
ax1.text(-0.05, 0.03, 'y', fontsize=16, rotation=300)
ax1.text(0.95, -0.05, 'z', fontsize=16, rotation=60)

#軸ラベル
for i in range(1,10):
    ax1.text(0.5+(10-i)/20.0, h*(1.0-(10-i)/10.0), '%d0' % i, fontsize=10)
    ax1.text((10-i)/20.0-0.04, h*(10-i)/10.0+0.02, '%d0' % i, fontsize=10, rotation=300)
    ax1.text(i/10.0-0.03, -0.025, '%d0' % i, fontsize=10, rotation=60)
#################################################################

plt.tight_layout()
plt.show()

こんな結果になります.
2-1.jpeg

線がはみ出してしまいますね...
これは$y$が負になるときの解も描いてしまっているのが原因です.

なんとかしたいので,無理やりはみ出た分を消してみました.

triangle2_2.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import math
from pylab import *

#xの範囲を適当に狭めてみる
x = np.arange(0,0.7,0.01)

fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.set_aspect('equal', 'datalim')
plt.tick_params(labelbottom=False, labelleft=False, labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

h = np.sqrt(3.0)*0.5

#ax+by+cz = E の係数
a = 2
b = 3
c = 6
E = 4 

#####formule for triangle plot######
z = ((E-b) -(a-b)*x)/(c-b)
X = z +x*0.5
Y = h*x

#はみ出た部分を外枠で隠すために,むりやり折り返してみる
y1 = np.sqrt(3.0)*X
y2 = np.sqrt(3.0)*(1.0-X)
y = minimum(y1, y2)
Y0 = minimum(Y, y)
ax1.plot(X,Y0, color='red')

################# triangle plot settings #######################
#内側目盛
for i in range(1,10):
    ax1.plot([i/20.0, 1.0-i/20.0],[h*i/10.0, h*i/10.0], color='gray', lw=0.5)
    ax1.plot([i/20.0, i/10.0],[h*i/10.0, 0.0], color='gray', lw=0.5)
    ax1.plot([0.5+i/20.0, i/10.0],[h*(1.0-i/10.0), 0.0], color='gray', lw=0.5)

#外周
ax1.plot([0.0, 1.0],[0.0, 0.0], 'k-', lw=2)
ax1.plot([0.0, 0.5],[0.0, h], 'k-', lw=2)
ax1.plot([1.0, 0.5],[0.0, h], 'k-', lw=2)

#頂点のラベル
ax1.text(0.52, h, 'x', fontsize=16)
ax1.text(-0.05, 0.03, 'y', fontsize=16, rotation=300)
ax1.text(0.95, -0.05, 'z', fontsize=16, rotation=60)

#軸ラベル
for i in range(1,10):
    ax1.text(0.5+(10-i)/20.0, h*(1.0-(10-i)/10.0), '%d0' % i, fontsize=10)
    ax1.text((10-i)/20.0-0.04, h*(10-i)/10.0+0.02, '%d0' % i, fontsize=10, rotation=300)
    ax1.text(i/10.0-0.03, -0.025, '%d0' % i, fontsize=10, rotation=60)
#################################################################

plt.tight_layout()
plt.show()

2-2.jpeg

安直にはみ出た部分を隠してみました.
あまりスマートなやり方ではないと思うので,もっと賢い方法をご存知でしたらご教示ください.

kumamupooh
惑星形成の研究をしていました. 理論計算をする人に役立つ記事を中心に書いていました. SEの卵です. 最近はwebがらみのシステムを開発しています.
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away