sympyを使ってみる(1)　入力された一変数関数のグラフをプロットする

「はじめに」の前に

この記事は何本かのシリーズに分けて投稿する予定(ここ大事)です。
あまり注目されることのないsympyに焦点を当てているので、併用してる他のライブラリなどの説明は面倒くさいので省く、あるいは後述となることとします。

はじめに

　この記事はPythonにおいて、sympyを用いてユーザーの入力した関数のグラフを表示するためのプログラムについてです。
　プログラムの流れとしては

• ユーザーにグラフの範囲を設定してもらう(デフォルト値は各象限10×10)
• ユーザー側から関数の入力を受ける(このプログラムでは一変数関数f(x)のみ対応している。)
• 指定した範囲で計算をする
• グラフに描画する

という感じです。

使用したライブラリ

• matplotlib : グラフ描画のライブラリ
• sympy : 代数計算のライブラリ
• numpy : 配列を扱うライブラリ(まぁ今回は無くてもよかった)

何故sympyなのか

　Pythonでのグラフ描画といえば、最もポピュラーなのはnumpy×matplotlibです。
　何故sympyを使ってやるのかというと一応理由がありまして、

• 直感的に書きやすい
  sympyは　変数の定義→立式→代入や連立　という流れで書いていくため、数学的な意味を考えながらのコーディングが出来ます。数学の証明問題のように進んでいくので、デバッグや改良がしやすいのも特徴の一つです。
• 不等式なども扱える
  sympyは不等式などのクソめんどい煩雑且つ精神を蝕むような計算にも対応しており、グラフなどを用いた数学的な使い方をする上で大変相性がいいです。
• 方程式が解ける
  ここがかなり美味しい所なのですが、sympyは色々な形の方程式を解くことが出来ます。
  例えば二つのグラフの交点などを知りたい時でも、グラフに点を打ったうえで点の座標も表示するなど、matplotlibと組み合わせることでグラフと数値の両方の面から考えることが(比較的短いコードで)出来ます。

実行環境

• Windows10(64bit)
• Anaconda(4.5.11)
  • Python(3.7.0)
  • matplotlib(3.0.2)
  • sympy(1.3)
  • numpy(1.15.4)
• Atom(1.33.1)
  • Hydrogen(2.7.0)

※個人的な嗜好により、Jupyter notebookではなくAtom上でJupyterを動かすプラグインHydrogenを利用しました。

ソースコード(Python3.7)

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp

Pi = sp.S.Pi # 円周率
E = sp.S.Exp1 # 自然対数の底
I = sp.S.ImaginaryUnit # 虚数単位

def reset(axis_data): #ユーザーでグラフの表示範囲を設定する
    print("setting:reset")
    for n in range(len(axis_data)):
        axis_data[n] = int(input(axis_data[n]+">>"))

def usedefult(axis_data): #デフォルト値(-10~10)を使う
    print("setting:defult")
    axis_data[0:4] = -10,10,-10,10

def setting(): #範囲指定をするかの選択
    print("グラフの設定をしますか？")
    axis_ = ["xmin","xmax","ymin","ymax"]

    if input("[yes/no]") == "yes":
        reset(axis_)
    else : #yes/noすら打てない人たまにいるので
        usedefult(axis_)

    global x_min,x_max
    x_min = axis_[0]
    x_max = axis_[1]+0.1

    plt.axis(axis_)
    plt.axhline(0,c='black',lw=2)
    plt.axvline(0,c='black',lw=2)
    plt.plot(0,0,'o',c='black')
    plt.grid()

def plot(step=0.1): #グラフをplotする
    x = sp.Symbol('x')
    f = sp.sympify(input("f>>"))
    y = []
    x_ = np.arange(x_min,x_max,step)
    for i in x_:
        f_ = f.subs(x,i)
        y.append(f_)
    plt.plot(x_,y,label=f)
    plt.legend()
    plt.show()

try:
    setting()
    plot()
except:
    print("ERROR")

コードの説明

それではコードを追っていきましょう。

前半：setting()でグラフの設定をする

setting()はグラフの設定をする関数です。
ユーザーにグラフの設定をしますか？とメッセージを表示し、yesが入力された場合はreset()を、それ以外が入力された場合はusedefult()を実行します。

入力値かデフォルト値か

• reset()
  では関数とグラフの範囲のリストであるaxis_の四要素それぞれに対し入力を受け、ユーザーの求める範囲のリストを返します。
• usedefult()
  では入力を受けず、 axis_には正負/xyの各方向に10ずつの範囲が指定されます。

定義域の決定

setting()の中心付近のコード

global x_min,x_max
x_min = axis_[0]
x_max = axis_[1]+0.1

ここはいわゆる関数の値域を定義しています。
この三行でグローバル変数x_min,x_maxを定義し、次のplot()が動くようにします。
まぁsetting()とplot()を分けなければいいだけなんですが個人的な趣味で分けました

後半：plot()でグラフを表示する

(ここでちょっとだけnumpyを使いました(使わなくてもいいですが))

式の入力と型変換

f = sp.sympify(input("f>>"))

ここで入力された式をsympyで使える形に変換します。

値を求める

y = []
x_ = np.arange(x_min,x_max,step)
for i in x_:
    f_ = f.subs(x,i)
    y.append(f_)

ここがsympy/numpyの使いどころで、さっき定義したx_min,x_maxの範囲を0.1刻みで式に代入していきます。

ここは内包表記というのを使うと

x_ = np.arange(x_min,x_max,step)
y = [f.subs(x,i) for i in x_]

とも書けますが見づらいので今回は使いません。
内包表記は早いですが、正直この程度の計算では誤差なので好きな方を使ってください。

終盤：例外処理で見なかったことにする

ここ最高に頭悪くて好き

try:
    setting()
    plot()
except:
    print("ERROR")

私は数学には疎いので、余計なことをして分からなくなるよりは、ということでまとめて例外に流してありますが、もし数学の問題を解いたりする際のツールとして使うなら、raiseを使って細かに場合分けをすることをお勧めします。

「グラフに関する説明は？」

無いです。
お疲れさまでした。