20240120_nabla.md

ナブラ記号とベクトル解析

~draw-rotation.py

def draw_rotation(xx, yy, vxx, vyy, d):
    rot_zz = vyy[1:-1,2:] - vyy[1:-1,:-2] - vxx[2:,1:-1] + vxx[:-2,1:-1]
    rot_zz /= 2 * d
    vmax = np.abs(rot_zz).max()
    if vmax < 1e-4:
        vmax = 1.0
    plt.title(r'${:.3f}\ \leqq\ \mathrm{{rot}}\ \vec{{v}}\ \leqq\ {:.3f}$'.format(rot_zz.min(), rot_zz.max()))
    plt.xlabel('$x$', fontsize=12)
    plt.ylabel('$y$', fontsize=12)
    plt.imshow(rot_zz, cmap='bwr', origin='lower', extent=(xx.min(), xx.max(), yy.min(), yy.max()), vmin=-vmax, vmax=vmax)
    plt.colorbar()

def draw_vector_field_and_rotation(xx, yy, vxx, vyy, d):
    plt.figure(figsize=(10, 4))
    plt.subplots_adjust(wspace=0.4, hspace=0.4)  # グラフ間の余白の大きさ
    plt.subplot(1, 2, 1)
    draw_vector_field(xx, yy, vxx, vyy)
    plt.subplot(1, 2, 2)
    draw_rotation(xx, yy, vxx, vyy, d)
    plt.show()


d = 0.01
xx, yy = np.meshgrid(np.arange(-1, 1+d, d), np.arange(-1, 1+d, d))

# 回転なしの例
draw_vector_field_and_rotation(xx, yy, vxx=xx*2, vyy=yy, d=d)
draw_vector_field_and_rotation(xx, yy, vxx=np.sin(xx), vyy=np.cos(yy), d=d)
draw_vector_field_and_rotation(xx, yy, vxx=yy, vyy=xx, d=d)

# 回転ありの例
draw_vector_field_and_rotation(xx, yy, vxx=-yy, vyy=xx, d=d)
draw_vector_field_and_rotation(xx, yy, vxx=np.sin(np.pi*yy), vyy=np.cos(np.pi*xx), d=d)
draw_vector_field_and_rotation(xx, yy, vxx=np.zeros(xx.shape), vyy=xx**2, d=d)

~draw-vector-field.py

import numpy as np
from matplotlib import pyplot as plt

def draw_vector_field(xx, yy, vxx, vyy, interval=20, scale=0.2):
    # ベクトル場の描画
    x = xx[::interval, ::interval].flatten()
    y = yy[::interval, ::interval].flatten()
    vx = vxx[::interval, ::interval].flatten() * scale
    vy = vyy[::interval, ::interval].flatten() * scale
    x_start = x - vx/2
    y_start = y - vy/2
    plt.title('Vector Field')
    plt.xlabel('$x$', fontsize='12')
    plt.ylabel('$y$', fontsize='12')
    plt.scatter(x, y, s=10, color='black')
    for x_start_, y_start_, vx_, vy_ in zip(x_start, y_start, vx, vy):
        plt.arrow(x_start_,y_start_,vx_,vy_, head_width=0.02, color='black', ec='black')

~draw_divergence.py

def draw_divergence(xx, yy, vxx, vyy, d):
    div = vxx[1:-1,2:] - vxx[1:-1,:-2] + vyy[2:,1:-1] - vyy[:-2,1:-1]
    div /= 2 * d
    vmax = np.abs(div).max()
    if vmax < 1e-4:
        vmax = 1.0
    plt.title(r'${:.3f}\ \leqq\ \mathrm{{div}}\ \vec{{v}}\ \leqq\ {:.3f}$'.format(div.min(), div.max()))
    plt.xlabel('$x$', fontsize=12)
    plt.ylabel('$y$', fontsize=12)
    plt.imshow(div, cmap='bwr', origin='lower', extent=(xx.min(), xx.max(), yy.min(), yy.max()), vmin=-vmax, vmax=vmax)
    plt.colorbar()

def draw_vector_field_and_divergence(xx, yy, vxx, vyy, d):
    plt.figure(figsize=(10, 4))
    plt.subplots_adjust(wspace=0.4, hspace=0.4)  # グラフ間の余白の大きさ
    plt.subplot(1, 2, 1)
    draw_vector_field(xx, yy, vxx, vyy)
    plt.subplot(1, 2, 2)
    draw_divergence(xx, yy, vxx, vyy, d)
    plt.show()


d = 0.01
xx, yy = np.meshgrid(np.arange(-1, 1+d, d), np.arange(-1, 1+d, d))

# 発散なしの例
draw_vector_field_and_divergence(xx, yy, xx, -yy, d=d)
draw_vector_field_and_divergence(xx, yy, -yy, xx, d=d)
draw_vector_field_and_divergence(xx, yy, np.sin(2*np.pi*yy), np.sin(2*np.pi*xx), d=d)

# 発散ありの例
draw_vector_field_and_divergence(xx, yy, xx, yy, d=d)
draw_vector_field_and_divergence(xx, yy, np.sin(2*np.pi*xx), np.sin(2*np.pi*yy), d=d)

~picture-of-divergence.py

import numpy as np
from matplotlib import pyplot as plt

plt.xlabel(r'$x$', fontsize=12)
plt.ylabel(r'$y$', fontsize=12)
plt.scatter(0, 0, label=r'$(x,y)$', color='black')
c = 'red'
plt.scatter(0.5, 0, color=c)
plt.arrow(0.5-0.3,0,0.6,0, label=r'$v_x(x+\Delta,y)$', width=0.01, color=c, ec=c)
c = 'blue'
plt.scatter(-0.5, 0, color=c)
plt.arrow(-0.5-0.3,0,0.6,0, label=r'$v_x(x-\Delta,y)$', width=0.01, color=c, ec=c)
c = 'orange'
plt.scatter(0, 0.5, color=c)
plt.arrow(0,0.5-0.3,0,0.6, label=r'$v_y(x,y+\Delta)$', width=0.01, color=c, ec=c)
c = 'green'
plt.scatter(0, -0.5, color=c)
plt.arrow(0,-0.5-0.3,0,0.6, label=r'$v_y(x,y-\Delta)$', width=0.01, color=c, ec=c)
ax = plt.gca()
ax.axes.xaxis.set_ticklabels([])
ax.axes.yaxis.set_ticklabels([])
plt.legend()
plt.grid()
plt.show()

~picture-of-rotation.py

import numpy as np
from matplotlib import pyplot as plt

plt.xlabel(r'$x$', fontsize=12)
plt.ylabel(r'$y$', fontsize=12)
plt.plot([1,1,-1,-1,1], [-1,1,1,-1,-1], color='black')
plt.scatter(0, 0, s=100, label=r'$(x,y,z)$', color='black')
c = 'red'
plt.scatter(1, 0, s=100, color=c)
plt.arrow(1,-0.2,0,0.4, label=r'$v_y(x+\Delta/2,y,z)$', width=0.02, color=c, ec=c)
c = 'blue'
plt.scatter(-1, 0, s=100, color=c)
plt.arrow(-1,-0.2,0,0.4, label=r'$v_y(x-\Delta/2,y,z)$', width=0.02, color=c, ec=c)
c = 'orange'
plt.scatter(0, 1, s=100, color=c)
plt.arrow(-0.2,1,0.4,0, label=r'$v_x(x,y+\Delta/2,z)$', width=0.02, color=c, ec=c)
c = 'green'
plt.scatter(0, -1, s=100, color=c)
plt.arrow(-0.2,-1,0.4,0, label=r'$v_x(x,y-\Delta/2,z)$', width=0.02, color=c, ec=c)
ax = plt.gca()
ax.axes.xaxis.set_ticklabels([])
ax.axes.yaxis.set_ticklabels([])
plt.legend(loc='upper right')
plt.grid()
plt.show()