42
38

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

呉高専Advent Calendar 2019

Day 12

Matplotlibを使って高校数学+大学数学に出てくる関数をビジュアル化

Last updated at Posted at 2019-12-11

この記事は呉高専 Advent Calendar 2019 12日目の記事です。

更新(2021/06/19): 「点の集まり(散布図)」を追加

はじめに

数学の勉強をしていると、さまざまな「関数」に出会います。我々はそれらの関数をしばしば微分したり積分したり、あるいはその関数の極値を求めたり変曲点を求めたりします。しかし、これらの操作は機械的な処理で計算できるため、もとの関数がどのようなどのような性質をもつか、どのような形であるかはあまり意識する必要がありません。それはそれで数学の偉大なところですが、今回は高校数学や大学数学で扱うような関数のグラフをきちんとイメージできるように、Pythonを使って様々な形で記述されるグラフを描画していきます。

本記事は 「コードをコピペしてあとはちょろっと書き換えるだけ」 を心がけて作成しています。そのため記事中にはコメントアウト多めのコードをたくさん記載しています。記事が長くなっているため、適宜 ジャンプ機能 を活用してお読みください。

環境

Python3系

  • 各自環境を用意してください

numpy

  • pip3 install numpyによりインストールできます

matplotlib

  • pip3 install matplotlibによりインストールできます

japanize-matplotlib

  • pip3 install japanize-matplotlib によりインストールできます
  • matplotlibで日本語を扱えるようにします

2Dグラフ

y = f(x)

以下は $y = x^2$ のグラフを描画する例です。

fig.png

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# 定義域を指定
x = np.linspace(-10, 10, 1000000)

# ここに関数を記述
y = x**2

plt.plot(x, y)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# 自動で軸をとりたい場合
plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合
# plt.axes().set_aspect('equal')
# plt.xlim([-10, 10])
# plt.ylim([-100, 100])

plt.show()

定義域(xの範囲)の指定

x = np.linspace(xの最小値, xの最大値, 分割数 = 1000000)により指定します。

分割数は基本1000000固定で良いですが、処理が重い場合は値を小さくしてください。

関数の記述

y = x**2の部分を適宜自分の描画したい関数に書き換えます。初等関数を利用する際はnpにある関数を利用します。mathモジュールは使わないようにしましょう。

例えば

  • $y = \sin{x}$ の場合はy = np.sin(x)

  • $y = e^x$ の場合はy = np.exp(x)

とします。

軸の設定

自動で軸を取る場合はそのままで良いですが、自分で軸の最大値最小値を設定したくなる場合もあります。

その場合は、

# 自動で軸をとりたい場合
# plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合
plt.axes().set_aspect('equal')
plt.xlim([xの最小値, xの最大値])
plt.ylim([yの最小値, yの最大値])

としましょう。

2つ以上のグラフを同時に描画

$y = \sin{x}$ と $y = \cos{x}$ を同時に描画する例です。

image.png

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# 定義域を指定
x = np.linspace(0, 2*PI, 1000000)

# ここに関数を記述
y1 = np.sin(x)
y2 = np.cos(x)

plt.plot(x, y1, label='sin')
plt.plot(x, y2, label='cos')

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# ラベルの描画
plt.legend()

# 自動で軸をとりたい場合
plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合
# plt.axes().set_aspect('equal')
# plt.xlim([-10, 10])
# plt.ylim([0, 100])

plt.show()

2つ以上の関数を同時に描画する場合は、label='ラベル名'をつけるようにしましょう。

3つ以上の場合でも同様に関数を増やしていけば良いです。グラフの色も自動で設定してくれます。

対数軸で描画

以下は $y = e^{-x} = \exp{(-x)}$ のグラフを対数軸で描画する例です。
電気の学生なら対数軸を扱うことも多くあると思います。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# 定義域の設定
x = np.linspace(-5, 5, 1000000)

# ここに関数を記述
y = np.exp(-x)

plt.plot(x, y, label='exp(-x)')

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# ラベルの描画
plt.legend()

# 対数軸に設定
plt.yscale('log')  # x軸も対数にする場合はplt.xscale('log')を追加

# 自動で軸をとりたい場合
# plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合

# plt.axes().set_aspect('equal')
# plt.xlim([-10, 10])
# plt.ylim([0, 100])

plt.show()

今回は対数軸なので「自動で軸をとる」設定は外しています。

陰関数の描画

陰関数とは $F(x, y) = 0$ の形で記述される関数です。

例えば、$x^2 + y^2 - 1 = 0$は陰関数と言えます。 $y = f(x)$のように $y$ を陽に記述できない場合によく使われますね。

以下は、$F(x, y) = x^2 - 2xy - y^2 + 7 = 0$ を描画する例です。どんなグラフになるか想像できるでしょうか?

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# x の範囲
x_range = np.linspace(-10, 10, 1000)
# y の範囲
y_range = np.linspace(-10, 10, 1000)

x, y = np.meshgrid(x_range, y_range)

# 陰関数 F = 0
F = x**2 - 2*x*y - y**2 + 7

# 軸を自動で設定する場合
plt.gca().set_aspect('equal', adjustable='box')
# 軸を手動で設定する場合
# plt.axis([-10, 10, -10, 10])

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.contour(x, y, F, [0], colors="blue")

plt.show()

細かい説明はしませんが、2変数関数 $z = F(x, y)$ の高さ $z = 0$ の地点の等高線を記述するという方法で陰関数を描画しています。

実際書き換えるのは関数 $F(x, y)$ の部分と $x$ と $y$ の範囲くらいでしょう。

曲線の媒介変数表示

媒介変数表示とは、関数を $y=f(x)$の形で表現せず、媒介変数と呼ばれる変数 $t$ によって表現する方法です。式で書くと以下のようになります。

\left\{
\begin{array}{ll}
x = \varphi(t) \\
y = \psi(t)
\end{array}
\right.

以下は

\left\{
\begin{array}{ll}
x = \left(1 + \cos{t} \right) \cos{t} \\
y = \left(1 + \cos{t} \right) \sin{t}
\end{array}
\right.

のグラフを描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# t の範囲
t = np.linspace(0, 2*PI, 1000000)

# 関数を記述
x = (1 + np.cos(t))*np.cos(t)
y = (1 + np.cos(t))*np.sin(t)

# 軸を自動で設定する場合
plt.gca().set_aspect('equal', adjustable='box')
# 軸を手動で設定する場合
# plt.axis([-10, 10, -10, 10])

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.plot(x, y)

plt.show()

例はいわゆるカージオイドです。

媒介変数表示を使えば興味深いグラフが簡単に得られますね。
参考 : 媒介変数表示された有名な曲線7つ - 高校数学の美しい物語

2次元のベクトル関数

高校数学までで扱うベクトルは向きと大きさが常に一定な「定ベクトル」でした。
一方、大学数学からは時間などの変数によってベクトルの大きさや向きが変わる「ベクトル関数」を扱います。

ベクトル関数は $\vec{r}(t) = (x(t),\ y(t))$ という形で書かれます。

以下は $\vec{r}(t) = (\cos^3{t},\ \sin^3{t})$ のグラフ(軌跡)を描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# t の範囲
t = np.linspace(0, 2*PI, 1000000)

# x(t) を記述
x = (np.cos(t))**3
# y(t) を記述
y = (np.sin(t))**3

# 軸を自動で設定する場合
plt.gca().set_aspect('equal', adjustable='box')
# 軸を手動で設定する場合
# plt.axes().set_aspect('equal')
# plt.axis([-10, 10, -10, 10])

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.plot(x, y)

plt.show()

プログラムは「曲線の媒介変数表示」とほぼ同じ、というか全く同じになりますね。

ベクトル軌跡(ナイキスト線図)

複素数 $z(t) = x(t) + y(t)\cdot i$ について、$t$を動かしたときの複素数の軌跡をプロットしたものを「ベクトル線図・ナイキスト線図」などといいます。

電気回路理論を例を出せば、インピーダンス $\dot{Z} = R + j \omega L$について、抵抗 $R$ が変化したらインピーダンス $\dot{Z}$ はどうなるかなーとか角周波数 $\omega$ を変化させたらどなるかなーとかっていうのを見るのに使います。

特に制御工学ではナイキスト線図という言い方が好まれます。伝達関数の周波数応答を見るために使います。(これ以上解説を詳しくすると長くなるのでやめておきます。)

以下は

z(\omega) = \dfrac{1}{1-\omega^2 + \omega i}

について、 $\omega$ を変化させる場合のベクトル軌跡を描画する例です。プログラム中では $\omega$ を t で表現しています。

image.png

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import japanise_matplotlib
PI = np.pi # 円周率をPIで使えるようにする
 
# 計算式
t = np.arange(0, 10, 0.001) # np.arange(tの最小値(普通 0 ), tの最大値, 刻み)

# 関数を記述
## 複素数は 4i などと `j` をつければ良い
## `変数*i` などという場合は `変数 * 1j` とする
z = 1 / (1 - t**2 + t*1j);

# 実部と虚部を取得
Re = z.real
Im = z.imag

# 横軸の変数。縦軸の変数。
plt.plot(Re, Im)

# グリッド補助線
plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# 軸を自動で設定する場合
plt.gca().set_aspect('equal', adjustable='box')
# 軸を手動で設定する場合
# plt.axes().set_aspect('equal')
# plt.axis([-10, 10, -10, 10])

# 描画実行
plt.show()

Pythonでは虚数単位は $i$ ではなく $j$ を用います(電気を勉強している人ならおなじみですね)。

例えば、

  • $4i$を表現したいときは、 4j

  • $2ti$を表現したいときは 2*t * 1j

とします。

有理化などをせず、そのまま記述すればプロットされるので嬉しいです。

二次元ベクトル

matplotlibにはベクトルを描画する関数 plt.quiver() がありますが、定ベクトルを書くには少し構文が厄介です。そこで plot_vector2D() という関数を自作し描画を簡単にしました。

以下は

\begin{align}
\vec{v_1} &= ( 5,\ 1) \\
\vec{v_2} &= ( 2,\ 4) \\
\vec{v_3} &= \vec{v_2} - \vec{v_1}
\end{align}

の3つのベクトルを描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

def plot_vector2D(pos, vec, color='blue', scale=1.0):
    plt.quiver(pos[0], pos[1], vec[0], vec[1], color=color, angles='xy', scale_units='xy', scale=scale)

# ベクトルを記述
v1 = np.array([5, 1])
v2 = np.array([2, 4])
v3 = v2 - v1

plot_vector2D([0, 0], v1, color='red')
plot_vector2D([0, 0], v2, color='green')
plot_vector2D(v1, v3)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# 軸の設定
## x の最大値最小値
plt.xlim([-6, 6])
## y の最大値最小値
plt.ylim([-6, 6])

plt.axes().set_aspect('equal')

plt.show()

plot_vector2D関数

定ベクトルをプロットするための関数 plot_vector2D() に以下のように引数を渡せばベクトルがプロットされます。

plot_vector2D(始点, 方向, color='色(省略可)', scale=単位ベクトルのスケール省略可)

例えば原点からベクトルを伸ばす場合は plot_vector2D([0, 0], v) とします。
colorscaleは、それぞれデフォルト値のblue1.0になります。

scale で指定した大きさですべてのベクトルが割られて描画されます。
例えば scale=5.0 と指定すると、すべてのベクトルの大きさが5分の1になります。

scale を変更しても始点の位置は調整されないことに注意してください。

基本的にはデフォルト値のscale=1.0 で良いかと思いますが、後述のベクトル場では適宜値を変更したほうが良いかもしれません。

軸の設定

plt.axes().set_aspect('equal', 'datalim') を使ってもうまく自動で調整されないことが多いので、手動で設定する必要があります。

plt.xlim([x軸の最小値, x軸の最大値])
plt.ylim([y軸の最小値, y軸の最大値])

のように設定してください。

二次元のベクトル場

水の流れのように、各場所・各点 $P(x,\ y)$ でベクトル $\boldsymbol{v}(P) = \boldsymbol{v}(x,\ y)$ が分布しているとき、それをベクトル場といいます。

以下は ベクトル場

\boldsymbol{v}(x,\ y) = \left(2y,\ 3x \right)

を描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
PI = np.pi  # 円周率をPIで使えるようにする

# x の範囲
x_range = np.arange(-3, 3, 0.5)
# y の範囲
y_range = np.arange(-3, 3, 0.5)
x, y = np.meshgrid(x_range, y_range)

# ベクトル場 v(x, y) = (u(x, y), v(x, y))
u = 2*y
v = 3*x

plt.quiver(x, y, u, v, color='blue', angles='xy', scale_units='xy', scale=5.0)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# 自動で軸をとりたい場合
plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合
# plt.axes().set_aspect('equal')
# plt.xlim([-10, 10])
# plt.ylim([0, 100])

plt.show()

x, y の範囲

いままで $x,\ y$ の範囲を指定する関数として np.linspace(最小値, 最大値, 分割数) を使っていましたが、ここでは np.arange(最小値, 最大値, 間隔) を使っています。 分割数でデータの細かさを表現するより、矢印をどのくらいの間隔で置くかで表現したほうがわかりやすいからです。

関数の記述

ベクトル場を $\boldsymbol{v}(x, y) = \left(u(x,\ y),\ v(x,\ y) \right)$ と表したときの $u(x,\ y),\ v(x,\ y)$ を書けばよいです。

変数 xy はベクトルの始点 $(x,\ y)$を表しています。

plt.quiver関数

plt.quiver関数でベクトルをプロットしていますが、

plt.quiver(x, y, u, v, color='blue', angles='xy', scale_units='xy', scale=5.0)

の引数のうち、scale=5.0 には注意してください。
本来のスケールで描画するなら scale=1.0 を指定すればよいです。
しかし、スケールが1だと矢印が長すぎて見えづらくなることもよくあります。そこで scale で数値を指定すれば、その大きさでベクトルが割られ、見やすくなります。

点の集まり(散布図)

点の集まり $x, y$ をそのまま plt.plot してしまうと、折れ線グラフが描画されてしまいます。

点の集まり(散布図)を描画するときは plt.scatter を使用します。

以下は、sin波っぽい(?)点の集まり書く例です。

scat.png

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# 点を記述
x = [-3.1416, -2.8274, -2.5133, -1.2566, 0.62832, 2.5133, 3.1416]
y = [-0.36717, -0.10182, -0.10289, -1.1583, 0.47613, 0.24801, -0.02356]

plt.scatter(x, y)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# 自動で軸をとりたい場合
plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合
# plt.axes().set_aspect('equal')
# plt.xlim([-10, 10])
# plt.ylim([-100, 100])

plt.show()

なお、$x, y$ を手打ちではなく、外部CSVファイルなどを読み込んで代入したい場合は、以下のようにすれば良いです。

data = np.loadtxt('sin_data.csv', delimiter=',')
x = data[:, 0]
y = data[:, 1]

3Dグラフ

z = f(x, y)

二変数関数のグラフは三次元空間に描画されます。

以下は $z = y^2 - x^2$ のグラフを描画する例です。

image.png

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import japanize_matplotlib
PI = np.pi # 円周率をPIで使えるようにする

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# x と y の範囲の設定
x_range = np.linspace(-20, 20, 1000)
y_range = np.linspace(-20, 20, 1000)
x, y = np.meshgrid(x_range, y_range)

# 関数を記述
z = y**2 - x**2

ax.plot_surface(x, y, z, cmap = "summer")
ax.contour(x, y, z, colors = "gray", offset = -1)  # 底面に等高線を描画

# 自動で軸を設定する場合は記述なし

# 手動で軸を設定する場合
# ax.set_xlim(-20, 20)
# ax.set_ylim(-20, 20)
# ax.set_zlim(-1500, 1500)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show() 

コードの雰囲気が少し変わりましたが、いじるところは2Dの場合と大差ありません。書き換えるとすれば $z = f(x, y)$ の記述と、軸の設定くらいでしょう。

軸の設定(スケール)には十分注意しましょう。今回は見栄え重視のため軸を自動で設定しました。そのため、x,y軸とz軸のスケールが異なっています(2変数関数はxとyのスケールに比べてzのスケールは大きくなりがちなので、手動でスケールを合わせたところで微妙な感じにはなりますが…)

とにかくグラフを眺めるときは軸のスケールの確認を忘れないようにしましょう。

三次元のベクトル関数

以下は $\vec{r}(t) = (\sin{t},\ \cos{t},\ t)$ のグラフを描画する例です。

image.png

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import japanize_matplotlib
PI = np.pi # 円周率をPIで使えるようにする

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# パラメータ t の範囲を設定
t = np.linspace(-6*PI, 6*PI, 1000)

# x(t) を記述
x = np.cos(t)
# y(t) を記述
y = np.sin(t)
# z(t) を記述
z = t

ax.plot(x, y, z)

# 自動で軸を設定する場合は記述なし

# 手動で軸を設定する場合
# ax.set_xlim(-20, 20)
# ax.set_ylim(-20, 20)
# ax.set_zlim(-1500, 1500)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show() 

物理の話をすれば、速度 $\vec{v}(t)$ や 位置 $\vec{x}(t)$ などが三次元ベクトル関数の例にあげられます。

三次元ベクトル

ベクトルを描画する関数 plt.quiver() を使いますが、二次元同様、定ベクトルを書くには少し構文が厄介です。
そこで plot_vector3D() という関数を自作し描画を簡単にしました。

以下は

\vec{v_1} = (1,\ 2,\ 3) \\
\vec{v_2} = (-4,\ 5,\ 1) \\
\vec{v_3} = \vec{v_2} - \vec{v_1}

の3つのベクトルを描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

def plot_vector3D(ax, pos, vec, color='blue', scale=1.0):
    ax.quiver(pos[0], pos[1], pos[2], vec[0], vec[1], vec[2], color=color, length=1/scale)

# ベクトルを記述
v1 = np.array([1, 2, 3])
v2 = np.array([-4, 4, 1])
v3 = v2 - v1

plot_vector3D(ax, [0, 0, 0], v1, color='red')
plot_vector3D(ax, [0, 0, 0], v2, color='green')
plot_vector3D(ax, v1, v3)

# 軸の設定
ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
ax.set_zlim([-5, 5])

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

plot_vector3D()関数に以下のように引数を渡せばベクトルがプロットされます。

plot_vector3D(ax, 始点, 方向, color='色(省略可)', scale=単位ベクトルのスケール)

scale を変更しても始点の位置は調整されないことに注意してください。

2次元の場合とほぼ同じですね。

三次元のベクトル場

三次元のベクトル場も二次元と同様に定義されます。

以下は ベクトル場

\boldsymbol{v}(t) = \left( x(x+y+z),\ y(x+y+z),\ z(x+y+z) \right)

を描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# x の範囲
x_range = np.arange(-2, 2, 1)
# y の範囲
y_range = np.arange(-2, 2, 1)
# z の範囲
z_range = np.arange(-2, 2, 1)
x, y, z = np.meshgrid(x_range, y_range, z_range)

# ベクトル場 v(x, y, z) = (u(x, y, z), v(x, y, z), w(x, y, z))
u = x*(x + y + z)
v = y*(x + y + z)
w = z*(x + y + z)

ax.quiver(x, y, z, u, v, w, color='blue', length=1/5.0, normalize=False)

ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
ax.set_zlim([-5, 5])

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

二次元の場合とほぼ同じですね。

単位ベクトルのスケールの指定がscaleではなくlengthになっていることに注意してください。
単位ベクトルが length 倍のスケールで描画されます。

uv曲面

3次元空間内の曲面は、位置ベクトル $\boldsymbol{r}$ が2変数 $u,\ v$の関数で、$\boldsymbol{r} = \boldsymbol{r}(u,\ v)$ で表すことができます。これを$uv$曲面と呼びます。

ベクトル解析の分野でよく使われますね。

以下は $\boldsymbol{r}(u,v) = (2u,\ 3v,\ u^2 + v^2)$ のグラフ(曲面)を描画する例です。

image.png

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import japanize_matplotlib
PI = np.pi # 円周率をPIで使えるようにする

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

# u の範囲を指定
urange = np.linspace(-20, 20, 1000)
# v の範囲を指定
vrange = np.linspace(-20, 20, 1000)
u, v = np.meshgrid(urange, vrange)

# 曲面 r(u, v) を記述
x = 2*u
y = 3*v
z = u**2 + v**2

ax.plot_surface(x, y, z, cmap = "summer")
ax.contour(x, y, z, colors = "gray", offset = -1)  # 底面に等高線を描画

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show() 

応用

グラフを装飾する

グラフを装飾する方法はたくさんありますが、ここではその一部を紹介します。

以下は オームの法則 $V = RI$ を描画した例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# 定義域を指定
I = np.linspace(-5, 5, 1000000)

# ここに関数を記述
R = 1.0
V = R*I

plt.plot(I, V)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

# タイトルをつける
plt.title(r'オームの法則(抵抗 $R = 1$ [${\rm \Omega}$])')

plt.xlabel('電流 $I$ [A]')
plt.ylabel('電圧 $V$ [V]')

# 自動で軸をとりたい場合
plt.axes().set_aspect('equal', 'datalim')
# 手動で軸を取りたい場合
# plt.axes().set_aspect('equal')
# plt.xlim([-10, 10])
# plt.ylim([-100, 100])

# 目盛りを 1 ごとつける
plt.xticks(np.arange(-5, 5+1, 1))
plt.yticks(np.arange(-5, 5+1, 1))

# x軸とy軸を実線にする
plt.hlines([0], -5, 5, color='gray')
plt.vlines([0], -5*R, 5*R, color='gray')

plt.show()

タイトルやラベルをつける

plt.title('')plt.xlabel('')plt.ylabel('') を用いることでグラフにタイトルやラベルをつけることができます。

$ で囲うことで、$\TeX$記法を使うこともできますが、バックスラッシュ\を使う場合は文字列を囲うクォーテーションの前に r をつける必要があります。

目盛りの細かさを変更する

plt.xticks(np.arange(-5, 5+1, 1))plt.yticks(np.arange(-5, 5+1, 1)) を使ってxとyの目盛りの細かさをそれぞれ変更できます。

np.arange(最小値, 最大値, 間隔(公差)) ですが、範囲が [最小値, 最大値) であることには注意しましょう。(最大値は含まないということです。)

x軸とy軸を実線にする

水平線を引く関数 plt.hlines() と 垂直線を引く関数 plt.vlines() を使ってx軸とy軸を引いてみました。

補助線を描画する

オームの法則を使った例ではx軸とy軸を引くために plt.hlines()plt.vlines() を使いましたが、これらの関数は補助線を描画するのに便利です。

$y = \sin{x}$ について

\begin{align}
x &= \pi,\ 2\pi \\
y &= -1,\ -\frac{\sqrt{2}}{2},\ \frac{\sqrt{2}}{2},\ 1
\end{align}

の部分に補助線を描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

x = np.linspace(0, 2*PI + PI/2, 1000000)

y = np.sin(x)

plt.plot(x, y)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.title(r'$y = \sin{x}$')

plt.axes().set_aspect('equal')
plt.xlim([0, 2*PI + PI/2])
plt.ylim([-2, 2])

# 水平な補助線を引く
plt.hlines([-1, -np.sqrt(2)/2, np.sqrt(2)/2, 1], 0, 2*PI + PI/2, color='red', linestyle='--')
# 垂直な補助線を引く
plt.vlines([PI, 2*PI], -2, 2, color='red', linestyle='--')

plt.hlines([0], 0, 2*PI + PI/2, color='gray')
plt.vlines([0], -2, 2, color='gray')

plt.show()

関数の引数の意味は以下の通りです。

plt.hlines([補助線を引き始めるyの値のリスト], xの最小値, xの最大値, color='', linestyle='-や--や:')
plt.vlines([補助線を引き始めるxの値のリスト], yの最小値, yの最大値, color='', linestyle='-や--や:')

2つ以上のグラフを並べて描画

2つ以上のグラフを並べて描画するには、fig.add_subplot()関数を使います。

以下は 曲面 $z = x^2 + y^2 - 1$ と 曲線 $F(x, y) = x^2 + y^2 - 1 = 0$ を同時に描画する例です。
image.png

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

x_range = np.linspace(-2, 2, 1000)
y_range = np.linspace(-2, 2, 1000)
x, y = np.meshgrid(x_range, y_range)

# 陰関数
z = x**2 + y**2 - 1

# ウィンドウの大きさを指定(グラフ1枚あたり 8 を想定)
fig = plt.figure(figsize=(16, 8))

# 描画の設定
ax1 = fig.add_subplot(1, 2, 1, projection="3d")
ax1.set_title('曲面z = f(x, y)')
ax1.set_xlabel("x", size = 16)
ax1.set_ylabel("y", size = 16)
ax1.set_zlabel("z", size = 16)
ax1.plot_surface(x, y, z, cmap = "summer")
ax1.contour(x, y, z, colors = "black", offset = -1)  # 底面に等高線を描画
ax1.set_aspect('equal', adjustable='datalim')

ax2 = fig.add_subplot(1, 2, 2)
ax2.set_title('曲線 f(x, y) = 0')
ax2.axis([-2, 2, -2, 2])
ax2.set_aspect('equal', adjustable='datalim')
ax2.grid(which='major', color='gray', linestyle='--')
ax2.grid(which='minor', color='gray', linestyle='--')
ax2.contour(x, y, z, [0], colors="blue")

plt.show()

2枚以上のグラフを扱う場合は、ax という変数を使ってそれぞれのグラフの設定・プロットを行うと考えてください。

まず、 fig.add_subplot()関数を使って、ひとつの ax を作成します。
「何あるグラフのうち、何番目にそのグラフを置くか」というように引数を指定します。

# 2Dの場合
ax = fig.add_subplot(, , 何番目)
# 3Dの場合
ax = fig.add_subplot(, 列数 何番目, projection="3d")

次に ax を使って設定、プロットを行います。 ax にある関数は plt にある関数とかなり似ていますが、関数名が全く同じというわけでもないので各自調べて下さい。

そして、この ax を表示したいグラフの枚数だけ作成します。

特別な関数

周期関数や条件分岐がある関数の描画を考えます。
結局は np.ndarray 型をいかに扱うかなのでわからないことがあれば np.ndarray をキーワードに検索をすれば良いでしょう。

いくつか例をあげます。

ReLU関数(ランプ関数)

ReLU関数(ランプ関数)は以下のように定義されます。

f(x) = \left\{
\begin{array}{ll}
x & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.

numpynp.maximum(a, b) を使うことでうまく関数を記述することができます。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

x = np.linspace(-5, 5, 1000000)

# ReLU関数
y = np.maximum(x, 0)

plt.plot(x, y)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.title('ReLU関数')

plt.axes().set_aspect('equal', 'datalim')

plt.xticks(np.arange(-5, 5+1, 1))
plt.yticks(np.arange(-5, 5+1, 1))

plt.hlines([0], -5, 5, color='gray')
plt.vlines([0], -5, 5, color='gray')

plt.show()

条件分岐がある関数

f(x) = \left\{
\begin{array}{ll}
\dfrac{3}{2}x & (0 \leq x < 2) \\
3 & (2 \leq x < 4) \\
-\dfrac{3}{2}x + 9 & (otherwise)
\end{array}
\right.

のグラフを描画することを考えましょう。

np.ndarray のなかで条件に合う要素を抜き出すのが np.where() 関数です。

np.where(条件式, Trueの場合の値, Falseの場合の値)

という使い方をします。複数の条件式を組み合わせる場合は andor の代わりに &| を使います。

以下に実装例を示します。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

x = np.linspace(0, 6, 1000000)

# 関数の定義
y = np.where((0 <= x) & (x < 2), 1.5*x, np.where((2 <= x) & (x < 4), 3, -1.5*x + 9))

plt.plot(x, y)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.axes().set_aspect('equal', 'datalim')

plt.xticks(np.arange(0, 6+1, 1))
plt.yticks(np.arange(0, 6+1, 1))

plt.hlines([0], 0, 6, color='gray')
plt.vlines([0], min(y)-1, max(y)+1, color='gray')

plt.show()

関数の実装部分がかなり読みにくいですが、これで実現できます。

こんな関数も面白いかもしれませんね~

# 関数の定義
y = np.sin(x)
y = np.where(np.abs(y) > np.sqrt(3) / 2, np.sign(y)*np.sqrt(3) / 2, y)

周期関数

周期の剰余(MOD)を計算するというアイデアで周期関数を記述することができます。

以下は $f(x) = x^2,\ \ f(x+2) = f(x)$ のグラフを描画する例です。

image.png

import math
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

x = np.linspace(-6, 6, 1000000)

# 周期 T の関数
T = 2
y = (x % T)**2


plt.plot(x, y)

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.axes().set_aspect('equal', 'datalim')

plt.xticks(np.arange(-6, 6+1, 1))
plt.yticks(np.arange(-1, 5+1, 1))

plt.hlines([0], -6, 6, color='gray')
plt.vlines([0], -1, 5, color='gray')

plt.show()

接線の方程式を描画する

ある点 $x = a$ における関数 $f(x)$ の接線の方程式を描画しましょう。

まず、接線の傾き $f'(a)$ を求めます。定義式から $f'(a)$ は

f'(a) = \lim_{h \to 0} \frac{f(a+h) - f(a)}{h}

で求めることができます。

さて、人間が計算するのであればこの定義式からゴチャゴチャ式を変形して……としますが、コンピュータ的には $h = 10^{-10}$ というように $h$ に微小な値を代入することで近似的に微分係数を計算することができます。ただし、もとの定義式のままだと、「 $a$ 点から $a+h$ 点までの差 / $h$」なので、正確に $a$ 点で微分することができません。そこで、

f'(a) = \lim_{h \to 0} \frac{f(a+h) - f(a-h)}{2h}

とすることで、(ある程度)正確に $a$ 点での微分を求めます。

次に、接線の方程式を求めます。通る点と直線の公式からずばり、

y = f'(a) (x - a) + f(a)

です。さあ、あとは実装するだけですね。

以下は $y = x^2 + 2x + 1$ のグラフと、 $x = 1$ における接線のグラフを描画する例です。

image.png

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
PI = np.pi  # 円周率をPIで使えるようにする

# 対象の関数
def f(x):
    return 0.5*x**2 + x + 0.5

# 微分係数を求める
def diff(fun, a):
    h = 10**(-10)
    return (fun(a + h) - fun(a - h)) / (2*h)

# 関数の定義
## 対象の関数
x1 = np.linspace(-6, 4, 1000000)
y1 = f(x1)

## 接線
x2 = np.linspace(-1, 5, 1000000)
a = 1
k = diff(f, a)
y2 = k*(x2 - a) + f(a)

plt.plot(x1, y1, label='$y=x^2 + 2x + 1$')
plt.plot(x2, y2, label='接線')

plt.grid(which='major', color='gray', linestyle='--')
plt.grid(which='minor', color='gray', linestyle='--')

plt.legend()

plt.axes().set_aspect('equal', 'datalim')

plt.show()

最後に

Matplotlibの関数はここで紹介した以外にもたくさんあります。ぜひいろいろ試してみてください。

気分によってはもうすこしサンプルを追加するかも。

42
38
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
42
38

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?