必要なモジュールの import

In [1]:

from sympy.abc import *
from sympy import *

init_printing()

偏微分

SymPy では，偏微分は１変数の場合の微分と同じ diff() 関数を使って
diff(f, x); や diff(f, y); などと書きます。

In [2]:

# 関数 f(x, y) の定義
f = Function('f')(x, y)

In [3]:

diff(f, x)

Out[3]:

\frac{\partial}{\partial x} f (x, y)

In [4]:

diff(f, y)

Out[4]:

\frac{\partial}{\partial y} f (x, y)

偏微分は交換可能

$\frac{\partial^{2}}{\partial x \partial y} f (x, y) = \frac{\partial^{2}}{\partial y \partial x} f (x, y)$

です。

$\frac{\partial^{2}}{\partial y \partial x} f (x, y) - \frac{\partial^{2}}{\partial x \partial y} f (x, y) = 0$

を示してみます。

In [5]:

Eq(Derivative(Derivative(f, x), y) - Derivative(Derivative(f, y), x),
   diff(diff(f, x), y) - diff(diff(f, y), x))

Out[5]:

\frac{\partial^{2}}{\partial y \partial x} f (x, y) - \frac{\partial^{2}}{\partial x \partial y} f (x, y) = 0

スカラー場の勾配（grad）

スカラー場 $φ (r)$ の勾配を定義する。

$\nabla φ \equiv grad φ = (\frac{\partial φ}{\partial x}, \frac{\partial φ}{\partial y}, \frac{\partial φ}{\partial z})$

In [6]:

def grad(phi):
    return Matrix([diff(phi,x), diff(phi,y), diff(phi,z)])

勾配（grad）の計算例

$r = \sqrt{x^{2} + y^{2} + z^{2}}$ のとき，以下の計算をしてみます。

$\nabla r = \frac{\partial r}{\partial x} i + \frac{\partial r}{\partial y} j + \frac{\partial r}{\partial z} k$

In [7]:

rvec = Matrix([x, y, z])
r = sqrt(rvec.dot(rvec))

In [8]:

grad(r)

Out[8]:

[\begin{matrix} \frac{x}{\sqrt{x^{2} + y^{2} + z^{2}}} \\ \frac{y}{\sqrt{x^{2} + y^{2} + z^{2}}} \\ \frac{z}{\sqrt{x^{2} + y^{2} + z^{2}}} \end{matrix}]

上の結果が

$\nabla r = \frac{x}{r} i + \frac{y}{r} j + \frac{z}{r} k = \frac{r}{r}$

となっていることを確認します。

左辺から右辺を引いてゼロベクトルになっていればいいですから…

In [9]:

grad(r) - rvec/r

Out[9]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

また，電磁気学では以下のような計算をする必要も出てくるので，やっておきます。

$\nabla (\frac{1}{r}) = - \frac{r}{r^{3}}$

In [10]:

grad(1/r)

Out[10]:

[\begin{matrix} - \frac{x}{{(x^{2} + y^{2} + z^{2})}^{\frac{3}{2}}} \\ - \frac{y}{{(x^{2} + y^{2} + z^{2})}^{\frac{3}{2}}} \\ - \frac{z}{{(x^{2} + y^{2} + z^{2})}^{\frac{3}{2}}} \end{matrix}]

In [11]:

_ + rvec/r**3

Out[11]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

ベクトル場の発散 (div)

ベクトル場の発散を定義する。

$\nabla \cdot v \equiv div v = \frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} + \frac{\partial v_{z}}{\partial z}$

In [12]:

def div(v):
    return diff(v[0], x) + diff(v[1], y) + diff(v[2], z)

発散 (div) の計算例

$\nabla \cdot r = \frac{\partial x}{\partial x} + \frac{\partial y}{\partial y} + \frac{\partial z}{\partial z} = 1 + 1 + 1 = 3$

（答えは空間の次元の数）

In [13]:

div(rvec)

Out[13]:

3

ベクトル場の回転 (rot)

ベクトルの回転を以下のようにして定義します。

$\begin{array}{rcl} \nabla \times v \equiv rot v & = & (\frac{\partial}{\partial y} v_{z} - \frac{\partial}{\partial z} v_{y}, \frac{\partial}{\partial z} v_{x} - \frac{\partial}{\partial x} v_{z}, \frac{\partial}{\partial x} v_{y} - \frac{\partial}{\partial y} v_{x}) \end{array}$

In [14]:

def rot(v):
    vx = v[0]
    vy = v[1]
    vz = v[2]
    return Matrix(
            [diff(vz, y) - diff(vy, z), 
             diff(vx, z) - diff(vz, x), 
             diff(vy, x) - diff(vx, y)]
            )

回転 (rot) の計算例

$\nabla \times (\frac{r}{r})$

In [15]:

rot(rvec/r)

Out[15]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

スカラー場やベクトル場の2階微分と恒等式

勾配 (grad) の発散 (div)

スカラー場 $φ$ の勾配 $\nabla φ$ の発散 $\nabla \cdot (\nabla φ)$

In [16]:

phi = Function('phi')(x, y, z)

In [17]:

div(grad(phi))

Out[17]:

\frac{\partial^{2}}{\partial x^{2}} ϕ (x, y, z) + \frac{\partial^{2}}{\partial y^{2}} ϕ (x, y, z) + \frac{\partial^{2}}{\partial z^{2}} ϕ (x, y, z)

ラプラス演算子の計算例

$\nabla^{2} (\frac{1}{r}) = 0 (r \neq 0)$

を確認します。

In [18]:

def Laplacian(phi):
    # スカラー関数 phi のみに適用
    return div(grad(phi))

In [19]:

Laplacian(1/r)

Out[19]:

\frac{3 x^{2}}{{(x^{2} + y^{2} + z^{2})}^{\frac{5}{2}}} + \frac{3 y^{2}}{{(x^{2} + y^{2} + z^{2})}^{\frac{5}{2}}} + \frac{3 z^{2}}{{(x^{2} + y^{2} + z^{2})}^{\frac{5}{2}}} - \frac{3}{{(x^{2} + y^{2} + z^{2})}^{\frac{3}{2}}}

In [20]:

simplify(_)

Out[20]:

0

勾配 (grad) の回転 (rot)

スカラー場 $φ$ の勾配 $\nabla φ$ の回転 $\nabla \times (\nabla φ)$ 。

恒等的に
$\nabla \times (\nabla φ) = 0$ である。

In [21]:

grad(phi)

Out[21]:

[\begin{matrix} \frac{\partial}{\partial x} ϕ (x, y, z) \\ \frac{\partial}{\partial y} ϕ (x, y, z) \\ \frac{\partial}{\partial z} ϕ (x, y, z) \end{matrix}]

In [22]:

rot(grad(phi))

Out[22]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

回転 (rot) の発散 (div)

ベクトル場 $a$ の回転 $\nabla \times a$ の発散 $\nabla \cdot (\nabla \times a)$

恒等的に
$\begin{array}{rcl} \nabla \cdot (\nabla \times a) & = & 0 \end{array}$
である。

In [23]:

a1 = Function('a1')(x,y,z)
a2 = Function('a2')(x,y,z)
a3 = Function('a3')(x,y,z)

In [24]:

a = [a1, a2, a3]

In [25]:

rot(a)

Out[25]:

[\begin{matrix} - \frac{\partial}{\partial z} a_{2} (x, y, z) + \frac{\partial}{\partial y} a_{3} (x, y, z) \\ \frac{\partial}{\partial z} a_{1} (x, y, z) - \frac{\partial}{\partial x} a_{3} (x, y, z) \\ - \frac{\partial}{\partial y} a_{1} (x, y, z) + \frac{\partial}{\partial x} a_{2} (x, y, z) \end{matrix}]

In [26]:

div(rot(a))

Out[26]:

0

回転 (rot) の回転 (rot)

ベクトル場 $a$ の回転 $\nabla \times a$ の回転 $\nabla \times (\nabla \times a)$

$\nabla \times (\nabla \times a) = \nabla (\nabla \cdot a) - \nabla^{2} a$

が成り立つことを確認する。

In [27]:

# ベクトルに作用する Laplacian を定義
def Lapv(v):
    return Matrix([Laplacian(v[0]),
                   Laplacian(v[1]),
                   Laplacian(v[2])])

In [28]:

rot(rot(a)) - (grad(div(a)) - Lapv(a))

Out[28]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

追記：外積の発散

一般に，2つのベクトル $E, B$ の外積 $E \times B$ の発散は，それぞれのベクトルの回転を使って以下のように書くことができる。

$\nabla \cdot (E \times B) = (\nabla \times E) \cdot B - E \cdot (\nabla \times B)$

この式は，電磁場のエネルギー密度・エネルギー流束のところで出てくる。

In [29]:

E1 = Function('E1')(x,y,z)
E2 = Function('E2')(x,y,z)
E3 = Function('E3')(x,y,z)
B1 = Function('B1')(x,y,z)
B2 = Function('B2')(x,y,z)
B3 = Function('B3')(x,y,z)

E = Matrix([E1, E2, E3])
B = Matrix([B1, B2, B3])

In [30]:

# 左辺から右辺を引いてゼロとなることを示す

div(E.cross(B)) - (rot(E).dot(B) - E.dot(rot(B)));

In [31]:

simplify(_)

Out[31]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

追記：スカラーとベクトルの積の発散

$\nabla \cdot (ϕ D) = (\nabla ϕ) \cdot D + (\nabla \cdot D) ϕ$

この式も，電磁場のエネルギー密度・エネルギー流束のところで出てくる。

In [32]:

D1 = Function('D1')(x,y,z)
D2 = Function('D2')(x,y,z)
D3 = Function('D3')(x,y,z)

D = Matrix([D1, D2, D3])

In [33]:

# 左辺から右辺を引いてゼロとなることを示す

div(phi * D) - (grad(phi).dot(D) + div(D) * phi);

In [34]:

simplify(_)

Out[34]:

[\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}]

参考：SymPy でスカラー場・ベクトル場の微分