Zenn
🟰

自動微分と数値積分を用いたグリーンの恒等式の検算

2025/01/09に公開
Julia
math
物理
tech

有限要素法などで用いられる公式

\iiint_V f \pmb{\nabla}^2 g ~ \mathrm{d}V + \iiint_V \pmb{\nabla} f \cdot \pmb{\nabla} g ~ \mathrm{d}V = \iint_{\partial V} f \pmb{\nabla} g \cdot \pmb{n} ~ \mathrm{d}S

とこれを特殊化した式をJulia言語の自動微分と数値積分を用いて検算します. 上式はWikipediaではグリーンの第一恒等式と呼ばれており, こちらのページではグリーンの公式と呼ばれています. また, 上式を含む事実は物理のかぎしっぽではグリーンの第一定理と呼ばれており, こちらのページではグリーンの定理と呼ばれています. この公式の特殊な場合は変分法における運動エネルギー項の計算で使われます.

パッケージ

次のパッケージを使用します.

パッケージの読み込み
# import Pkg; Pkg.add("ForwardDiff")
# import Pkg; Pkg.add("Zygote")
# import Pkg; Pkg.add("QuadGK")
# import Pkg; Pkg.add("Cubature")
# import Pkg; Pkg.add("Antique")
import ForwardDiff
import Zygote
import QuadGK
import Cubature
import Antique

1階微分と2階微分は次のように計算できます. 1変数ではForwardDiff.jlが最速です. 詳しくはこちらを参照してください.

1階微分と2階微分の例
import ForwardDiff
f(x) = 2*x^3
df(x) = ForwardDiff.derivative(f, x)
d²f(x) = ForwardDiff.derivative(df, x)
@show df(1.0)
@show d²f(1.0)
# df(1.0) = 6.0
# d²f(1.0) = 12.0

勾配は gradient, ラプラシアンは hessian のトレースを取ることで計算できます. 多変数のヘッシアンではZygote.jlが最速です. 詳しくはこちらを参照してください.

勾配とラプラシアンの例
import Zygote
f(x) = x[1]^3 + 2*x[2]^3
∇f(x) = first(Zygote.gradient(f, x))
∇²f(x) = sum(first(Zygote.diaghessian(f, x)))
@show ∇f([1.0, 1.0])
@show ∇²f([1.0, 1.0])
# ∇f([1.0, 1.0]) = [3.0, 6.0]
# ∇²f([1.0, 1.0]) = 18.0

１次元

グリーンの公式は１次元における部分積分の公式の多次元への拡張になっているとのことです.

\int_a^b f(x) g^{\prime \prime}(x) \mathrm{d}x = \left. f(x) g^{\prime}(x) \right|_a ^b - \int_a^b f^{\prime}(x) g^{\prime}(x) \mathrm{d}x

下記のように左辺と右辺が誤差の範囲で一致することが確かめられます.

1次元の場合
import ForwardDiff
import QuadGK

a = -2
b = 5

f(x) = 2*sin(3*x)
df(x) = ForwardDiff.derivative(f, x)
d²f(x) = ForwardDiff.derivative(df, x)

g(x) = 3*x^2 + 2*x + 1
dg(x) = ForwardDiff.derivative(g, x)
d²g(x) = ForwardDiff.derivative(dg, x)

first(QuadGK.quadgk(x -> f(x) * d²g(x), a, b)) |> println # 左辺
f(b)*dg(b) - f(a)*dg(a) - first(QuadGK.quadgk(x -> df(x) * dg(x), a, b)) |> println # 右辺

# 6.879432798036747
# 6.879432798036696

上記の式において (a,b) = (-\infty,\infty) および f(-\infty) = f(\infty) = 0 または g'(-\infty) = g'(\infty) = 0 などの条件から右辺の第1項が0になり

\int_{-\infty}^{\infty} f(x) \frac{\mathrm{d}^2}{{\mathrm{d}x}^2} g(x) ~\mathrm{d} x = - \int_{-\infty}^{\infty} \frac{\mathrm{d} f(x)}{\mathrm{d} x} \frac{\mathrm{d} g(x)}{\mathrm{d} x} ~\mathrm{d}x

が得られます. 束縛状態ではこの条件を満たすため, 変分法における運動エネルギー項の計算に使えます. 証明は分かりませんが, より一般の演算子に対する記述がL. D. ランダウ, E. M. リフシッツ 著, 好村 滋洋, 井上 健男 訳『量子力学 ランダウ＝リフシッツ物理学小教程』(筑摩書房, 2008) p.39にあります.

1次元の特別な場合
import ForwardDiff
import QuadGK

f(x) = exp(-x^2)
df(x) = ForwardDiff.derivative(f, x)
d²f(x) = ForwardDiff.derivative(df, x)

g(x) = exp(-(x-4)^2)
dg(x) = ForwardDiff.derivative(g, x)
d²g(x) = ForwardDiff.derivative(dg, x)

 first(QuadGK.quadgk(x -> f(x) * d²g(x), -Inf, Inf)) |> println # 左辺
-first(QuadGK.quadgk(x -> df(x) * dg(x), -Inf, Inf)) |> println # 右辺

# 0.006306600811352378
# 0.006306600811372876

一次元調和振動子の波動関数での例も計算しておきましょう.

1次元の調和振動子
import ForwardDiff
import QuadGK

ψ(x) = π^(-1//4)*exp(-x^2/2)
dψ(x) = ForwardDiff.derivative(ψ, x)
d²ψ(x) = ForwardDiff.derivative(dψ, x)

first(QuadGK.quadgk(x -> ψ(x)' * d²ψ(x), -Inf, Inf)) |> println # 左辺
first(QuadGK.quadgk(x -> -dψ(x) * dψ(x), -Inf, Inf)) |> println # 右辺

# -0.5000000000000135
# -0.4999999999999719

正しくエネルギーが計算できました. Antique.jlを使うと次のように書けます.

1次元の調和振動子
import ForwardDiff
import QuadGK
import Antique
HO = Antique.HarmonicOscillator(k=1.0, m=1.0, ℏ=1.0)

ψ(x) = Antique.ψ(HO, x, n=0)
dψ(x) = ForwardDiff.derivative(x -> ψ(x), x)
d²ψ(x) = ForwardDiff.derivative(x -> dψ(x), x)

first(QuadGK.quadgk(x -> ψ(x)' * d²ψ(x), -Inf, Inf)) |> println # 左辺
first(QuadGK.quadgk(x -> -dψ(x) * dψ(x), -Inf, Inf)) |> println # 右辺

# -0.5000000000000135
# -0.4999999999999719

２次元

杉浦 光夫, 清水 英男, 金子 晃, 岡本 和夫『解析演習』(東京大学出版会, 1989) 7.4 系（平面におけるグリーンの公式）

\iint_D f \Delta g ~ dx dy + \iint_D \pmb{\nabla} f \cdot \pmb{\nabla} g ~ dx dy = \int_{\partial D} f \pmb{\nabla} g \cdot \pmb{n} ~ ds
2次元の場合
import Zygote
import Cubature
import LinearAlgebra

f(x) = x[1]^3 + 2*x[2]^3
∇f(x) = first(Zygote.gradient(f, x))
∇²f(x) = sum(first(Zygote.diaghessian(f, x)))

g(x) = exp(-(x[1]-4)^2) + x[2]^3
∇g(x) = first(Zygote.gradient(g, x))
∇²g(x) = sum(first(Zygote.diaghessian(g, x)))

# 左辺
first(Cubature.hcubature(x -> f(x) * ∇²g(x), [-2, 0], [3, 5])) +
first(Cubature.hcubature(x -> LinearAlgebra.dot(∇f(x), ∇g(x)), [-2, 0], [3, 5])) |> println

# 右辺
👉 = first(QuadGK.quadgk(y -> LinearAlgebra.dot(f([ 3,y]) * ∇g([ 3,y]), [ 1, 0]),  0, 5))
👈 = first(QuadGK.quadgk(y -> LinearAlgebra.dot(f([-2,y]) * ∇g([-2,y]), [-1, 0]),  0, 5))
👆 = first(QuadGK.quadgk(x -> LinearAlgebra.dot(f([ x,5]) * ∇g([ x,5]), [ 0, 1]), -2, 3))
👇 = first(QuadGK.quadgk(x -> LinearAlgebra.dot(f([ x,0]) * ∇g([ x,0]), [ 0,-1]), -2, 3))
👉 + 👈 + 👆 + 👇 |> println

# 95298.00210507686
# 95298.00209984844

３次元

杉浦 光夫, 清水 英男, 金子 晃, 岡本 和夫『解析演習』(東京大学出版会, 1989) より 7.6 系（空間におけるグリーンの公式）

\iiint_V f \Delta g ~ dx dy dz + \iiint_V \pmb{\nabla} f \cdot \pmb{\nabla} g ~ dx dy dz = \iint_{\partial V} f \pmb{\nabla} g \cdot \pmb{n} ~ dS
3次元の場合
import Zygote
import Cubature
import LinearAlgebra

f(x) = x[1]^3 + 2*x[2]^3 + 3*x[3]^4
∇f(x) = first(Zygote.gradient(f, x))
∇²f(x) = sum(first(Zygote.diaghessian(f, x)))

g(x) = exp(-(x[1]-4)^2) + x[2]^3 + x[3]^3
∇g(x) = first(Zygote.gradient(g, x))
∇²g(x) = sum(first(Zygote.diaghessian(g, x)))

# 左辺
first(Cubature.hcubature(x -> f(x) * ∇²g(x), [0,0,0], [3,4,5])) +
first(Cubature.hcubature(x -> LinearAlgebra.dot(∇f(x), ∇g(x)), [0,0,0], [3,4,5])) |> println

# 右辺
⚀ = first(Cubature.hcubature(x -> LinearAlgebra.dot(f([3,x[1],x[2]]) * ∇g([3,x[1],x[2]]), [ 1, 0, 0]), [0,0], [4,5]))
⚅ = first(Cubature.hcubature(x -> LinearAlgebra.dot(f([0,x[1],x[2]]) * ∇g([0,x[1],x[2]]), [-1, 0, 0]), [0,0], [4,5]))
⚁ = first(Cubature.hcubature(x -> LinearAlgebra.dot(f([x[1],4,x[2]]) * ∇g([x[1],4,x[2]]), [ 0, 1, 0]), [0,0], [3,5]))
⚄ = first(Cubature.hcubature(x -> LinearAlgebra.dot(f([x[1],0,x[2]]) * ∇g([x[1],0,x[2]]), [ 0,-1, 0]), [0,0], [3,5]))
⚂ = first(Cubature.hcubature(x -> LinearAlgebra.dot(f([x[1],x[2],5]) * ∇g([x[1],x[2],5]), [ 0, 0, 1]), [0,0], [3,4]))
⚃ = first(Cubature.hcubature(x -> LinearAlgebra.dot(f([x[1],x[2],0]) * ∇g([x[1],x[2],0]), [ 0, 0,-1]), [0,0], [3,4]))
⚀ + ⚁ + ⚂ + ⚃ + ⚄ + ⚅ |> println

# 2.0957813798088743e6
# 2.0957813797704456e6

右辺の積分を直方体の6面に場合分けする作業が手間でした. 上記の特別な場合である

\iiint_{\mathbb{R}^3} f^*(\pmb{r}) \nabla^2 g(\pmb{r}) ~\mathrm{d}\pmb{r} = - \iiint_{\mathbb{R}^3} \nabla f^*(\pmb{r}) \cdot \nabla g(\pmb{r}) ~\mathrm{d}\pmb{r}

も検算しておきます. 束縛状態の波動関数は無限遠で0であり勾配も0なので, 上記のグリーンの公式の右辺は0になります. とはいえ無限遠まで数値的に積分することはできませんので, 適当にカットオフを加えましょう.

3次元の特別な場合
import Zygote
import Cubature
import LinearAlgebra

f(x) = exp(-x[1]^2) * exp(-2*x[2]^2) * exp(-3*x[3]^2)
∇f(x) = first(Zygote.gradient(f, x))
∇²f(x) = sum(first(Zygote.diaghessian(f, x)))

g(x) = exp(-(x[1]-1)^2) * exp(-2*x[2]^2) * exp(-3*x[3]^2)
∇g(x) = first(Zygote.gradient(g, x))
∇²g(x) = sum(first(Zygote.diaghessian(g, x)))

 first(Cubature.hcubature(x -> f(x) * ∇²g(x), [-10,-10,-10], [10,10,10])) |> println
-first(Cubature.hcubature(x -> LinearAlgebra.dot(∇f(x), ∇g(x)), [-10,-10,-10], [10,10,10])) |> println

# -2.437400824735525
# -2.4374008246263603

まとめ

グリーンの第一恒等式を1, 2, 3次元の場合について数値的に検算できました. 同様の手法を使って色々な公式が検算できます.

バージョン情報

バージョン情報
Julia v1.11.2
ForwardDiff v0.10.38
Zygote v0.6.75
QuadGK v2.11.1
Cubature v1.5.1
Antique v0.11.2

https://gist.github.com/ohno/5e1954431bd9de9b2d3b783c429a10fe

Discussion