SymPy Plotting Backends で2曲線の描画範囲を個別に設定して plot と塗りつぶし

In [1]:
from sympy import *
# 1文字変数の Symbol の定義が省略できる
from sympy.abc import *
# π
from sympy import pi
Pi = float(pi)

# SymPy Plotting Backends (SPB)
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

同一の描画範囲で2本の曲線を plot()

SymPy の plot() で2つの陽関数の曲線を同時にグラフにする例。

以下の例では $y=\cos x$ と $y=\sin x$ を 同一の描画範囲 $0 \leq x \leq 2\pi$ で plot() します。書式は…

plot(expr1, expr2, ..., range, **kwargs)
In [2]:
plot(cos(x), sin(x), (x, 0, 2*pi));

2本の曲線の描画範囲を個別に設定して plot()

2つの曲線の描画範囲をそれぞれ別に,また曲線の描画範囲とグラフの表示範囲を独立に設定することもできます。

以下の例では,グラフの表示範囲を横軸 (-0.2, %pi+0.2),縦軸 (-0.1, 1.1) にし,

  • $0 \leq x \leq \frac{\pi}{4}$ の範囲で $y=\cos x$
  • $0 \leq x \leq \pi$ の範囲で $y=\sin x$

plot() します。個別の描画範囲,凡例,線色の設定と,グラフの表示範囲の設定の例。書式は…

plot((expr1, range1), (expr2, range2), ..., **kwargs)
plot((expr1, range1, label1, rendering_kw1), (expr2, range2, label2, rendering_kw2), ..., **kwargs)
In [3]:
plot((cos(x), (x, 0, pi/4), "cos $x$", {"color":"red"}), 
     (sin(x), (x, 0, pi), "sin $x$", {"color":"blue"} ), 
     xlim = (0, pi), ylim = (-0.1, 1.1));

陽関数で表された2曲線の間を塗りつぶす

$0 \leq x \leq \frac{\pi}{4}$ の範囲で $y=\cos x$ と $y=\sin x$ に囲まれた領域を塗りつぶしてみます。

オプション fill に塗りつぶす $x$ の範囲 'x',領域の一方の境界を示す陽関数 'y1',同じく他方の境界を示す陽関数 'y2' の数値リストを辞書型で指定する(らしい)。数値リストは np.linspace()np.arange() を使うのが定番だが,ここでは NumPy に頼らずに試してみる。

In [4]:
# NumPy を import して使う場合
# import numpy as np
# x_array = np.linspace(0, np.pi/4, 100)
# y1_array = np.cos(x_array)
# y2_array = np.sin(x_array)

# NumPy に頼らない場合
# 0 ~ Pi/4 を 100 等分した数値リストを作成
x_list = [(Pi/4 * i/100) for i in range(101)]
# x_array に対応した y1 の数値リストを作成
y1_list = lambdify(x, cos(x))(x_list)
# x_array に対応した y2 の数値リストを作成
y2_list = lambdify(x, sin(x))(x_list)

plot((cos(x), (x, 0, pi/4), "cos $x$", {"color":"red"}), 
     (sin(x), (x, 0, pi), "sin $x$", {"color":"blue"}),
     (x, 0, pi), xlim = (0, pi), ylim = (-0.1, 1.1), 
     fill={'x': x_list,
           'y1':y1_list,
           'y2':y2_list,
           'color':'yellow'});

参考:plot_implicit() による塗りつぶし

plot_implicit() を使い,$\sin x \leq y \leq \cos x$ の領域を p1 のように塗りつぶすことができます。境界を滑らかにするために depth = 3としています。depth の値を大きく (4 が最大値) すると,それなりに滑らかになりますが,時間がかかります。また,ファイルサイズもデカくなります。

In [5]:
p1 = plot_implicit((sin(x) - y)*(y - cos(x)) > 0, (x, 0, pi/4), 
                   adaptive = True, depth = 3, 
                   xlim = (0, pi), ylim = (-0.1, 1.1));

横軸の目盛位置と目盛ラベルの設定

今度は,以下のページにある例:
$\displaystyle \frac{\pi}{4} \leq x \leq \frac{5\pi}{4}$ の範囲で $y=\sin x$ と $y=\cos x$ で囲まれた部分を塗りつぶします。

横軸の目盛を $\frac{\pi}{4}$ ごとにつける設定をしてみます。

「理工系の数学B」の積分の練習問題で,「図の塗りつぶした部分の面積を求めよ。」なんていうときに使います。

In [6]:
# Pi/4 ~ 5 Pi/4 を 100 等分した数値リストを作成
x_list = [(Pi/4 + Pi * i/100) for i in range(101)]
# x_list に対応した y1 の数値リストを作成
y1_list = lambdify(x, sin(x))(x_list)
# x_array に対応した y2 の数値リストを作成
y2_list = lambdify(x, cos(x))(x_list)

p = plot((cos(x), "cos $x$", {"color":"red"}), 
         (sin(x), "sin $x$", {"color":"blue"}), 
         (x, 0, 2*pi), xlim = (0, 2*pi), ylim = (-1.1, 1.1), 
         fill={'x': x_list,
               'y1':y1_list,
               'y2':y2_list,
               'color':'yellow'}, show = False);
ax = p.ax
# x 軸の目盛は π/4 ごとに
ax.set_xticks(
    [(i*Pi/4) for i in range(9)])
# x 軸の目盛のラベル。既約分数表記で。
ax.set_xticklabels(
    [(str(Rational(i,4))+" $\pi$") for i in range(9)]);
# 手動で設定する場合
# ax.set_xticklabels(
#    ["$0$", " $\pi/4$", "$\pi/2$", "$3\pi/4$", "$\pi$",
#     "$5\pi/4$", "$3\pi/2$", "$7\pi/4$", "$2\pi$"]);
# y 軸の目盛は 0.5 ごとに
ax.set_yticks(
    [(0.5*i) for i in range(-2,3)]);

縦線を引く

2曲線の交点である $x = \frac{\pi}{4}$ と $x = \frac{5\pi}{4}$ のところに vlines() で縦線を引いてみます。明示的に Matplotlib.pyplot を import していませんが,SPB を import すると使えるようです。

In [7]:
p = plot((cos(x), "cos $x$", {"color":"red"}), 
         (sin(x), "sin $x$", {"color":"blue"}), 
         (x, 0, 2*pi), xlim = (0, 2*pi), ylim = (-1.1, 1.1), 
         fill={'x': x_list,
               'y1':y1_list,
               'y2':y2_list,
               'color':'yellow'}, show = False);
ax = p.ax
ax.set_xticks(
    [(i*Pi/4) for i in range(9)])
# 手動で設定する場合
ax.set_xticklabels(
   ["$0$", "$\pi/4$", "$\pi/2$", "$3\pi/4$", "$\pi$",
    "$5\pi/4$", "$3\pi/2$", "$7\pi/4$", "$2\pi$"]);
ax.set_yticks(
    [(0.5*i) for i in range(-2,3)]);

# 2本の縦線を引く
ax.vlines(x=Pi/4, ymin = -1.2, ymax=sin(Pi/4), color = 'black')
ax.vlines(x=5*Pi/4, ymin = -1.2, ymax=sin(5*Pi/4), color = 'black');

異なる範囲での塗りつぶしの例

もう一つの例。$y=f(x)$ と $x$ 軸の間を,$0.5 \leq x \leq 2$ までは黄色で,$2 \leq x \leq 2.5$ までは灰色で塗りつぶしてみます。

In [8]:
import numpy as np

f = 0.6*x + 0.4*cos(x)

# 0.5 ~ 2 を 100 等分した ndarray を作成
x1_array = np.linspace(0.5, 2, 100)
# 2 ~ 2.5 を 100 等分した ndarray を作成
x2_array = np.linspace(2, 2.5, 100)
# x1_array に対応した f の ndarray を作成
f1_array = lambdify(x, f)(x1_array)
# x2_array に対応した f の ndarray を作成
f2_array = lambdify(x, f)(x2_array)

plot(f, (x, 0.2, 2.8), "",
     xlim = (-0.2, 3.5), ylim = (-0.2, 1.5), 
     fill = [{'x': x1_array, 'y1':f1_array, 'color':'yellow'}, 
             {'x': x2_array, 'y1':f2_array, 'color':'lightgray'}]);

list ではなく ndarray を使ったほうがいいかも

In [9]:
y = 0.6*x + 0.4*cos(x)
In [10]:
x_list = [0.1*i for i in range(6)]
type(x_list)
Out[10]:
list

$y = 0.6 x + 0.4 \cos x$ という簡単な関数であっても,list から lambdify()x_list に対応した y の list を作成しようとすると,以下のようになぜかエラーで作成できない。

In [11]:
# x_list に対応した y1 の array を作成しようとするとエラー
y1_list = lambdify(x, y)(x_list)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-8b6ff2fe8f18> in <module>
      1 # x_list に対応した y1 の array を作成しようとするとエラー
----> 2 y1_list = lambdify(x, y)(x_list)

<lambdifygenerated-19> in _lambdifygenerated(x)
      1 def _lambdifygenerated(x):
----> 2     return 0.6*x + 0.4*cos(x)

TypeError: can't multiply sequence by non-int of type 'float'

一方,ndarray の x_array から lambdify()y の ndarray を作成すると,特に問題なく作成できる。

In [12]:
x_array = np.linspace(0, 0.5, 6)
type(x_array)
Out[12]:
numpy.ndarray
In [13]:
# x_array に対応した y の array を作成
y1_array = lambdify(x, y)(x_array)
# エラーなし,問題なく lambdify できる。

グラフを pdf ファイルに保存

日本語を含んでいても,特に問題なく pdf ファイルとして保存できる。書式は…

plot(f(x),...).save("filename.pdf");

または,

p = plot(f(x),...);
p.save("filename.pdf");
In [14]:
p= plot(f, (x, 0.2, 2.8), "関数 $f(x)$", title = "グラフ",
     xlim = (-0.2, 3.5), ylim = (-0.2, 1.5), 
     fill = [{'x': x1_array, 'y1':f1_array, 'color':'yellow'}, 
             {'x': x2_array, 'y1':f2_array, 'color':'lightgray'}]);
p.save("spbnuri.pdf");