向量微積分
梯度
\nabla_{\boldsymbol{x}} f(\boldsymbol{x})\overset{\underset{\mathrm{def}}{}}{=} \left[ \frac{\partial f(\boldsymbol{x})}{\partial x_1}, \frac{\partial f(\boldsymbol{x})}{\partial x_2},\cdots,\frac{\partial f(\boldsymbol{x})}{\partial x_n} \right]^T=\frac{\partial f(\boldsymbol{x})}{\partial \boldsymbol{x}}
grad 是最斜的方向
grad(f,p) 就是函數 f 在 p 點上最斜的方向。
散度
{\displaystyle \operatorname {div} \mathbf {A} =\nabla \cdot \mathbf {A} ={\frac {\partial A_{x}}{\partial x}}+{\frac {\partial A_{y}}{\partial y}}+{\frac {\partial A_{z}}{\partial z}}}
以散度 divA 就是單位體積向量場 A 之通量 (flux)
散度 ⇐⇒ 通量(曲面積分) ⇐⇒ Gauss 散度定理
旋度
{\displaystyle \mathbf {curl\,} \ \mathbf {A} ={\boldsymbol {\nabla }}\times \mathbf {A} =\left({\frac {\partial A_{z}}{\partial y}}-{\frac {\partial A_{y}}{\partial z}}\right)\mathbf {i} +\left({\frac {\partial A_{x}}{\partial z}}-{\frac {\partial A_{z}}{\partial x}}\right)\mathbf {j} +\left({\frac {\partial A_{y}}{\partial x}}-{\frac {\partial A_{x}}{\partial y}}\right)\mathbf {k} }
旋度的行列式表示法
{\displaystyle \mathbf {curl\,} \mathbf {A} ={\begin{vmatrix}\mathbf {i} &\mathbf {j} &\mathbf {k} \\{\frac {\partial }{\partial x}}&{\frac {\partial }{\partial y}}&{\frac {\partial }{\partial z}}\\A_{x}&A_{y}&A_{z}\end{vmatrix}}}
因為 |curl A · n| ≤ |curl A|, 所以旋度 curl A 可以詮釋為單位面積向量場 A 之最大環流。
因為線積分與座標無關, 所以旋度與所選取之座標無關,
旋度 ⇐⇒ 環流(線積分) ⇐⇒ Stokes定理
Laplace 算子
{\displaystyle \Delta f=\nabla ^{2}f=\nabla \cdot \nabla f}
線積分
\int_C F(r) dr = \int_a^b F(r(t)) dr/dt dt
完整程式
import numpy as np
pi = np.pi
step = 0.001
# 函數 f 對變數 p[k] 的偏微分: df / dp[k]
def df(f, p, k):
p1 = p.copy()
p1[k] = p[k]+step
return (f(p1) - f(p)) / step
# 單變數微分
def df1(f, x):
return df(f, np.array([x]), 0)
# 梯度：函數 f 在點 p 上的梯度
def grad(f, p):
gp = p.copy()
for k in range(len(p)):
gp[k] = df(f, p, k)
return gp
# 散度：函數 f 在點 p 上的散度
# 限制：f 為 n 維向量函數 f(Rn)->Rn
def div(f, p):
r = 0
for i in range(p.size):
r += df(f, p, i)[i]
return r
# 旋度：函數 f 在點 p 上的旋度
# 限制：f 為 3 維向量函數 f(R3)->R3
def curl(f, p):
rx = df(f, p, 1)[2]-df(f, p, 2)[1]
ry = df(f, p, 2)[0]-df(f, p, 0)[2]
rz = df(f, p, 0)[1]-df(f, p, 1)[0]
return np.array([rx, ry, rz])
# 拉普拉斯算子 ∆ = div∇
def laplace(f, p):
return div(lambda p:grad(f, p), p)
def lintegrate(F, r, a, b, dt=step):
area = 0.
t = a
while (t < b):
p = r(t)
dr = df1(r, t)
# print('t={} p={} F(p)={}'.format(t, p, F(p)))
area += np.dot(F(p), dr)*dt
t += dt
return area
測試程式
梯度
from vfield import *
# 我們想找函數 f 的最低點
def f(p):
[x,y] = p
print('x,y=', x, y)
return x*x + y*y
p = np.array([1.0,3])
print('df(f, p, 0) = ', df(f, p, 0))
print('df(f, p, 1) = ', df(f, p, 1))
print('grad(f)=', grad(f, p))
執行結果
$ python3 grad1.py
x,y= 1.001 3.0
x,y= 1.0 3.0
df(f, p, 0) = 2.0009999999999195
x,y= 1.0 3.001
x,y= 1.0 3.0
df(f, p, 1) = 6.000999999999479
x,y= 1.001 3.0
x,y= 1.0 3.0
x,y= 1.0 3.001
x,y= 1.0 3.0
grad(f)= [2.001 6.001]
散度
from vfield import *
# f = (3xz, 2xy, -yzz)
# div(f) = (3z+2x-2yz)
def f(p):
x,y,z = p
return np.array([3*x*z, 2*x*y, -y*z*z])
p = np.array([1.0, 0, 3])
print('div(f,p)=', div(f,p)) # 3z+2x+0 = 3*3+2*1+0 = 11
執行結果
$ python3 div1.py
div(f,p)= 10.999999999998565
f = (3xz, 2xy, -yz^2)
div f = 3z+2x-2yz
在 (x,y,z)=(1,0,3) 時，答案應為 (33+21+0)=11, 非常接近！
旋度
from vfield import *
# f = (yz, 3zx, z)
# curl(f) = (-3x, y, 2x)
def f(p):
x,y,z = p
return np.array([y*z, 3*z*x, z])
p = np.array([1.0, 0, 3])
print('f=(yz, 3zx, z) p=', p)
print('curl(f,p)=', curl(f,p)) # [-3,0,6]
執行結果
$ python3 curl1.py
f=(yz, 3zx, z) p= [1. 0. 3.]
curl(f,p)= [-3. 0. 6.]
f = (yz, 3zx, z)
curl f = (-3x, y, 2z)
在 (1,0,3) 的旋度應為 （-3, 0, 6)
拉普拉斯算子
from vfield import *
from numpy import linalg as LA
def f(p):
return 1.0/LA.norm(p)
p = np.array([1.,2.,3.])
print(laplace(f, p))
執行結果
$ python3 laplace1.py
2.1029733510147253e-05
引力場的拉普拉斯算子為零，程式結果 0.000021 非常接近零。
線積分
from vfield import *
from math import *
def F(p):
x,y=p
return np.array([-y, -x*y])
def r(t):
return np.array([cos(t), sin(t)])
print("lintegrate(0,pi/2)=", lintegrate(F, r, 0., pi/2))
本範例參考 紙本 kreyszig 工程數學精華本（歐亞書局） 8ed (424 頁)，答案為 pi/4-1/3 = 0.4521。
$ python3 lintegrate1.py
lintegrate(0,pi/2)= 0.4521851778637378