python - 可視化函數

2020-04-22
#code
Project

昨天一個心血來潮，想自己做看看可視化函數。

我一開始想用C++寫，但是想到python在這方面似乎比較簡單，所以我就決定使用python來做(懶)

我是用Anaconda做的，畢竟裡面都已經裝好了，不用白不用。

首先我第一個做的是二維平面的圖:

code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
import matplotlib.pyplot as plt

# function
f = lambda x: x * x

# 產生一段從-10到10的數字，間隔 0.125
x = np.arange(-10, 10, .125)
y = f(x)

plt.figure()
plt.xlabel("x")
plt.ylabel("y")
plt.plot(x, y)
plt.show()

非常之簡單。

接著是兩個圖交疊的:

code:

1
2
3
4
5
6
7
8
9
10
11
# 將上面的code改一下就好了
f = lambda x: 10 * np.sin(x)
g = lambda x: x * x

x = np.arange(-10, 10, .125)
y1 = f(x)
y2 = g(x)

plt.plot(x, y1)
plt.plot(x, y2)
plt.show()

似乎有點過於簡單了?

接著是3D的圖形，根據我上網找的資料，需要用到mpl_toolkits.mplot3d裡面的Axes3D。
因為現在是x, y會產生一個面，然後再根據在面上的某個位置，決定畫在Z軸的哪裡。
所以我們需要np.meshgrid(x, y)這個函式來產生這個面。

三維:

code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()

# ax = Axes3D(fig)
# 2022.5.18 更新
# 上述寫法已經被拋棄，應該改成
ax = plt.axes(projection = "3d")

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

x = np.arange(-5, 5, .125)
y = np.arange(-5, 5, .125)

# 產生輸入面
X, Y = np.meshgrid(x, y)

Z = X + Y

ax.plot_surface(X, Y, Z)
plt.show()

覺得好像只是單一顏色好像太單調了。
順便也把圖形弄得複雜一點。

$f(x, y) = x^2+y^2$:

$f(x, y) = \sqrt{x^2+y^2}$:

$f(x, y) = cos(\sqrt{x^2+y^2})$:

$f(x, y) = sin(\sqrt{x^2+y^2})$:

$f(x, y) = sin(x \times y)$ :

$f(x, y) = sin(exp(x^2 + y^2))$: (xy平面上視角)

code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
fig = plt.figure()
# ax = Axes3D(fig)
ax = plt.axes(projection = "3d")

ax.view_init(45, 45)

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

ax.set_zlim3d(-3, 3)

x = np.arange(-5, 5, .125)
y = np.arange(-5, 5, .125)

X, Y = np.meshgrid(x, y)

Z = f(X, Y)

# rstride, cstride 相當於row, column的解析度，數值越高解析度越低，可以自己調調看
# cmap則是顏色映射表，有許多預設種類
ax.plot_surface(X, Y, Z, rstride = 1, cstride = 1, cmap = 'rainbow')
plt.show()

我們還有一個複數域沒用到，但是當我試著想畫出複數域時，問題開始變得棘手了…
當我們輸入進一個複數$a+bi$時需要兩個變數$(a, b)$去記錄，而函數輸出是$c+di$也需要兩個變數$(c, d)$去記錄。
很明顯的我們需要四個軸來表示，但是我們就算畫出來了，我們也無法理解(吧，還是只有我無法理解?

我上網找到幾個方法:

  1. 只畫實部
  2. 只畫虛部
  3. 畫在複數平面上，對於原點的距離
  4. 畫複數平面上的方位角
  5. 用顏色代表第四維度
  6. Conformal map，共形映射

前四個方法都需要犧牲一點東西。
只畫實部、只畫虛部就不多講了。值得注意的是，通常會把方位角代進sin()，以此讓它變成連續的。
PS. 可以使用mpmath的arg()算方位角

$f(x) = x^2 + 1$(未用sin):

$f(x) = x^2 + 1$(使用sin):

巴尼斯G函數的複數平面上的方位角(未用sin):

巴尼斯G函數的複數平面上的方位角(使用sin):

第五個方法需要用到colormap，彩色映射的方法。

code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import mpmath

# f(x) = x ** 2 + 1
f = lambda x: np.power(x, 2) + 1

fig = plt.figure()
ax = plt.axes(projection = "3d")

# 設定視角
ax.view_init(45, 20)

# f(a + bi) = c + di
ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_zlabel('c')

a = np.arange(-5, 5, .125)
b = np.arange(-5, 5, .125)

A, B = np.meshgrid(a, b)
an, bn = A.shape
C, D = A * 0, A * 0

for i in range(an):
    for j in range(bn):
        try:
            z = mpmath.mp.mpc(A[i, j], B[i, j])
            w = f(z)
            
            C[i, j] = w.real
            D[i, j] = w.imag
            
        except (ValueError, TypeError, ZeroDivisionError):
            # do something...
            pass

# 用顏色代表虛部
color_dimension = D
minn, maxx = color_dimension.min(), color_dimension.max()
norm = matplotlib.colors.Normalize()
m = plt.cm.ScalarMappable(norm = norm, cmap = 'jet')
fcolors = m.to_rgba(color_dimension)

surf = ax.plot_surface(A, B, C, rstride = 1, cstride = 1, facecolors = fcolors, vmin = minn, vmax = maxx)

# 右邊那條
plt.colorbar(m)
plt.show()

$f(x) = x^2 + 1$

$f(x) = x^\pi$

$f(x) = sin(x) \times cos(x)$

至於第六個方法，待研究…