Study math with Python: Draw a sympy (scipy) graph with matplotlib
Visualization is important
It is often difficult to intuitively understand what the result will look like even if it is obtained by a function. Therefore, drawing a graph is the first step for analysis and thinking. For the Japanese memo of matplotlib, see here. The sympy plot is here. Official I struggled while looking at it.
import
import sympy as sym
sym.init_printing()
Pi = sym.S.Pi #Pi
E = sym.S.Exp1 #The bottom of the natural logarithm
I = sym.S.ImaginaryUnit #Imaginary unit
#Definition of variables to use(All lowercase letters are symbols)
(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) = sym.symbols('a b c d e f g h i j k l m n o p q r s t u v w x y z')
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
Do your best in color
The color list was pre-generated by another program so that it would be a ring with constant brightness in the L * a * b * color space. Designed with a safety margin for contrast rather than vividness of color, color blindness is not considered. Each list consists of 37 HEX strings. The number 37 was selected because it is the prime number closest to 36, which is the number that advances 360 degrees in one lap by 10 degrees. (If it is 37, it advances 9.73 degrees per one.) Since it is a prime number, no matter how you choose the initial color and the number of steps (it must be a primitive root), it will always be a ** cyclic group ** that passes through 37 colors one by one. This is an intuitive match that the phase advances by about 10 steps * 10 degrees. As much as possible, the same color does not appear and the cycle is long.
def generator_p(start): #Generator to generate colors in sequence if arguments are specified, otherwise
prime = 37 #Number of lists(Since it is a prime number, it will always be the maximum period no matter what you specify for primitive.)
primitive = 8 #Source(Hue per step is primitive* 37/Go 360 degrees)
g0 = start #Starting source(What is actually returned is g0* primitive)
while True:
val = yield g0
if val: #If an argument is passed to the generator, it will be reset as the starting source
g0 = val
g0 += primitive
if g0 >= prime: #Take the surplus
g0 -= prime
gen_hexnum = generator_p(0) #Put the initial color in the argument
gen_hexnum.__next__() #Generator initialization
def hexNum(num, type): #List of colors
#Basic color(A calm color that can be used universally on a white background)
color_basic = ['#9c1954', '#9e1a46', '#9d1e38', '#9a252a', '#962d1c', '#8f350d', '#873c00',
'#7d4300', '#734900', '#674f00', '#5b5300', '#4d5700', '#3e5b00', '#2b5d00',
'#0f6009', '#00611c', '#00632b', '#00643a', '#006449', '#006558', '#006567',
'#006575', '#006482', '#00648e', '#006298', '#0060a0', '#005ea6', '#005baa',
'#0056aa', '#0051a9', '#2f4ca4', '#50459d', '#673d94', '#78358a', '#862d7d', '#902470', '#981e62']
#Saturated colors(It stands out, but you shouldn't use it a lot.)
color_vivid = ['#ffadc7', '#ffadbc', '#ffaeb2', '#ffb0a8', '#ffb29e', '#ffb596', '#f9b88f',
'#f1bc8a', '#e9bf86', '#dfc385', '#d5c685', '#caca87', '#becd8b', '#b2cf90',
'#a6d298', '#99d4a0', '#8dd5aa', '#80d7b4', '#74d7bf', '#6ad8ca', '#61d8d5',
'#5bd7e0', '#59d6e9', '#5dd5f2', '#65d3f9', '#71d1ff', '#7fceff', '#8ecbff',
'#9ec8ff', '#afc4ff', '#bec1ff', '#cdbdfd', '#dab9f7', '#e6b6ef', '#f0b3e6', '#f9b0dc', '#ffaed2']
#Light color(Not used on a white background. Expected to be used on a black background.)
color_thin = ['#ffc3d5', '#ffc3cd', '#ffc3c5', '#ffc4be', '#ffc6b7', '#ffc8b1', '#fbcaac',
'#f5cca9', '#efcfa6', '#e8d2a5', '#e0d4a5', '#d8d7a6', '#cfd9a9', '#c6dbad',
'#bdddb3', '#b5deb9', '#ace0c0', '#a5e0c7', '#9ee1cf', '#98e1d7', '#94e1df',
'#92e1e6', '#92e0ed', '#94dff4', '#99def9', '#9fdcfd', '#a7daff', '#b1d8ff',
'#bbd5ff', '#c5d3ff', '#cfd0ff', '#dacdfc', '#e3cbf7', '#ecc8f2', '#f3c6eb', '#f9c5e4', '#fec4dc']
if (type == 'thin'):
hex_list = color_thin
elif (type == 'vivid'):
hex_list = color_vivid
else:
hex_list = color_basic
if num is not None:
gen_hexnum.send(num)
return hex_list[gen_hexnum.__next__()]
Creating a 2D graph plot function: 1-variable explicit function
def plot(func, xrange, ylim=None, xlog=None, ylog=None, cn=None, ct='basic', close=None, figs=None, simbol=x):
if figs:
plt.figure(figsize=figs)
X = np.arange(xrange[0], xrange[1], xrange[2])
if hasattr(func, "subs"): # sympy
Y = [func.subs([(simbol, K)]) for K in X]
else: # scipy
Y = [func(K) for K in X]
if ylim:
plt.ylim(ylim)
if xlog:
plt.xscale('log')
else:
plt.xscale('linear')
if ylog:
plt.yscale('log')
else:
plt.yscale('linear')
plt.plot(X, Y, color=hexNum(cn, ct))
if close:
plt.show()
plt.close()
The plain matplotlib.plot is confusing, so I created a function that wraps only frequently used arguments. It's easier to use if you give the function, width, and color, and then ask for it.
F = sym.sin(x)/x
G = sym.cos(x)/x
plot(F, (-10, 10, 0.1), cn=0, figs=(10, 6))
plot(G, (-10, 10, 0.1))
plot(F+0.2, (-10, 10, 0.1))
plot(F+0.4, (-10, 10, 0.1))
plot(F+0.6, (-10, 10, 0.1))
plot(F+0.8, (-10, 10, 0.1))
plot(F+1.0, (-10, 10, 0.1))
def A(x):
return np.sin(x)/x
plot(A, (-10, 10, 0.1))
plot(lambda X: (X/9.0)**2, (-10, 10, 0.1), (-0.5, 1.5), close=True)
To use sympy comfortably, I want to use all lowercase letters as symbols, so
Function names and loop variables are intentionally capitalized.
Color-coded effect by color, automatic judgment of sympy and numpy functions by hasattr,
With the close argument, all the options have to be written together at the end.
I am satisfied because I can plot it quite intuitively.
However, looking at the graph, the poles of $ \ frac {\ cos (x)} {x} $ are strange.
The atmosphere is that the line in the middle of $-\ infinty \ longrightarrow \ infinty $ is drawn.
Is there a solution?
Creation of 2D graph plot function: Parameter display
def plotp(func, t_range=(-4, 4, 100), axis=[-4, 4, -4, 4],
xlog=None, ylog=None, cn=None, ct='basic', close=None, figs=None, simbol=t):
if figs:
plt.figure(figsize=figs)
X, Y = func
t_ndarray = np.linspace(t_range[0], t_range[1], t_range[2]) # ndarray
if hasattr(X, "subs"): # sympy
Xp = np.array([float(sym.N(X.subs([(simbol, T)]))) for T in t_ndarray])
Yp = np.array([float(sym.N(Y.subs([(simbol, T)]))) for T in t_ndarray])
else: # scipy
Xp = [X(K) for K in t_ndarray]
Yp = [Y(K) for K in t_ndarray]
plt.axis(axis)
if xlog:
plt.xscale('log')
else:
plt.xscale('linear')
if ylog:
plt.yscale('log')
else:
plt.yscale('linear')
plt.plot(Xp, Yp, color=hexNum(cn, ct))
if close:
plt.show()
plt.close()
The difference from plot () is that Func is passed as a list of two functions, with the addition of an axis that specifies the x and y regions. By default, the symbol defaults to t. Also, as the third range of t, it is specified by the number of divisions instead of the step size.
X = t**2 * sym.sin(t)
Y = t**2 * sym.cos(t)
plotp([X, Y], (-4 * np.pi, 4 * np.pi, 400), [-100, 100, -100, 100], figs=(6, 6), close=True)
plotp([lambda T: T * np.sin(T), lambda T: T * np.cos(T)], (-4 * np.pi, 4 * np.pi, 400), [-12, 12, -10, 10],
cn=14, figs=(6, 6), close=True)
It worked for both sympy and scipy.
Creating a 2D graph plot function: 2-variable implicit function
It is a function expressed in the format of $ F (x, y) = 0 $. It seems difficult to handle implicit functions directly. So use sympy's plot_implicit as it is.
F = x**3 + y**2 -1
sym.plot_implicit(F, (x, -2, 2), (y, -2, 2), line_color=hexNum(2, 'basic'))
This is an elliptic curve equation.
The color seems to be usable, but it is light.
F = x**3 - 2 * x + y**2 > 0
sym.plot_implicit(F, (x, -2, 2), (y, -2, 2), line_color=hexNum(2, 'thin'))
The advantage is that you can use inequalities.
Creating a 3D graph plot function: 2-variable explicit function
A function of the form $ z = f (x, y) $.
from matplotlib import cm
def plot3D(func, xrange=(-2, 2, 0.1), yrange=(-2, 2, 0.1), ptype='plot', simbol=[x, y]):
fig = plt.figure()
ax = Axes3D(fig)
if ptype == 'plot':
ax.plot(func[0], func[1], func[2])
plt.show()
plt.close()
return
X = np.arange(xrange[0], xrange[1], xrange[2])
Y = np.arange(yrange[0], yrange[1], yrange[2])
Xmesh, Ymesh = np.meshgrid(X, Y)
if hasattr(func, "subs"): # sympy
Z = np.array([np.array([float(sym.N(func.subs([(simbol[0], K), (simbol[1], L)]))) for K in X]) for L in Y])
else: # scipy
Z = func(Xmesh, Ymesh)
if ptype == 'surface':
ax.plot_surface(Xmesh, Ymesh, Z, cmap=cm.bwr)
elif ptype == 'wireframe':
ax.plot_wireframe(Xmesh,Ymesh, Z)
elif ptype == 'contour':
ax.contour3D(X, Y, Z)
plt.show()
plt.close()
plot3D(lambda T, U: T * U * np.sin(T), (-4, 4, 0.25), (-4, 4, 0.25), 'surface')
plot3D((x ** 2 - 3 - y ** 2), (-4, 4, 0.25), (-4, 4, 0.25), 'surface')
plot3D((x ** 2 - 3 - y ** 2), (-4, 4, 0.25), (-4, 4, 0.25), 'wireframe')
plot3D((x ** 2 - 3 - y ** 2), (-4, 4, 0.25), (-4, 4, 0.25), 'contour')
The default color is used for 3D because it is better to see the ups and downs with a difference in brightness. I plotted the function of the horse saddle that often appears in the figure of the universe with negative curvature.
Creation of 3D graph plot function: Parameter display
def plotp3(func, simbol=[t], rangeL=[(-4, 4, 100)]):
X, Y, Z = func
t1 = simbol[0]
if len(simbol) >= 2:
u1 = simbol[1]
t_rangeL = rangeL[0]
t_ndarray = np.linspace(t_rangeL[0], t_rangeL[1], t_rangeL[2]) # ndarray
if len(rangeL) >= 2:
u_rangeL = rangeL[1]
u_ndarray = np.linspace(u_rangeL[0], u_rangeL[1], u_rangeL[2]) # ndarray
if len(rangeL) == 1:
Xp = np.array([float(sym.N(X.subs([(t1, T)]))) for T in t_ndarray])
Yp = np.array([float(sym.N(Y.subs([(t1, T)]))) for T in t_ndarray])
Zp = np.array([float(sym.N(Z.subs([(t1, T)]))) for T in t_ndarray])
else:
Xp = np.array([float(sym.N(X.subs([(t1, T), (u1, U)]))) for T in t_ndarray for U in u_ndarray])
Yp = np.array([float(sym.N(Y.subs([(t1, T), (u1, U)]))) for T in t_ndarray for U in u_ndarray])
Zp = np.array([float(sym.N(Z.subs([(t1, T), (u1, U)]))) for T in t_ndarray for U in u_ndarray])
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(Xp, Yp, Zp)
plt.show()
plt.close()
X = 2*sym.cos(t)
Y = 5*sym.sin(t)
Z = t
plotp3([X, Y, Z], [t], [(-5 * np.pi, 5 * np.pi, 200)])
X = (5+2*sym.cos(t))*sym.cos(u)
Y = (5+2*sym.sin(u))*sym.sin(u)
Z = 2*sym.sin(t)-u/sym.oo
plotp3([X, Y, Z], [t, u], [(-np.pi, np.pi, 40), (-np.pi, np.pi, 40)])
I'm getting tired from this area. 3D is sympy default, so it's okay.
X = 2*sym.cos(t)
Y = 5*sym.sin(t)
Z = t
sym.plotting.plot3d_parametric_line(X, Y, Z, (t, -5*np.pi, 5*np.pi))
X = (5+2*sym.cos(t))*sym.cos(u)
Y = (5+2*sym.sin(u))*sym.sin(u)
Z = 2*sym.sin(t)-u/sym.oo
sym.plotting.plot3d_parametric_surface(X, Y, Z, (t, -np.pi, np.pi), (u, -np.pi, np.pi))
Recommended Posts