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

はじめに

　突然ですが，数式や幾何にミステリアスな魅力を感じることはありませんか？
　かっこいい幾何図形，きれいな曲線を描画したいとき，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/アトラクター）

1. いろいろな曲線の確認 (http://izumi-math.jp/S_Yoshida/matome/sc_iroironakyokusen.pdf) ↩

2. 媒介変数表示された有名な曲線７つ(https://manabitimes.jp/math/898) ↩

3. Penroseの不可能三角形　(https://ja.wikipedia.org/wiki/ペンローズの三角形) ↩

4. アトラクタ（https://ja.wikipedia.org/wiki/アトラクター） ↩