LoginSignup
15
7

More than 1 year has passed since last update.

数値計算Advent Calendar 2020Day 2
@sharow

抵抗格子問題の数値計算による近似解

Last updated at Posted at 2020-12-02

無限に広がる1Ω抵抗の格子の抵抗値を問う問題、Nerd Snipingの代名詞になっている問題です。

xkcd/356 Creative Commons Attribution-NonCommercial 2.5 License.

Googleの入社試験で出たということで有名になったようですが、それ以前からある問題です。解析解で詳しくやっているのはこれです。重要なことは、私には無理ってことです。

でも数値計算による近似解なら私にもできそうです。解析解で求められた一般解を使わずにです(使ったら負け)。

高度な数学はいりません。高度な回路理論もいりません。必要なのはオームの法則とキルヒホッフの法則、それから回路網を線型方程式にする方法、あとJuliaです。

オームの法則とキルヒホッフの法則

基本となる2つの法則をおさえてきましょう。オームの法則は2点間の電位差が$V=IR$を満たすというものです。$R$は電気抵抗で、単位はΩです。

この関係から、二つの値が分かれば、残りの値が計算できることがわかります。

  • $V=IR$
  • $I=\frac {V}{R}$
  • $R=\frac {V}{I}$

これが電圧、電流、抵抗の2点間の関係です。

一方キルヒホッフの法則は、2点間の関係ではなく、閉路と節点(ノード)に関する保存則です。電流則と電圧則があります。電流則は、あるノードに入る電流の総和と出て行く電流の総和が一致するという法則です。

\displaystyle \sum_{i=1}^N I_i=0

よく交差点に例えられます。交差点に入る車と出ていく車の数は一致する、という比喩です。それくらいの、「常識が通用する」くらいに捉えてかまいません。

電圧則はここでは使わないので説明しません、詳細はWikipediaなどみてください。電圧則は閉路に関する法則なので、今回使う回路解析においては閉路(ループ)ごとに式が出てくるので向いてません。

超簡単な例: 3ノード

簡単すぎて、慣れている人は暗算で10秒かからないでしょう。普通は合成抵抗を求めてオームの法則で終わりですが、ここではあえて電流則から線型方程式を組み立てて解いてみます。すごく面倒くさい無駄なことしてるように見えるでしょうけど、後にこの方法が数秒の暗算よりもよくスケールするという点をお見せできればと思います。

電源の電圧は1Vとします。

a0_1.png

ノードAの電圧を$V_A$、Bの電圧を$V_B$、同$V_C$としましょう。

それらのノードに関して電流則を書きます。その際、それぞれのノードに対して、そこから電流が湧き出していると仮定します。もちろん、実際に湧き出したりはしません、言い換えると $\sum湧き出し=0$ です(電流則)。1

waki.png

俺から湧き出している方向に流れる電流は、湧き出してった先とその間の抵抗に生じる電位差から求められます。

俺からおまえに流れる電流(矢印)= \frac{俺-おまえ}{R}

それから電圧源があれば電流が逆向きに流れると仮定して、湧き出している方向と同じであれば足せばよいです。逆方向であれば引きます。それらすべてを足すと、電流則によりゼロになります。つまり、

  • ノードA: $\frac{V_A-V_B}{1}+\frac{V_A-V_C}{1}+I_s= 0$
  • ノードB: $\frac{V_B-V_A}{1}+\frac{V_B-V_C}{1}=0$
  • ノードC: $\frac{V_C-V_B}{1}+\frac{V_C-V_A}{1}-I_s=0$

分母がすべて1なのは、すべて1Ωだからです。
$I_s$は電圧源の +端子から-端子に向かって通る電流です(図参照)。方向が直観とは逆である点に注意してください。その理由は、2点間の関係は$I=\frac{V_{drop}}{R}$、言い換えれば$\frac{V_{drop}}{R}-I=0$を満たします。そのうえで電流則も満たすので$\sum \frac{V_{drop}}{R} - \sum I = 0$です。どちらかの符号を逆にしないと成立しません。電圧則でいうと、電圧源の電圧降下が逆の符号であることと同じです2。
ざっくり言えば、電圧源を通る不明な$I_?$は+端子から-端子への電流を前提とする2、ということです。デメリットは、電流の解が逆の符号で出てくるだけです。

全部1Ωなので、$\frac{V_A-V_B}{1} = V_A - V_B$にしちゃいましょう。それから電源の電位、$V_A-V_C = 1$を加えました。

\begin{cases}
 V_A - V_B + V_A - V_C + I_s = 0 \\
 V_B-V_A + V_B - V_C = 0 \\
 V_C-V_B + V_C - V_A - I_s = 0 \\
 V_A-V_C = 1
\end{cases}

これを行列で表せば、

\begin{pmatrix}
 2 & -1 & -1 & 1 \\
 -1 & 2 & -1 & 0 \\
 -1 & -1 & 2 & -1 \\
 1 & 0 & -1 & 0
\end{pmatrix}
\begin{pmatrix}
 V_A \\
 V_B \\
 V_C \\
 I_s
\end{pmatrix} =
\begin{pmatrix}
 0 \\
 0 \\
 0 \\
 1
\end{pmatrix}

ここまでくると左の行列について逆行列を求める一般的な線型方程式になります。
Wolfram先生に訊いてみました。(※3)

Screenshot_2020-11-30 inverse-Wolfram Alpha-fs8.png
ということで、

\begin{pmatrix}
 V_A \\
 V_B \\
 V_C \\
 I_s
\end{pmatrix} =
\frac{1}{18}
\begin{pmatrix}
 1 & -2 & 1 & 9 \\
 -2 & 4 & -2 & 0 \\
 1 & -2 & 1 & -9 \\
 -9 & 0 & 9 & -27
\end{pmatrix}
\begin{pmatrix}
 0 \\
 0 \\
 0 \\
 1
\end{pmatrix}

より、

\begin{cases}
 V_A = \frac{9}{18} = 0.5 [V] \\
 V_B = \frac{0}{18} = 0 [V] \\
 V_C = -\frac{9}{18} = -0.5 [V] \\
 I_s = -\frac{27}{18} = -\frac{3}{2} = -1.5 [A]
\end{cases}

ここでは$V_C$をリファレンス(基準点、図では薄いグランド記号)として考えてますから、全てのノードから$V_C$を引きます($I_s$からは引いちゃダメです)。

\begin{cases}
 V_A = 1 [V] \\
 V_B = 0.5 [V] \\
 V_C = 0 [V] \\
 I_s = -1.5 [A]
\end{cases}

電圧源からは+端子から-端子に向かって-1.5[A]・・言い換えれば、-端子から+端子に向かって1.5[A]流れます。1[V]の電圧源から1.5[A]流れるのですから、$R_{total}=\frac{1}{1.5}=0.6666...$[Ω]になります。

オームの法則から言っても、2[Ω]と1[Ω]の並列接続ですから、
$\frac{1}{\frac{1}{1}+\frac{1}{2}}=\frac{2}{3}=0.6666...$[Ω]となります。

補足。ノードCはグランドなので削除することができます。先の係数行列Aは

\begin{pmatrix}
 2 & -1 & 1 \\
 -1 & 2 & 0 \\
 1 & 0 &  0
\end{pmatrix}

となり、(擬似でない)逆行列で求めることができます。電位の値はグランドとの電位になります。

>>> a
array([[ 2, -1,  1],
       [-1,  2,  0],
       [ 1,  0,  0]])
>>> np.linalg.solve(a, [0,0,1])
array([ 1. ,  0.5, -1.5])

逆にいうとこの方法はフローティングな回路を解析できません。・・まぁできないこともないのですが、係数行列が正則にならないのでコストが無駄に高くなります。

3x3ノード

まだ人間にも解けます。テブナンの定理とか使えば。

a1.png

先の簡単な問題とスケールが違うだけで、本質的には何も変わりません。単にノードが多いだけです。ここでは同じ説明を繰り返すのではなく、簡略化したModified Nodal Analysis(MNA)について一般的な点をおさえつつ解いていきます。

先の例でみたように、抵抗・電圧源のネットワークは線型方程式の形で以下のように表せます。

\mathbf {Ax} = \mathbf b

回路網解析の文脈ではこれらを
係数行列(coefficient matrix): $\mathbf A$
不明/変数行列(unknown matrix/vector): $\mathbf x$
励起行列(excitation matrix/vector): $\mathbf b$

ここで、$\mathbf A$行列は$\mathbf G$と$\mathbf i$から成ります。

\mathbf A = 
\begin{bmatrix}
\begin{array}{c|c}
\mathbf G & \mathbf i \\
\hline
\mathbf i^{\mathrm{T}} & \mathbf 0
\end{array}
\end{bmatrix}

先の例で見てみます。

G.png

$\mathbf G$はコンダクタンス行列といい、値の単位はΩの逆数であるジーメンス[S]です。なので、たとえば100Ωの値は$\frac{1}{100}$として扱わなければなりません。$\frac{V}{R}=I$はコンダクタンスで表すと$GV=I$になります(係数行列内の値は、当たり前ですが、係数(掛け算)でなければなりませんで)。今回は全部1[Ω]なので、その逆数も1[S]になり、書くのが面倒なので省略しました。
$\mathbf i$は電流の向きを表す行列で、(電圧源と抵抗しか無いから)今回は1, 0, -1しかとらず、右の行列と下の行列は転置したものになります。電圧源の数によって大きさが変わります。残りの隙間はゼロです。

$\mathbf x$には不明なノードの電圧と電圧源の電流が入ります。不明なので変数です。
$\mathbf b$には電圧源と電流源が入ります(今回電流源は使ってないのでそこはみんなゼロです)。

これらを準備したうえで、計算機で

\mathbf {x} = \mathbf {A}^{-1}\mathbf b

を求めれば一発でDC動作点が求まります。解が求まれば逆行列でなくともかまいません。消去法でもOKです。

では3x3の回路網にとりかかりましょう。今回も全部1[Ω]（=1[S]）なのでRは省略します。

\begin{cases}
 V_A-V_B + V_A-V_D = 0 \\
 V_B-V_A + V_B-V_C + V_B-V_E + I_s = 0 \\
 V_C-V_B + V_C-V_F = 0 \\
 V_D-V_A + V_D-V_E + V_D-V_G = 0 \\
 V_E-V_B + V_E-V_D + V_E-V_F + V_E-V_H = 0 \\
 V_F-V_C + V_F-V_E + V_F-V_I = 0 \\
 V_G-V_D + V_G-V_H - I_s = 0 \\
 V_H-V_G + V_H-V_E + V_H-V_I = 0 \\
 V_I-V_F + V_I-V_H = 0 \\
 V_B-V_G = 1
\end{cases}

しどいなこれ。

\begin{pmatrix}
 2 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -1 & 3 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & 1 \\
 0 & -1 & 2 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
 -1 & 0 & 0 & 3 & -1 & 0 & -1 & 0 & 0 & 0 \\
 0 & -1 & 0 & -1 & 4 & -1 & 0 & -1 & 0 & 0 \\
 0 & 0 & -1 & 0 & -1 & 3 & 0 & 0 & -1 & 0 \\
 0 & 0 & 0 & -1 & 0 & 0 & 2 & -1 & 0 & -1 \\
 0 & 0 & 0 & 0 & -1 & 0 & -1 & 3 & -1 & 0 \\
 0 & 0 & 0 & 0 & 0 & -1 & 0 & -1 & 2 & 0 \\
 0 & 1 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
\end{pmatrix}
\begin{pmatrix}
 V_A \\
 V_B \\
 V_C \\
 V_D \\
 V_E \\
 V_F \\
 V_G \\
 V_H \\
 V_I \\
 I_s
\end{pmatrix} =
\begin{pmatrix}
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 0 \\
 1
\end{pmatrix}

ノードが増えると疎行列になっていくことが分かります。

・・なぜかWolfram先生は行列が大きくなるとダメだったのでJuliaで。

julia> A
10×10 Array{Int64,2}:
  2  -1   0  -1   0   0   0   0   0   0
 -1   3  -1   0  -1   0   0   0   0   1
  0  -1   2   0   0  -1   0   0   0   0
 -1   0   0   3  -1   0  -1   0   0   0
  0  -1   0  -1   4  -1   0  -1   0   0
  0   0  -1   0  -1   3   0   0  -1   0
  0   0   0  -1   0   0   2  -1   0  -1
  0   0   0   0  -1   0  -1   3  -1   0
  0   0   0   0   0  -1   0  -1   2   0
  0   1   0   0   0   0  -1   0   0   0

julia> A \ [0 0 0 0 0 0 0 0 0 1]'
10×1 Array{Float64,2}:
 -0.05172413793103442
  0.22413793103448276
  0.05172413793103449
 -0.3275862068965516
 -0.15517241379310343
 -0.12068965517241378
 -0.775862068965517
 -0.39655172413793094
 -0.2586206896551724
 -0.8275862068965516

$V_G$(-0.77586...)を基準点(0[V])として計算すると、

\begin{cases}
 V_A = 0.7241 [V] \\
 V_B = 1.0 [V] \\
 V_C = 0.8275 [V] \\
 V_D = 0.4482 [V] \\
 V_E = 0.6206 [V] \\
 V_F = 0.6551 [V] \\
 V_G = 0 [V] \\
 V_H = 0.3793 [V] \\
 V_I = 0.5172 [V] \\
 I_s = 0.8275 [A]
\end{cases}

falstad.png

falstadでのシミュレーションと一致しています。

$V_B$、$V_G$間からみた抵抗$R_{total}$は$\frac{1}{0.82758621}=1.20833$[Ω]でした。

本題: 巨大な2次元抵抗格子

元々の2次元抵抗格子問題とは無限に広がる2次元の1Ω抵抗器の格子上で、ある点と、そこから(1, 2)離れた点の抵抗値を求めよ、という問題でした。数値計算では無限は無理ですが、10x10や100x100の解をみることで収束していく点を推測することは可能です。

inf.png

計算にあたり、コンダクタンス行列$\mathbf G$と抵抗ネットワークの関係について整理しておきましょう。いま4ノードの何も接続されていない回路網(？)があるとします。

lattice0.png

$\mathbf G$はノード数xノード数の大きさになります。$\mathbf G$の右上がAノードで、一個右がB, C、...とつづきます。縦軸もA, B, C、 ...です。

abcd.png

ここでA～B間に$R_1$を接続した場合を考えます。

lattice1.png

  • ノードA: $\frac{V_A-V_B}{R_1}=0 \Longleftrightarrow \frac{V_A}{R_1}-\frac{V_B}{R_1}=0$
  • ノードB: $\frac{V_B-V_A}{R_1}=0 \Longleftrightarrow \frac{V_B}{R_1}-\frac{V_A}{R_1}=0$

ですから、不明行列$\mathbf x$に移る$V_?$を忘れて係数だけ考えれば、

\mathbf G =
\begin{pmatrix}
 \frac{1}{R_1} & -\frac{1}{R_1} & 0 & 0 \\
 -\frac{1}{R_1} & \frac{1}{R_1} & 0 & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0
\end{pmatrix}

さらに$R_2$を繋げてみます

lattice2.png

  • ノードA: $\frac{V_A-V_B}{R_1}=0 \Longleftrightarrow \frac{V_A}{R_1}-\frac{V_B}{R_1}=0$
  • ノードB: $\frac{V_B-V_A}{R_1}+\frac{V_B-V_C}{R_2}=0 \Longleftrightarrow \frac{V_B}{R_1}-\frac{V_A}{R_1}+\frac{V_B}{R_2}-\frac{V_C}{R_2}=0$
  • ノードC: $\frac{V_C-V_B}{R_2}=0 \Longleftrightarrow \frac{V_C}{R_2}-\frac{V_B}{R_2}=0$
\mathbf G =
\begin{pmatrix}
 \frac{1}{R_1} & -\frac{1}{R_1} & 0 & 0 \\
 -\frac{1}{R_1} & \frac{1}{R_1} + \frac{1}{R_2} & -\frac{1}{R_2} & 0 \\
 0 & -\frac{1}{R_2} & \frac{1}{R_2} & 0 \\
 0 & 0 & 0 & 0
\end{pmatrix}

そうです、$\mathbf G$は重み(コンダクタンス)をもった無向グラフの隣接行列表現です。
あるノード間を抵抗$R_x$で接続することは、$\mathbf G$上でそのノード間に

\begin{bmatrix}
\begin{array}{c|c}
\frac{1}{R_x} & -\frac{1}{R_x} \\
\hline
-\frac{1}{R_x} & \frac{1}{R_x}
\end{array}
\end{bmatrix}

をそれぞれ足すことに等しいです。コンダクタンスで表せば、

\begin{bmatrix}
\begin{array}{c|c}
G_x & -G_x \\
\hline
-G_x & G_x
\end{array}
\end{bmatrix}

電流測の式など考えず、行列上にスタンプのようにポンと押せば(足せば)いいのです。
材料は全て揃いました。あとは計算するだけです。

汚れ仕事はJuliaにやらせよう

すいません、Juliaは初心者なのでベタ書きっぽくてアレかもです(汚したのは私です)。

説明で「解が求まれば逆行列でなくともよい」と言いました。計算途中で逆行列全部は必要ないのです。Juliaの \演算子3 はそれをうまいことやってくれます(たぶんLU分解です)。

julia version 1.5.3を使いました。

using LinearAlgebra
using SparseArrays

# 格子の一辺
lattice_size = 10
if size(ARGS)[1] >= 1
    lattice_size = parse(Int, ARGS[1])
end
println("lattice_size: ", lattice_size)

nodes = lattice_size * lattice_size
connected = Set()
G = spzeros(nodes, nodes)  # コンダクタンス行列

function connect(n1, n2)
    if isnothing(n1) || isnothing(n2) || n1 == n2 || (n1, n2) in connected
        return
    end
    G[n1, n1] += 1
    G[n2, n2] += 1
    G[n1, n2] -= 1
    G[n2, n1] -= 1
    union!(connected, [(n1, n2), (n2, n1)])
end

function lattice2index(y1, x1, y2, x2)
    outside = n -> (n < 1 || n > lattice_size)
    if outside(x1) || outside(x2) || outside(y1) || outside(y2)
        return (nothing, nothing)
    else
        return ((y1-1) * lattice_size + x1, (y2-1) * lattice_size + x2)
    end
end

for y = 1:lattice_size, x = 1:lattice_size
    # 上下左右の隣ノードに接続
    connect(lattice2index(y, x, y, x+1)...)
    connect(lattice2index(y, x, y, x-1)...)
    connect(lattice2index(y, x, y+1, x)...)
    connect(lattice2index(y, x, y-1, x)...)
end

# 電圧源の配置
if lattice_size <= 5
    src_idx = (3, 1)  # 狭すぎる場合
else
    t = Int(floor(lattice_size / 2))  # 真ん中らへんに配置
    src_idx = (t, t)
end
I = spzeros(nodes, 1)  # 電流行列
src = (src_idx[1]-1) * lattice_size + src_idx[2]
I[src] = -1  # ここから
I[src - (2*lattice_size) + 1] = 1 # ここへ接続

# 係数行列
A = sparse([
    G I
    I' [0]
])
b = zeros(nodes+1, 1)  # 励起行列
b[end] = 1  # 電圧源 1V

x = A \ b

src_current = -x[end]  # 電圧源から出た電流
total_resistance = 1 / src_current
expected = 4/pi - 1/2
println("resistor_count: ", length(connected)/2)
println("total_resistance: ", total_resistance)
println("expected: ", expected)
println("error: ", max(total_resistance, expected) - min(total_resistance, expected))

答えは解析解で$\frac{4}{\pi} - \frac{1}{2}$ ≒ $0.7732395447351628...$Ωになることが分かっています。

20x20あたりから結構いい値になってきます。

lattice_size: 10
resistor_count: 180.0
total_resistance: 0.805318134007435
expected: 0.7732395447351628
error: 0.03207858927227225

lattice_size: 20
resistor_count: 760.0
total_resistance: 0.7803087426450676
expected: 0.7732395447351628
error: 0.00706919790990479

lattice_size: 30
resistor_count: 1740.0
total_resistance: 0.776322664496643
expected: 0.7732395447351628
error: 0.0030831197614802353

lattice_size: 50
resistor_count: 4900.0
total_resistance: 0.7743392987522498
expected: 0.7732395447351628
error: 0.0010997540170870623

格子サイズ100～500。

lattice_size: 100
resistor_count: 19800.0
total_resistance: 0.773513442503032
expected: 0.7732395447351628
error: 0.0002738977678692356

lattice_size: 200
resistor_count: 79600.0
total_resistance: 0.7733079548515003
expected: 0.7732395447351628
error: 6.841011633751393e-5

lattice_size: 400
resistor_count: 319200.0
total_resistance: 0.7732566432550614
expected: 0.7732395447351628
error: 1.7098519898617326e-5

lattice_size: 500
resistor_count: 499000.0
total_resistance: 0.7732504874801841
expected: 0.7732395447351628
error: 1.0942745021336187e-5

1000x1000。抵抗器200万個。

$ /usr/bin/time julia solver.jl 1000
lattice_size: 1000
resistor_count: 1.998e6
total_resistance: 0.7732422803188629
expected: 0.7732395447351628
error: 2.735583700119726e-6
635.72user 2.98system 10:29.80elapsed 101%CPU (0avgtext+0avgdata 3486580maxresident)k
280inputs+0outputs (2major+1506650minor)pagefaults 0swaps

解析解との差が2.74μΩ。この辺でやめておきました。

この規模($\mathbf G$が100万x100万)の行列の解を普通のPCで、現実的な時間・現実的なメモリ空間で解ける理由はJuliaのSparseArraysの恩恵が大きいと思います。標準機能でSparse Arrayのまま\演算子が使えるのはチートではないでしょうか。

おまけ1

格子で特に斜め方向の考慮もしてないのでダイヤ状に広がるのかと思ってました。ダイヤ状になるのは電圧源のごく周辺だけのようです？ ちょっと意外。

図は電位で、電圧源を接続する距離を少し離した場合です。

cont.png

heatmap.png

heatmap2.png

電圧源の片側を動かしてみました。

anim.gif

こうしてみると静電場に似ていますね。

おまけ2

中心からの抵抗値で図にしてみました。

heatmap.png

surface.png

Link

参考文献

  • Leon O Chua, Pen-Min Lin "COMPUTER-AIDED ANALYSIS OF ELECTRONIC CIRCUITS" Prentice-Hall 1975

  1. 電圧則を使うと「全てのループを右回転で一巡した時云々」という仮定になります.  ↩

  2. Feynman Lectures On Physics / AC Circuits "we must remember that the voltage drop across a generator is the negative of the emf in that direction." ↩

  3. How to use the backslash operator in Julia? ↩

15
7
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
15
7