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()

こんな結果になります．

線がはみ出してしまいますね．．．
これは$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()

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