はじめに

　突然ですが，数式や幾何にミステリアスな魅力を感じることはありませんか？
　かっこいい幾何図形，きれいな曲線を描画したいとき，Web情報¹²やPython Matplotlibが頼りです．スキマ時間の「ちょこっとコーディング」，気楽なちょこっとプログラミングとしても手頃なお題です．図形や空間認識能力を心地よく刺激する効果も期待できます．コードが描くきれいな図形，それを作り出す数式の意味に自然に興味が出たら，とても良い感じです．某情報学部の各フロアにちりばめられた「あの曲線」も含まれます．

いろいろな曲線

　図形の名前でどんな形状か，ビジュアルが頭の中に浮かんできますか？

Lemniscate
Heart Curve
Cardioid
Lissajous Curve
Suddle Curve
Rose Curve
Golden Spiral
Penrose Triangle
Strange Attractor (Lorenz)

動作環境

MacOS sonoma 14.1.2
Python3.11.6
pip 23.3.2
VSCode 1.85.1
JupyterLab 3.4.5
numpy 1.24.4
japanize-matplotlib 1.1.3
matplotlib 3.5.3
matplotlib-inline 0.1.6
plotly 5.13.0

必要に応じて，上記のライブラリをpipでインストールしておきましょう：

    $ python3 -m pip install update pip
    $ python3 -m pip install --upgrade pip
    $ python3 -m pip install numpy matplotlib plotly

図形とコード

　以下，図形の視覚化例と，そのコード例を示します：

Lemniscate

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2 * np.pi, 1000)

def lemniscate(a,theta):
    x = a * np.sqrt(2)*np.cos(theta)/(1 + np.sin(theta)**2)
    y = x * np.sin(theta)
    return x , y

x,y = lemniscate(2,theta)

plt.plot(x, y,linewidth=12,color='gold')
plt.axis('equal')
plt.title('Lemniscate Curve')
plt.show()

Heart Curve

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 2 * np.pi, 100)
x = 16 * np.sin(t) ** 3
y = 13 * np.cos(t) - 5 * np.cos(2 * t) - 2 * np.cos(3 * t) - np.cos(4 * t)

plt.plot(x, y, linewidth=16, color='gold')
plt.axis('equal')
plt.title('Heart Curve')
plt.show()

Cardioid

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0,2*np.pi,300)
r =  3 + 5 * np.cos(theta)

plt.polar(theta, r, 'gold',linewidth=10)
plt.show()

Lissajous

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 2*np.pi, 300)
x = np.sin(t + np.pi/4)
y = np.cos(2*t)

theta = 12
theta = np.deg2rad(theta) 
R = np.array([[np.cos(theta), -np.sin(theta)], 
              [np.sin(theta), np.cos(theta)]])
x_kp, y_kp = np.dot(R, [x, y])

plt.figure(figsize=(4,4))
plt.plot(x_kp, y_kp, 'orange', linewidth=16)
plt.axis('equal')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 2*np.pi, 600)
x = np.cos(2*t + np.pi/2)
y = np.sin(3*t)

theta = 0
theta = np.deg2rad(theta) 
R = np.array([[np.cos(theta), -np.sin(theta)], 
              [np.sin(theta), np.cos(theta)]])
x_kp, y_kp = np.dot(R, [x, y])

plt.figure(figsize=(4,4))
plt.plot(x_kp, y_kp, 'orange', linewidth=8)
plt.axis('equal')
plt.show()

Suddle Point

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

X = np.linspace(-5, 5, 50)
Y = np.linspace(-5, 5, 50)
X, Y = np.meshgrid(X, Y)
Z = X**2 - Y**2

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(X, Y, Z, color='gold', cmap='copper')

plt.show()

Suddle Point (plotly)

Plotlyにより，対話的に視点を変えたり，拡大縮小回転しながら眺めることができます．（JupyterLab, Colaboratory経由でブラウザ表示）

import numpy as np
import plotly.graph_objects as go

def f(x,y):
    return x**2 - y**2

x = np.linspace(-5, 5, 50) 
y = np.linspace(-5, 5, 50) 
X, Y = np.meshgrid(x, y) 

Z = f(X, Y) 

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_layout(title='Saddle Point for x^2 - y^2') 
fig.show()

Rose Curve

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2*np.pi, 1000)

N = 5
r = 7 * np.cos(N * theta)

plt.polar(theta, r, 'gold',linewidth=7)
plt.show()

Golden Spiral

import matplotlib.pyplot as plt
import numpy as np

def polar_to_cartesian(r, theta):
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    return x, y

theta_min = 0
theta_max = 6 * np.pi
theta_step = 0.02
theta = np.arange(theta_min, theta_max, theta_step)

phi = (1 + np.sqrt(5)) / 2  # golden ratio
r = phi ** (theta / np.pi)

x, y = polar_to_cartesian(r, theta)

plt.figure(figsize=(5,5))
plt.plot(x,y,'gold',linewidth=7)
plt.axis('equal')
plt.show()

Penrose Triangle（Penrose三角形³）

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import math
# 重心（中心）に対して，一定角度で回転した座標
def rotate_th_degrees(th,X,Y):
    N = len(X)
    def center(x,y):
        xs,ys = sum(x),sum(y)
        return xs/N, ys/N
    rth = math.radians(th)
    
    Rx,Ry = [],[]
    G = center(X,Y)
    for i in range(N):
        x = X[i] - G[0]
        y = Y[i] - G[1]
        xt = x * math.cos(rth) - y * math.sin(rth)
        yt = x * math.sin(rth) + y * math.cos(rth)
        Rx += [xt + G[0]]
        Ry += [yt + G[1]]
    return Rx,Ry
    
# Penroseの不可能三角形
def penrose():
    k = math.sqrt(3)/2
    S = 2.0
    L = S/8
    H = L * k
    x0 = S/2
    y0 = k*x0
    def p(x,y):
        return list( [ L*x-x0, H*y-y0 ] )
    p1 = p(1.0,0.0)
    p2 = p(7.0,0.0)
    p3 = p(7.5,1.0)
    p4 = p(4.5,7.0)
    p5 = p(3.5,7.0)
    p6 = p(0.5,1.0)
    p7 = p(4.5,5.0)
    p8 = p(3.0,2.0)
    p9 = p(2.0,2.0)
    p10 = p(5.0,2.0)
    p11 = p(5.5,1.0)
    p12 = p(4.0,4.0)
    plist = [p1,p2,p3,p4,p5,p6,p1,p2,p7,p8,p9,p4,p9,p10,p12,p11,p6] 
    return [ p[0] for p in plist ],[ p[1] for p in plist ]

x,y = penrose()
px,py = rotate_th_degrees(180,x,y)

fig = plt.figure()
fig, ax = plt.subplots()
ax.plot(x,y,linewidth=4)
ax.plot(px,py,linewidth=4)
ax.set_aspect('equal')
plt.show()

Penrose三角形を，重心に対して45度ずつ，回転表示

fig = plt.figure()
fig, ax = plt.subplots()

x,y = penrose()
ax.plot(x,y,linewidth=2)
ax.set_aspect('equal')    

for th in range(45,360,45):
    x,y = rotate_th_degrees(th,x,y)
    ax.plot(x,y,linewidth=2)

plt.show()

Strange Attractor (Lorenz)⁴

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import colors

def lorenz(x, y, z, s=8, r=22, b=8/3):
    xd = s*(y - x)
    yd = r*x - y - x*z
    zd = x*y - b*z
    return xd, yd, zd

steps,dt = 18000,0.005
xs = np.empty(steps+1)
ys = np.empty(steps+1)
zs = np.empty(steps+1)
xs[0], ys[0], zs[0] = 1, 1, 1

for i in range(steps):
    xd, yd, zd = lorenz(xs[i], ys[i], zs[i])
    xs[i+1] = xs[i] + xd*dt
    ys[i+1] = ys[i] + yd*dt
    zs[i+1] = zs[i] + zd*dt

ax = plt.figure(figsize=(8,8), dpi=300).add_subplot(projection='3d')
ax.view_init(elev=25, azim=120)
cn = colors.Normalize(ys.min(),ys.max())

for i in range(steps):
    ax.plot(xs[i:i+2], ys[i:i+2], zs[i:i+2], 
            color=plt.cm.viridis(cn(xs[i])), lw=0.3)

ax.set_title('Lorenz "strange" Attractor', size=8)
plt.show()

おわりに

　幾何や数式に興味を抱くきっかけと，スキマ時間にちょこっとコーディングのおすすめ部材として，きれいな曲線，幾何図形を描いてみるのも，ちょっとおしゃれな趣味ではないでしょうか？

参考情報：

いろいろな曲線の確認 (http://izumi-math.jp/S_Yoshida/matome/sc_iroironakyokusen.pdf)
媒介変数表示された有名な曲線７つ(https://manabitimes.jp/math/898)
Penroseの不可能三角形　(https://ja.wikipedia.org/wiki/ペンローズの三角形)
ローレンツ・アトラクタ（https://ja.wikipedia.org/wiki/アトラクター）

いろいろな曲線の確認 (http://izumi-math.jp/S_Yoshida/matome/sc_iroironakyokusen.pdf) ↩
媒介変数表示された有名な曲線７つ(https://manabitimes.jp/math/898) ↩
Penroseの不可能三角形　(https://ja.wikipedia.org/wiki/ペンローズの三角形) ↩
アトラクタ（https://ja.wikipedia.org/wiki/アトラクター） ↩

ちょこっとコーディング：幾何図形や曲線を描く