Pythonでグラフ計算機｜Pythonで数学を学ぼう！第18回

はじめに｜レポート

私が毎月開催している桜井進の算数・数学教室（Zoomオンラインセミナー）ですが、2/20（日）桜井進のPython・UNIX・Math教室のテーマは「Pythonで自分だけのグラフ計算機をつくる」でした。

今回の目標は、陰関数の2Dおよび3Dのグラフを出力することです。

陰関数

関数と言えば、1次関数y=3x、2次関数y=x²+x+1、三角関数y=sin x、指数関数y=e^xのようなyがxの関数であるy=f(x)の形式がよく知られています。

円x²+y²=1、双曲線x²-y²=1、だ円3x²+2y²+xy=1のようなyがxの式では表すことができない f(x, y)=0 の形をした関数が陰関数（implicit function）です。

通常の関数 y=f(x) も f(x)-y=0 と変形することで陰関数としてあつかうことができます。

座標系と座標

座標系（互いに直交する座標軸の取り方。点の位置の決め方）に対して座標（点の位置を表す数の組）が定まります。

よく使われる座標系・座標が次の2つです。

直交座標系

座標を考案したとされるデカルト（1596-1650）にちなんでデカルト座標系（Cartesian coordinate system）とも呼ばれます。

y=f(x)、f(x, y)=0、f(x, y, z) の座標 (x, y, z) が直交座標です。

極座標系（polar coordinate system）

平面 (r,θ）（円座標ともいう）、空間 (r,θ,φ) （球面座標ともいう）
で表されるのが極座標（極形式）です。原点からの距離rと2つの角θとφによって点の位置を表す方法です。

どちらも平面（空間）の点を表す表し方の違いなので、互いに変換ができます。

プログラム「2dgraph.py」

>> # 2D平面グラフ（陰関数）
>> # 2dgraph.py
>> # python -m pip install japanize-matplotlib
>> import numpy as np
>> from numpy import sin, cos, tan
>> import matplotlib.pyplot as plt
>> import japanize_matplotlib
>>
>> print('f(x, y)=0 のf(x, y)を入力します')
>> print('例 y=x^2 -> f(x, y) = x**2 - y (= 0)を入力')
>> print('例 x^2+y^2=1 -> f(x, y) = x**2 + y**2 -1 (= 0)を入力')
>> print('例 sin(xy)=1 -> f(x, y) = sin(x*y) -1 (= 0)を入力')
>> f = input('f(x, y) >>> ')
>> a = int(input('-a≦x≦a, -a≦y≦a, a >>> '))
>>
>>
>> fig = plt.figure()
>> axis = fig.add_subplot(1,1,1)
>> x, y = np.mgrid[-a: a: 4001j, -a: a: 4001j]
>>
>> axis.set_title('陰関数 (implicit funtion) '+f+'=0')
>> axis.set_xlabel(r'$x$', fontsize=12)
>> axis.set_ylabel(r'$y$', fontsize=12)
>> axis.tick_params(axis='both', length=10, which='major')
>> axis.tick_params(axis='both', length=5, which='minor')
>> axis.minorticks_on()
>> axis.set_aspect('equal', 'box')
>> axis.set(xlim=(-a, a), ylim=(-a, a))
>>
>> axis.contour(x, y, eval(f),[0], colors="blue")
>>
>> plt.grid()
>> plt.show()
>> fig.savefig('2dplot.png', transparent=False)

このプログラムファイルは次からダウンロードできます。
https://drive.google.com/file/d/1aFIIXJDD5JqJK-WdJejVpmde2dfuqCr2/view?usp=sharing

準備｜ライブラリインストール

プログラム「2dgraph.py」にはjapanize-matplotlibライブラリが必要になります。これはグラフ描画ライブラリmatplotlibを日本語表示させるためのものです。

matplotlibをimportした後、japanize_matplotlibをimportします。

ターミナルで次を入力します。このコマンド文はプログラム「2dgraph.py」の冒頭に記載してあります。

「2dgraph.py」をダブルクリックしてテキストエディタで開き、該当テキストをコピーします。

ターミナル上でペーストして実行すればOK。
$ python -m pip install japanize-matplotlib

実行結果

① 2次関数y=x²

ターミナルに次を入力します。
$ python 2dgraph.py

f(x, y)とaの入力を求められます。
f(x, y)=0のf(x, y)にはx²-y（またはy-x²）を入力します。

次にaに対して2を入力することで、変数xとyの範囲を-2≦x≦2、-2≦y≦2と設定します。

グラフが別ウィンドウに表示され、画像ファイル 2dplot.png も同時に出力されます。（「2dgraph.py」があるディレクトリ）

②楕円曲線 y² = x³-4x

③x²+y²+sin(6x)+sin(6y)=1

プログラム「3dgraph.py」

>> # 3Dグラフ（陰関数）
>> # implicit3d3.py
>> # ライブラリscikit-imageインストール
>> # python -m pip install scikit-image
>> import numpy as np
>> from numpy import sin, cos, tan
>> from skimage import measure
>> import matplotlib.pyplot as plt
>> import japanize_matplotlib
>> from mpl_toolkits.mplot3d import axes3d
>>
>> print('f(x, y, z)=0 のf(x, y, z)を入力します')
>> print('例 x^2+y^2+z^2=1 -> f(x, y, z) = x**2 + y**2 + z**2 -1 (= 0)を入力')
>> print('例 z = cos(x*y) -> f(x, y, z) = cos(x*y) - z (= 0)を入力')
>> print('例 x^2+y^2+z^2+sin(4x)+sin(4y)+sin(4z)=1 -> f(x, y, z)=x**2+y**2+z**2+sin(4*x)+sin(4*y)+sin(4*z)-1 (= 0)を入力')
>> f = input('f(x, y, z) >>> ')
>>
>> def implicit_function(x,y,z):
>> global f
>> return eval(f)
>>
>> xl = np.linspace(-2, 2, 100)
>> X, Y, Z = np.meshgrid(xl, xl, xl)
>> F = implicit_function(X, Y, Z)
>>
>> verts, faces, normals, values = measure.marching_cubes(F, 0, spacing=[np.diff(xl)[0]]*3)
>> verts -= 3
>>
>> fig = plt.figure()
>> axis = fig.add_subplot(111, projection='3d')
>> axis.set_title('陰関数 (implicit funtion) '+f+'=0')
>> axis.set_xlabel(r'$x$', fontsize=12)
>> axis.set_ylabel(r'$y$', fontsize=12)
>> axis.set_zlabel(r'$z$', fontsize=12)
>>
>> axis.plot_trisurf(verts[:, 0], verts[:, 1], faces, verts[:, 2], cmap='hsv', lw=0, alpha=0.7,shade=True, antialiased=True) # cmap='hsv' カラーマップ
>> # カラーマップ https://matplotlib.org/stable/tutorials/colors/colormaps.html
>>
>> plt.show()
>> fig.savefig('3dplot.png', transparent=False)