●数で遊ぼう
今回は簡単な数理パズルを出題します。Python でも簡単にプログラムを作ることができますが、できるだけ NumPy を使ってみてください。
●多角数
点を多角形の形に並べたとき、その総数を「多角数 (polygonal number)」といいます。詳しい説明は拙作のページ Puzzle DE Programming 多角数 をお読みください。
- 三角数と四角数を求めるプログラムを作ってください。
- 三角錐数と四角錐数を求めるプログラムを作ってください。
- p 角数とを p 角錐数求めるプログラムを作ってください。
- 三角数かつ四角数である数を「平方三角数 (square triangular number)」といいます。1,500,000 以下の平方三角数をすべて求めてください。
解答
●フィボナッチ数
フィボナッチ (fibonacci) 数はイタリアの数学者レオナルド・フィボナッチにちなんで名付けられた数です。詳しい説明は拙作のページ Puzzle DE Programming フィボナッチ数 をお読みください。
- n 以下のフィボナッチ数列を求めるプログラムを作ってください。
- 300,000,000 以下で最も大きいフィボナッチ数 Fn を求めてください。
- 初項 F0 から n 番目までのフィボナッチ数の和 F0 + F1 + F2 + ... + Fn を考えます。300,000,000 以下のフィボナッチ数の総和を求めてください。
- 300,000,000 以下のフィボナッチ数で、偶数になる項の総和を求めてください。
- 300,000,000 を 1 つ以上の連続しない相異なるフィボナッチ数の和として表してください。たとえば、7 = 1 + 1 + 2 + 3 と表すことができますが、1 を 2 回使っていることと、1, 2, 3 は連続したフィボナッチ数 (F2, F3, F4) なので条件を満たしていません。7 = 2 + 5 とすると条件を満たします。
解答
●素数
「素数 (prime number)」は正の約数が 1 と自分自身しかない自然数のことです。1 は素数に含めません。詳しい説明は拙作のページ Puzzle DE Programming 素数 をお読みください。
- 3,000,000 以下の素数の個数とその最大値を求めてください。
- 3,000,001 を素因数分解してください。素因数分解とは、素数でない整数 (合成数) を素数の積の形に書き表すことです。たとえば、12 は \(2^2 \times 3\) と素因数分解することができます。
- 3,000,000 以下のフィボナッチ素数 (フィボナッチ数で素数) の個数とその最大値を求めてください。
- 3,000,000 以下の素数で、差の最大値とその素数の組を求めてください。
- 差が 2 である素数の組を「双子素数 (twin prime)」といいます。3,000,000 以下の双子素数の個数とその最大値を求めてください。
解答
●約数
約数 (divisor) は整数 n を割り切る整数のことです。詳しい説明は拙作のページ Puzzle DE Programming 約数 をお読みください。
- 自然数 n の約数の個数を求めるプログラムを作ってください。
- 自然数 n の約数の和を求めるプログラムを作ってください。
- 自然数 n の約数をすべて求めるプログラムを作ってください。
解答
●ハミング数
7 以上の素数で割り切れない正の整数を「ハミング数 (Hamming Numbers)」といいます。ハミング数は素因子が 2, 3, 5 しかない自然数のことで、素因数分解したとき \(2^i \times 3^j \times 5^k \ (i, j, k \geq 0)\) の形式になります。たとえば、100 以下のハミング数は次のようになります。
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50,
54, 60, 64, 72, 75, 80, 81, 90, 96, 100
N 以下のハミング数を求めるプログラムを作ってください。
解答
●魔方陣
3 行 3 列の魔方陣を解くプログラムを作ってください。
上図の A から I の場所に 1 から 9 までの数字をひとつずつ配置します。縦横斜めの合計が等しくなるように数字を配置してください。
解答
●多角数の解答
リスト : 多角数の解答
import numpy as np
# 三角数と四角数
def triangular_number(n):
return np.arange(1, n + 1, dtype=np.int64).cumsum()
def square_number(n):
xs = np.arange(1, n + 1, dtype=np.int64)
return xs ** 2
# 三角錐数と四角錐数
def triangular_pyramidal(n):
return triangular_number(n).cumsum()
def square_pyramidal(n):
return square_number(n).cumsum()
# 多角数と多角錐数
def polygonal_number(p, n):
xs = np.arange(1, n + 1, dtype=np.int64)
return ((p - 2) * xs ** 2 - (p - 4) * xs) // 2
def polygonal_pyramidal(p, n):
return polygonal_number(p, n).cumsum()
# 平方三角数
def square_triangular_number(n):
r = np.roots(np.poly1d([1, 1, -2 * n]))
x = np.floor(r[1])
y = np.floor(np.sqrt(n))
return np.intersect1d(triangular_number(x), square_number(y))
print(triangular_number(20))
print(square_number(20))
print(triangular_pyramidal(20))
print(square_pyramidal(20))
print(polygonal_number(5, 20))
print(polygonal_pyramidal(5, 20))
print(polygonal_number(6, 20))
print(polygonal_pyramidal(6, 20))
print(square_triangular_number(1500000))
[ 1 3 6 10 15 21 28 36 45 55 66 78 91 105 120 136 153 171 190 210] # 三角数
[ 1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256 289 324 361 400] # 四角数
[ 1 4 10 20 35 56 84 120 165 220 286 364 455 560 680 # 三角錐数
816 969 1140 1330 1540]
[ 1 5 14 30 55 91 140 204 285 385 506 650 819 1015 # 四角錐数
1240 1496 1785 2109 2470 2870]
[ 1 5 12 22 35 51 70 92 117 145 176 210 247 287 330 376 425 477 532 590] # 五角数
[ 1 6 18 40 75 126 196 288 405 550 726 936 1183 1470 # 五角錐数
1800 2176 2601 3078 3610 4200]
[ 1 6 15 28 45 66 91 120 153 190 231 276 325 378 435 496 561 630 703 780] # 六角数
[ 1 7 22 50 95 161 252 372 525 715 946 1222 1547 1925 # 六角錐数
2360 2856 3417 4047 4750 5530]
[ 1 36 1225 41616 1413721] # 平方三角数
●フィボナッチ数の解答
リスト : フィボナッチ数の解答
import numpy as np
# フィボナッチ数
def fibo(n):
if n == 0 or n == 1:
return n
else:
f = np.matrix([[1, 1], [1, 0]], 'int64')
return (f ** (n - 1))[0, 0]
# n 以下のフィボナッチ数列
def fibo_seq(n):
fn = np.frompyfunc(fibo, 1, 1)
phi = (1 + np.sqrt(5)) / 2
x = np.floor(np.log(n * np.sqrt(5)) / np.log(phi))
return fn(np.arange(x + 1, dtype='int64'))
# n 以下で最大のフィボナッチ数
def fibo_max(n):
return fibo_seq(n)[-1]
# n 以下のフィボナッチ数の総和
def fibo_sum(n):
return np.sum(fibo_seq(n))
# n 以下のフィボナッチ数の偶数の総和
def fibo_even_sum(n):
xs = fibo_seq(n)
return np.sum(xs[xs % 2 == 0])
# 整数を 1 つ以上の連続しない相異なるフィボナッチ数の和で表す
# ゼッケンドルフの表現
def solver(n):
xs = []
ys = fibo_seq(n + 1)
while n > 0:
m = ys[-1]
xs.append(m)
n -= m
ys = ys[ys <= n]
return xs
print(fibo_seq(30000))
print(fibo_sum(300000000))
print(fibo_even_sum(300000000))
print(solver(300000000))
print(solver(267914296))
print(solver(267914296 - 1))
[0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657]
701408732
350704366
[267914296, 24157817, 5702887, 2178309, 46368, 233, 89, 1]
[267914296]
[165580141, 63245986, 24157817, 9227465, 3524578, 1346269, 514229, 196418, 75025, 28657,
10946, 4181, 1597, 610, 233, 89, 34, 13, 5, 2]
●素数の解答
リスト : 素数の解答
import numpy as np
# 素数
def sieve(n):
ps = np.ones((n - 3) // 2 + 1, 'u1')
x = 3
while x * x <= n:
y = (x - 3) // 2
if ps[y]: ps[y + x::x] = 0
x += 2
xs = np.where(ps == 1)[0]
return np.insert(xs * 2 + 3, 0, 2)
# 素因数分解
def factorization(n):
xs = []
ys = []
for p in sieve(np.int64(np.sqrt(n))):
if n < p * p: break
c = 0
while n % p == 0:
c += 1
n //= p
if c > 0:
xs.append(p)
ys.append(c)
if n > 1:
xs.append(n)
ys.append(1)
return np.array(xs), np.array(ys, dtype=np.int64)
# フィボナッチ数
def fibo(n):
if n == 0 or n == 1:
return n
else:
f = np.matrix([[1, 1], [1, 0]], 'int64')
return (f ** (n - 1))[0, 0]
# n 以下のフィボナッチ数列
def fibo_seq(n):
fn = np.frompyfunc(fibo, 1, 1)
phi = (1 + np.sqrt(5)) / 2
x = np.floor(np.log(n * np.sqrt(5)) / np.log(phi))
return fn(np.arange(x + 1, dtype='int64'))
# フィボナッチ素数
def fibo_prime(n):
xs = fibo_seq(n)
ys = sieve(n)
return np.intersect1d(xs, ys)
# 素数の個数と最大値
xs = sieve(3000000)
print(len(xs), np.max(xs))
# 素因数分解
print(factorization(3000001))
# フィボナッチ素数
a = fibo_prime(3000001)
print(a)
print(len(a), a[-1])
# 素数の差
ys = np.diff(xs)
m = np.argmax(ys)
print(np.max(ys), xs[m], xs[m + 1])
print(np.histogram(ys, np.arange(0, 151, 10)))
# 双子素数
zs = xs[np.where(ys == 2)[0]]
print(len(zs), np.max(zs))
216816 2999999
(array([ 853, 3517]), array([1, 1]))
[2 3 5 13 89 233 1597 28657 514229]
9 514229
148 2010733 2010881
(array([91347, 74257, 28769, 13778, 5391, 1855, 853, 360, 121,
50, 21, 6, 4, 2, 1], dtype=int64), array([ 0, 10, 20, 30, 40, 50,
60, 70, 80, 90, 100, 110, 120, 130, 140, 150]))
20932 2999831
●約数の解答
リスト : 約数の解答
import numpy as np
# 素数
def sieve(n):
ps = np.ones((n - 3) // 2 + 1, 'u1')
x = 3
while x * x <= n:
y = (x - 3) // 2
if ps[y]: ps[y + x::x] = 0
x += 2
xs = np.where(ps == 1)[0]
return np.insert(xs * 2 + 3, 0, 2)
# 素因数分解
def factorization(n):
xs = []
ys = []
for p in sieve(np.int64(np.sqrt(n))):
if n < p * p: break
c = 0
while n % p == 0:
c += 1
n //= p
if c > 0:
xs.append(p)
ys.append(c)
if n > 1:
xs.append(n)
ys.append(1)
return np.array(xs), np.array(ys, dtype=np.int64)
# 約数の個数
def divisor_num(n):
xs = factorization(n)[1] + 1
return np.prod(xs)
# 約数の合計値
def divisor_sum(n):
xs, ys = factorization(n)
a = 1
for x, y in zip(xs, ys):
a *= np.sum(x ** np.arange(y + 1))
return a
# 約数
def divisor(n):
xs, ys = factorization(n)
zs = xs[0] ** np.arange(ys[0] + 1)
for x, y in zip(xs[1:], ys[1:]):
zs = np.outer(zs, x ** np.arange(y + 1)).flatten()
zs.sort()
return zs
print(divisor_num(24))
print(divisor_num(12345678))
print(divisor_num(123456789))
print(divisor_num(1234567890))
print(divisor_num(1111111111))
print(divisor_sum(24))
print(divisor_sum(12345678))
print(divisor_sum(123456789))
print(divisor_sum(1234567890))
print(divisor_sum(1111111111))
print(divisor(24))
print(divisor(12345678))
print(divisor(123456789))
print(divisor(1234567890))
print(divisor(1111111111))
8
24
12
48
16
60
27319968
178422816
3211610688
1246404096
[ 1 2 3 4 6 8 12 24]
[ 1 2 3 6 9 18 47 94
141 282 423 846 14593 29186 43779 87558
131337 262674 685871 1371742 2057613 4115226 6172839 12345678]
[ 1 3 9 3607 3803 10821 11409
32463 34227 13717421 41152263 123456789]
[ 1 2 3 5 6 9
10 15 18 30 45 90
3607 3803 7214 7606 10821 11409
18035 19015 21642 22818 32463 34227
36070 38030 54105 57045 64926 68454
108210 114090 162315 171135 324630 342270
13717421 27434842 41152263 68587105 82304526 123456789
137174210 205761315 246913578 411522630 617283945 1234567890]
[ 1 11 41 271 451 2981
9091 11111 100001 122221 372731 2463661
4100041 27100271 101010101 1111111111]
●ハミング数の解答
ハミング数を求める一番簡単な方法は 1 から n までの数列を生成し、その要素が 2, 3, 5 だけで割り切れるか試すことです。プログラムはとても簡単ですが、引数 n の値が大きくなると時間がかかるようになります。
n に比べてハミング数の個数は少ないようなので、式 \(2^i \times 3^j \times 5^k \ (i, j, k \geq 0)\) を使ってハミング数を生成したほうがよさそうです。この式の計算は NumPy の meshgrid() を使うと簡単です。
プログラムと結果を以下に示します。
リスト : ハミング数の解答
import numpy as np
def hamming(n):
def init(m):
xs = [1]
k = m
while k <= n:
xs.append(k)
k *= m
return np.array(xs, dtype=np.float64)
x, y, z = np.meshgrid(init(2), init(3), init(5))
a = (x * y * z).flat
return np.sort(a[a <= n])
print(hamming(100))
for x in [1000, 10000, 100000, 1000000, 10000000]:
print(x, len(hamming(x)))
[ 1. 2. 3. 4. 5. 6. 8. 9. 10. 12. 15. 16. 18. 20.
24. 25. 27. 30. 32. 36. 40. 45. 48. 50. 54. 60. 64. 72.
75. 80. 81. 90. 96. 100.]
1000 86
10000 175
100000 313
1000000 507
10000000 768
●魔方陣の解答
下記のプログラムは単純な生成検定法なので、実行時間はちょっと遅いです。
リスト : 魔方陣 (3 * 3 盤)
import numpy as np
import itertools
def magic_square():
for xs in itertools.permutations(range(1, 9 + 1)):
a = np.array(xs).reshape((3, 3))
n = a.trace()
if a[:, ::-1].trace() != n: continue
b = a.sum(axis=0)
if not np.all(b == n): continue
b = a.sum(axis=1)
if not np.all(b == n): continue
print(a)
magic_square()
[[2 7 6]
[9 5 1]
[4 3 8]]
[[2 9 4]
[7 5 3]
[6 1 8]]
[[4 3 8]
[9 5 1]
[2 7 6]]
[[4 9 2]
[3 5 7]
[8 1 6]]
[[6 1 8]
[7 5 3]
[2 9 4]]
[[6 7 2]
[1 5 9]
[8 3 4]]
[[8 1 6]
[3 5 7]
[4 9 2]]
[[8 3 4]
[1 5 9]
[6 7 2]]
●連立一次方程式の解法
今回は連立一次方程式を解く基本的なアルゴリズムを取り上げます。NumPy には連立一次方程式を解く関数 linalg.solve() が用意されているので、実用的にはそちらを使えばいいのですが、Python (NumPy) とアルゴリズムのお勉強ということで、あえてプログラムを作ってみましょう。
●ガウスの消去法
ガウスの消去法 (Gaussian elimination) は人が連立方程式を解くときの方法とほとんど同じです。次の図を見てください。
a1 * x + a2 * y + a3 * z = d1 a1 * x + a2 * y + a3 * z = d1 a1 * x + a2 * y + a3 * z = d1
b1 * x + b2 * y + b3 * z = d2 ==> b2' * y + b3' * z = d2' ==> b2' * y + b3' * z = d2'
c1 * x + c2 * y + c3 * z = d3 c2' * y + c3' * z = d3' c3'' * z = d3''
図 : 前進消去
1 番目の式を b1 / a1 倍して 2 番目の式から引き算すると、x の項を消去することができます。同様の方法で 3 番目の式の x の項を消去することができます。次に、2 番目の式を c2' / b2' 倍して 3 番目の式から引き算すると y の項を消去することができます。これを「前進消去」といいます。
前進消去を行うと、最後の式には変数が z しかありません。z の値は d3'' / c3'' に決定することができます。次はひとつ前の式に戻り、 変数 z に d3'' / c3'' を代入すると変数 y の値を求めることができます。これを繰り返して、先頭の式に戻ると変数 x の値を求めることができます。これを「後退代入」といいます。
連立一次方程式は行列を使って Ax = b と表すことができます。A を係数行列といい、A と b を連結した行列を拡大係数行列といいます。
[[a1, a2, a3], [[a1, a2, a3, d1],
[b1, b2, b3], [b1, b2, b3, d2],
[c1, c2, c3]] [c1, c2, c3, d3]]
係数行列 拡大係数行列
ガウスの消去法は、拡大係数行列を使うと簡単にプログラムを作ることができます。Python でプログラムすると、次のようになります。
リスト : ガウスの消去法
def gauss(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = [x + [y] for x, y in zip(xs, ys)]
# 前進消去
for i in range(n - 1):
for j in range(i + 1, n):
temp = zs[j][i] / zs[i][i]
# 係数を 0 にする計算は省略可
# その場合は range(i + 1, n + 1) とする
for k in range(i, n + 1):
zs[j][k] -= temp * zs[i][k]
# for debug
print(zs)
# 後退代入
# 解は zs の n 列目にセットする
for i in range(n - 1, -1, -1):
for j in range(i + 1, n): zs[i][n] -= zs[i][j] * zs[j][n]
zs[i][n] /= zs[i][i]
return [z[-1] for z in zs]
# x + y = 100
# 2x + 4y = 272
print(gauss([[1, 1], [2, 4]], [100, 272]))
# x + y + z = 10
# 2x + 4y + 6z = 38
# 2x + 4z = 14
print(gauss([[1, 1, 1], [2, 4, 6], [2, 0, 4]], [10, 38, 14]))
# a + b + c + d = -5
# -a + b - c + d = -7
# 8a + 4b + 2c + d = -31
# -8a + 4b + -2c + d = -35
print(gauss([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]], [-5, -7, -31, -35]))
# a - b + c - d + e = 1
# 12a - 6b + 2c = 0
# a + b + c + d + e = 8
# 12a + 6b + 2c = 0
# 4a + 3b + 2c + d = 1
print(gauss([[1,-1,1,-1,1], [12,-6,2,0,0], [1,1,1,1,1], [12,6,2,0,0], [4,3,2,1,0]], [1,0,8,0,1]))
関数 gauss() の引数 xs が係数行列、ys が右辺の定数を格納したリストです。最初に、拡大係数行列を生成して変数 zs にセットします。あとはアルゴリズムをそのままプログラムしただけなので、特に難しいところはないと思います。実行結果は次のようになります。なお、拡大行列の表示は手作業で整形しています。
[[1, 1, 100],
[0.0, 2.0, 72.0]]
[64.0, 36.0]
[[1, 1, 1, 10],
[0.0, 2.0, 4.0, 18.0],
[0.0, -2.0, 2.0, -6.0]]
[[1, 1, 1, 10],
[0.0, 2.0, 4.0, 18.0],
[0.0, 0.0, 6.0, 12.0]]
[3.0, 5.0, 2.0]
[[1, 1, 1, 1, -5],
[0.0, 2.0, 0.0, 2.0, -12.0],
[0.0, -4.0, -6.0, -7.0, 9.0],
[0.0, 12.0, 6.0, 9.0, -75.0]]
[[1, 1, 1, 1, -5],
[0.0, 2.0, 0.0, 2.0, -12.0],
[0.0, 0.0, -6.0, -3.0, -15.0],
[0.0, 0.0, 6.0, -3.0, -3.0]]
[[1, 1, 1, 1, -5],
[0.0, 2.0, 0.0, 2.0, -12.0],
[0.0, 0.0, -6.0, -3.0, -15.0],
[0.0, 0.0, 0.0, -6.0, -18.0]]
[0.0, -9.0, 1.0, 3.0]
[[1, -1, 1, -1, 1, 1],
[0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
[0.0, 2.0, 0.0, 2.0, 0.0, 7.0],
[0.0, 18.0, -10.0, 12.0, -12.0, -12.0],
[0.0, 7.0, -2.0, 5.0, -4.0, -3.0]]
[[1, -1, 1, -1, 1, 1],
[0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
[0.0, 0.0, 3.333333333333333, -2.0, 4.0, 11.0],
[0.0, 0.0, 20.0, -24.0, 24.0, 24.0],
[0.0, 0.0, 9.666666666666668, -9.0, 10.0, 11.0]]
[[1, -1, 1, -1, 1, 1],
[0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
[0.0, 0.0, 3.333333333333333, -2.0, 4.0, 11.0],
[0.0, 0.0, 0.0, -11.999999999999998, -3.552713678800501e-15, -42.000000000000014],
[0.0, 0.0, 0.0, -3.1999999999999984, -1.6000000000000032, -20.90000000000001]]
[[1, -1, 1, -1, 1, 1],
[0.0, 6.0, -10.0, 12.0, -12.0, -12.0],
[0.0, 0.0, 3.333333333333333, -2.0, 4.0, 11.0],
[0.0, 0.0, 0.0, -11.999999999999998, -3.552713678800501e-15, -42.000000000000014],
[0.0, 0.0, 0.0, -4.440892098500626e-16, -1.6000000000000023, -9.70000000000001]]
[0.3124999999999991, 0.0, -1.874999999999997, 3.5, 6.062499999999997]
NumPy を使ったプログラムと実行結果を以下に示します。
リスト : ガウスの消去法 (NumPy バージョン)
import numpy as np
def gauss_np(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
# 前進消去
for i in range(n - 1):
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i:] -= temp * zs[i, i:]
# for debug
print(zs)
# 後退代入
for i in range(n - 1, -1, -1):
zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
zs[i, n] /= zs[i, i]
return zs[:, n]
print(gauss_np(np.array([[1, 1], [2, 4]]), np.array([100, 272])))
print(gauss_np(np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]]), np.array([10, 38, 14])))
print(gauss_np(np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]]), np.array([-5, -7, -31, -35])))
print(gauss_np(np.array([[1,-1,1,-1,1], [12,-6,2,0,0], [1,1,1,1,1], [12,6,2,0,0], [4,3,2,1,0]]), np.array([1,0,8,0,1])))
[[ 1. 1. 100.]
[ 0. 2. 72.]]
[ 64. 36.]
[[ 1. 1. 1. 10.]
[ 0. 2. 4. 18.]
[ 0. -2. 2. -6.]]
[[ 1. 1. 1. 10.]
[ 0. 2. 4. 18.]
[ 0. 0. 6. 12.]]
[ 3. 5. 2.]
[[ 1. 1. 1. 1. -5.]
[ 0. 2. 0. 2. -12.]
[ 0. -4. -6. -7. 9.]
[ 0. 12. 6. 9. -75.]]
[[ 1. 1. 1. 1. -5.]
[ 0. 2. 0. 2. -12.]
[ 0. 0. -6. -3. -15.]
[ 0. 0. 6. -3. -3.]]
[[ 1. 1. 1. 1. -5.]
[ 0. 2. 0. 2. -12.]
[ 0. 0. -6. -3. -15.]
[ 0. 0. 0. -6. -18.]]
[ 0. -9. 1. 3.]
[[ 1. -1. 1. -1. 1. 1.]
[ 0. 6. -10. 12. -12. -12.]
[ 0. 2. 0. 2. 0. 7.]
[ 0. 18. -10. 12. -12. -12.]
[ 0. 7. -2. 5. -4. -3.]]
[[ 1. -1. 1. -1. 1. 1. ]
[ 0. 6. -10. 12. -12. -12. ]
[ 0. 0. 3.33333333 -2. 4. 11. ]
[ 0. 0. 20. -24. 24. 24. ]
[ 0. 0. 9.66666667 -9. 10. 11. ]]
[[ 1.00000000e+00 -1.00000000e+00 1.00000000e+00 -1.00000000e+00
1.00000000e+00 1.00000000e+00]
[ 0.00000000e+00 6.00000000e+00 -1.00000000e+01 1.20000000e+01
-1.20000000e+01 -1.20000000e+01]
[ 0.00000000e+00 0.00000000e+00 3.33333333e+00 -2.00000000e+00
4.00000000e+00 1.10000000e+01]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -1.20000000e+01
-3.55271368e-15 -4.20000000e+01]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -3.20000000e+00
-1.60000000e+00 -2.09000000e+01]]
[[ 1.00000000e+00 -1.00000000e+00 1.00000000e+00 -1.00000000e+00
1.00000000e+00 1.00000000e+00]
[ 0.00000000e+00 6.00000000e+00 -1.00000000e+01 1.20000000e+01
-1.20000000e+01 -1.20000000e+01]
[ 0.00000000e+00 0.00000000e+00 3.33333333e+00 -2.00000000e+00
4.00000000e+00 1.10000000e+01]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -1.20000000e+01
-3.55271368e-15 -4.20000000e+01]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -4.44089210e-16
-1.60000000e+00 -9.70000000e+00]]
[ 0.3125 0. -1.875 3.5 6.0625]
●ガウス・ジョルダン法
ガウスの消去法は、係数行列の部分を上三角行列に変換することで、連立一次方程式を解きました。ここで、係数行列を単位行列に変換できれば、後退代入することなく解を求めることができます。これを「ガウス・ジョルダン法 (Gauss - Jordan elimination)」といいます。
簡単な動作例を下図に示します。
次の連立方程式をガウス・ジョルダン法で解く
x + y + z = 10
2x + 4y + 6z = 38
2x + 4z = 14
拡大係数行列を zs とする
(1) zs[0, 0] を 1 にする
[[ 1. 1. 1. 10.], # zs[0, 0:] /= 1
[ 2. 4. 6. 38.],
[ 2. 0. 4. 14.]]
(2) zs[1, 0], zs[2, 0] を 0 にする
[[ 1. 1. 1. 10.]
[ 0. 2. 4. 18.] # zs[1, 0:] -= zs[0, 0:] * 2
[ 0. -2. 2. -6.]] # zs[2, 0:] -= zs[0, 0:] * 2
(3) zs[1, 1] を 1 にする
[[ 1. 1. 1. 10.]
[ 0. 1. 2. 9.] # zs[1, 1:] /= 2
[ 0. -2. 2. -6.]]
(4) zs[0, 1], zs[2, 1] を 0 にする
[[ 1. 0. -1. 1.] # zs[0, 1:] -= zs[1, 1:] * 1
[ 0. 1. 2. 9.]
[ 0. 0. 6. 12.]] # zs[2, 1:] -= zs[1, 1:] * -2
(5) zs[2, 2] を 1 にする
[[ 1. 0. -1. 1.]
[ 0. 1. 2. 9.]
[ 0. 0. 1. 2.]] # zs[2, 2:] /= 6
(6) zs[2, 0], zs[2, 1] を 0 にする
[[ 1. 0. 0. 3.] # zs[0, 2:] -= zs[2, 2:] * -1
[ 0. 1. 0. 5.] # zs[1, 2:] -= zs[2, 2:] * 2
[ 0. 0. 1. 2.]]
このように、ガウス・ジョルダン法の基本的な考え方は実にシンプルで、プログラムも簡単になるのですが、速度はガウスの消去法よりも少し遅くなるといわれています。以下にプログラムと実行例を示します。
リスト : ガウス・ジョルダン法
import numpy as np
def gauss_jordan(xs, ys):
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
for i in range(n):
# zs[i, i] を 1 にする
temp = zs[i, i]
zs[i, i:] /= temp
# zs[i] 以外の行の i 番目の係数を 0 にする
for j in range(n):
if i != j:
temp = zs[j, i]
zs[j, i:] -= zs[i, i:] * temp
# for debug
print(zs)
return zs[:, n]
print(gauss_jordan(np.array([[1, 1], [2, 4]]), np.array([100, 272])))
print(gauss_jordan(np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]]), np.array([10, 38, 14])))
print(gauss_jordan(np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]]), np.array([-5, -7, -31, -35])))
print(gauss_jordan(np.array([[1,-1,1,-1,1], [12,-6,2,0,0], [1,1,1,1,1], [12,6,2,0,0], [4,3,2,1,0]]), np.array([1,0,8,0,1])))
[[ 1. 1. 100.]
[ 0. 2. 72.]]
[[ 1. 0. 64.]
[ 0. 1. 36.]]
[ 64. 36.]
[[ 1. 1. 1. 10.]
[ 0. 2. 4. 18.]
[ 0. -2. 2. -6.]]
[[ 1. 0. -1. 1.]
[ 0. 1. 2. 9.]
[ 0. 0. 6. 12.]]
[[ 1. 0. 0. 3.]
[ 0. 1. 0. 5.]
[ 0. 0. 1. 2.]]
[ 3. 5. 2.]
[[ 1. 1. 1. 1. -5.]
[ 0. 2. 0. 2. -12.]
[ 0. -4. -6. -7. 9.]
[ 0. 12. 6. 9. -75.]]
[[ 1. 0. 1. 0. 1.]
[ 0. 1. 0. 1. -6.]
[ 0. 0. -6. -3. -15.]
[ 0. 0. 6. -3. -3.]]
[[ 1. 0. 0. -0.5 -1.5]
[ 0. 1. 0. 1. -6. ]
[ -0. -0. 1. 0.5 2.5]
[ 0. 0. 0. -6. -18. ]]
[[ 1. 0. 0. 0. 0.]
[ 0. 1. 0. 0. -9.]
[ 0. 0. 1. 0. 1.]
[-0. -0. -0. 1. 3.]]
[ 0. -9. 1. 3.]
[[ 1. -1. 1. -1. 1. 1.]
[ 0. 6. -10. 12. -12. -12.]
[ 0. 2. 0. 2. 0. 7.]
[ 0. 18. -10. 12. -12. -12.]
[ 0. 7. -2. 5. -4. -3.]]
[[ 1. 0. -0.66666667 1. -1. -1. ]
[ 0. 1. -1.66666667 2. -2. -2. ]
[ 0. 0. 3.33333333 -2. 4. 11. ]
[ 0. 0. 20. -24. 24. 24. ]
[ 0. 0. 9.66666667 -9. 10. 11. ]]
[[ 1. 0. 0. 0.6 -0.2 1.2]
[ 0. 1. 0. 1. 0. 3.5]
[ 0. 0. 1. -0.6 1.2 3.3]
[ 0. 0. 0. -12. 0. -42. ]
[ 0. 0. 0. -3.2 -1.6 -20.9]]
[[ 1. 0. 0. 0. -0.2 -0.9]
[ 0. 1. 0. 0. 0. 0. ]
[ 0. 0. 1. 0. 1.2 5.4]
[-0. -0. -0. 1. -0. 3.5]
[ 0. 0. 0. 0. -1.6 -9.7]]
[[ 1. 0. 0. 0. 0. 0.3125]
[ 0. 1. 0. 0. 0. 0. ]
[ 0. 0. 1. 0. 0. -1.875 ]
[-0. -0. -0. 1. 0. 3.5 ]
[-0. -0. -0. -0. 1. 6.0625]]
[ 0.3125 0. -1.875 3.5 6.0625]
●ピボット選択
ところで、ガウスの消去法とガウス・ジョルダン法は、拡大係数行列の対角成分 zs[i, i] が 0 になると計算できなくなります。たとえば、次の連立方程式は解くことができますが、zs[0, 0] が 0 になっているため、プログラムの実行中にエラー (divide by zero) が発生します。
2y + 4z = 14 [[0. 2. 4. 14.]
x + y + z = 10 [1. 1. 1. 10.]
4x + 2y + 6z = 38 [4. 2. 6. 38.]]
連立方程式 拡大係数行列 zs
このとき、拡大係数行列の行を交換すると連立一次方程式が解ける場合があります。また、zs[i, i] が 0 に等しくなくても 0 に近い値だと解の誤差が大きくなります。そこで、zs[i, i] の絶対値がなるべく大きくなるように行を交換します。これを「部分ピボット選択」といいます。
なお、参考文献 1 によると、『行だけでなく列も交換する完全ピボット選択という方法もあるが、通常は部分ピボット選択で十分である』 とのことです。
部分ピボット選択のプログラムは NumPy を使うと簡単です。次のリストを見てください。
リスト : 部分ピボット選択
import numpy as np
def select_pivot(xs, i):
k = np.abs(xs[i:,i]).argmax() + i
if k != i:
temp = xs[i].copy()
xs[i] = xs[k]
xs[k] = temp
# ガウスの消去法
def gauss1(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
# 前進消去
for i in range(n - 1):
# ピボット選択
select_pivot(zs, i)
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i:] -= temp * zs[i, i:]
# for debug
print(zs)
# 後退代入
for i in range(n - 1, -1, -1):
zs[i, n] -= zs[i, i+1:n] @ zs[i+1:, n]
zs[i, n] /= zs[i, i]
return zs[:, n]
# ガウス・ジョルダン法
def gauss_jordan1(xs, ys):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), ys.astype(np.float_)]
for i in range(n):
# ピボット選択
select_pivot(zs, i)
temp = zs[i, i]
zs[i, i:] /= temp
for j in range(n):
if i != j:
temp = zs[j, i]
zs[j, i:] -= zs[i, i:] * temp
# for debug
print(zs)
return zs[:, n]
print(gauss1(np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]]), np.array([14, 10, 38])))
print(gauss_jordan1(np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]]), np.array([14, 10, 38])))
ピボット選択は関数 select_pivot() で行います。np.abs() で xs[i:, i] の絶対値を求め、その中から最大値の位置を argmax() で求めます。実際の位置はそれに i を加算したものになります。これを変数 k にセットします。k と i が異なる値であれば i 行と k 行と交換します。
Python のリストであれば xs[i], xs[k] = xs[k], xs[i] で交換できるのですが、NumPy の配列は交換できません。NumPy のスライス操作はビューを生成するだけなので、copy() を使って行の値をコピーする必要があります。ご注意くださいませ。
あとは select_pivot() を呼び出すだけです。それでは実際に試してみましょう。
[[ 4. 2. 6. 38. ]
[ 0. 0.5 -0.5 0.5]
[ 0. 2. 4. 14. ]]
[[ 4. 2. 6. 38. ]
[ 0. 2. 4. 14. ]
[ 0. 0. -1.5 -3. ]]
[ 5. 3. 2.]
[[ 1. 0.5 1.5 9.5]
[ 0. 0.5 -0.5 0.5]
[ 0. 2. 4. 14. ]]
[[ 1. 0. 0.5 6. ]
[ 0. 1. 2. 7. ]
[ 0. 0. -1.5 -3. ]]
[[ 1. 0. 0. 5.]
[ 0. 1. 0. 3.]
[ 0. 0. 1. 2.]]
[ 5. 3. 2.]
正しく解を求めることができました。
●逆行列
連立方程式 \(Ax = b\) を解く場合、A の逆行列 \(A^{-1}\) がわかれば、\(x = A^{-1}b\) で解を求めることができます。実をいうと、逆行列はガウス・ジョルダン法を使って簡単に求めることができるのです。次の図を見てください。
逆行列の定義 : A @ A-1 = I
[[a11, a12, a13], [[b11, b12, b13], [[1, 0, 0],
[a21, a22, a23], @ [b21, b22, b23], = [0, 1, 0],
[a31, a32, a33]] [b31, b32, b33]] [0, 0, 1]]
係数行列 逆行列 単位行列 I
上記の定義により逆行列 A-1 は次の方程式を解けば求めることができる
A @ [b11, b21, b31] = [1, 0, 0]
A @ [b12, b22, b32] = [0, 1, 0]
A @ [b13, b23, b33] = [0, 0, 1]
これは次の拡大係数行列にガウス・ジョルダン法を適用すること同じ
[[a11, a12, a13, 1, 0, 0],
[a21, a22, a23, 0, 1, 0],
[a31, a32, a33, 0, 0, 1]]
拡大係数行列の右半分に逆行列が求まる
プログラムも簡単です。次のリストを見てください。
リスト : ガウス・ジョルダン法で逆行列を求める
import numpy as np
def select_pivot(xs, i):
k = np.abs(xs[i:,i]).argmax() + i
if k != i:
temp = xs[i].copy()
xs[i] = xs[k]
xs[k] = temp
def matinv_gj(xs):
# 拡大係数行列の生成
n = len(xs)
zs = np.c_[xs.astype(np.float_), np.eye(n)]
for i in range(n):
# ピボットの選択
select_pivot(zs, i)
temp = zs[i, i]
zs[i, i:] /= temp
for j in range(n):
if i != j:
temp = zs[j, i]
zs[j, i:] -= zs[i, i:] * temp
return zs[:, n:]
# 簡単なテスト
a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
a5 = np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]])
for a in [a1, a2, a3, a4, a5]:
b = matinv_gj(a)
print(b)
print(a @ b)
拡大係数行列を作るとき、右半分に単位行列 eye(n) を連結します。あとはガウス・ジョルダン法で解を求め、最後に拡大係数行列の右半分を返します。それでは実行してみましょう。
[[ 2. -0.5]
[-1. 0.5]]
[[1. 0.]
[0. 1.]]
[[ 1.33333333 -0.33333333 0.16666667]
[ 0.33333333 0.16666667 -0.33333333]
[-0.66666667 0.16666667 0.16666667]]
[[1.00000000e+00 2.77555756e-17 2.77555756e-17]
[2.22044605e-16 1.00000000e+00 5.55111512e-17]
[0.00000000e+00 0.00000000e+00 1.00000000e+00]]
[[-0.16666667 0.16666667 0.08333333 -0.08333333]
[-0.16666667 -0.16666667 0.16666667 0.16666667]
[ 0.66666667 -0.66666667 -0.08333333 0.08333333]
[ 0.66666667 0.66666667 -0.16666667 -0.16666667]]
[[ 1.00000000e+00 0.00000000e+00 0.00000000e+00 -2.77555756e-17]
[ 0.00000000e+00 1.00000000e+00 0.00000000e+00 2.77555756e-17]
[-1.11022302e-16 0.00000000e+00 1.00000000e+00 -1.11022302e-16]
[ 0.00000000e+00 -1.11022302e-16 -8.32667268e-17 1.00000000e+00]]
[[-6.25000000e-02 4.16666667e-02 6.25000000e-02 8.33333333e-02
-1.25000000e-01]
[ 0.00000000e+00 -8.33333333e-02 0.00000000e+00 8.33333333e-02
0.00000000e+00]
[ 3.75000000e-01 -3.46944695e-18 -3.75000000e-01 -2.50000000e-01
7.50000000e-01]
[-5.00000000e-01 8.33333333e-02 5.00000000e-01 -8.33333333e-02
-5.55111512e-17]
[ 1.87500000e-01 -4.16666667e-02 8.12500000e-01 1.66666667e-01
-6.25000000e-01]]
[[ 1.00000000e+00 2.08166817e-17 0.00000000e+00 2.77555756e-17
1.11022302e-16]
[ 1.11022302e-16 1.00000000e+00 -1.11022302e-16 1.11022302e-16
0.00000000e+00]
[ 0.00000000e+00 6.93889390e-18 1.00000000e+00 5.55111512e-17
1.11022302e-16]
[ 1.11022302e-16 2.08166817e-17 -1.11022302e-16 1.00000000e+00
0.00000000e+00]
[-5.55111512e-17 -1.38777878e-17 0.00000000e+00 5.55111512e-17
1.00000000e+00]]
[[-0.33333333 0.33333333 0.16666667]
[ 0.16666667 1.33333333 -0.33333333]
[ 0.16666667 -0.66666667 0.16666667]]
[[1.00000000e+00 0.00000000e+00 0.00000000e+00]
[2.77555756e-17 1.00000000e+00 2.77555756e-17]
[5.55111512e-17 2.22044605e-16 1.00000000e+00]]
行列と逆行列の積を計算すると、浮動小数点数の計算による誤差がありますが、単位行列になることがわかります。
ところで、実際に連立一次方程式を解くとき、逆行列を使用することはほとんどありません。わざわざ逆行列を求めるよりも、ガウスの消去法を使ったほうがより速く解を求めることができるからです。同じ係数行列の連立方程式を何度も解く場合でも、逆行列を求めるのではなく、次に説明する LU 分解を使うことが多いようです。
●LU分解
連立一次方程式 \(Ax = b\) を解くとき、あらかじめ係数行列 \(A\) を下三角行列 \(L\) (lower) と上三角行列 \(U\) (Upper) に分解しておくと、簡単に解くことができます。\(Ax = b\) は \(LUx = b\) と表すことができるので、\(Ux = y\) とおくと \(Ly = b\) を解くことで \(y\) を求めることができます。次に、求めた \(y\) を使って \(Ux = y\) を解けば、\(Ax = b\) の解 \(x\) を求めることができます。簡単な例を示しましょう。
連立方程式
x + y + z = 10
2x + 4y + 6z = 38
2x + 4z = 14
[[1. 1. 1.] [[1. 0. 0.] [[1. 1. 1.]
[2. 4. 6.] = [2. 1. 0.] @ [0. 2. 4.]
[2. 0. 4.]] [2. -1. 1.]] [0. 0. 6.]]
(1) LU 分解
[[1. 0. 0. 10.] y1 = 10
[2. 1. 0. 38.] y2 = 38 - 2 * 10 = 18
[2. -1. 1. 14.]] y3 = 14 - 2 * 10 + 18 = 12
(2) 前進代入
[[1. 1. 1. 10.] x1 = 10 - x2 - x3 = 3
[0. 2. 4. 18.] x2 = (18 - 4 * x3) / 2 = 5
[0. 0. 6. 12.]] x3 = 12 / 6 = 2
(3) 後退代入
LU 分解はガウスの消去法と同じアルゴリズムで求めることができます。\(U\) は前進消去で得た行列で、\(L\) の要素は前進消去するとき係数を 0 にするための倍率になります。実際に L @ U を計算すると元の行列に戻ります。\(Ly = b\) を解く場合、\(L\) は下三角行列なので前進代入で \(y\) を求めることができます。次に \(Ux = y\) を解く場合、\(U\) は上三角行列なので後退代入で \(x\) を求めることができます。
ガウスの消去法の計算量は \(O(N^3)\) ですが、前進代入と後進代入の計算量は \(O(N^2)\) です。同じ係数行列を何度も使う場合、あらかじめ係数行列を LU 分解しておけば、ガウスの消去法よりも効率的に連立方程式を解くことができるわけです。
実際には \(A\) を \(L\) と \(U\) の二つに分けるのではなく、\(U\) の 0 の部分に \(L\) の要素をセットした行列 \(A'\) を作ることにします。プログラムは次のようになります。
リスト : LU 分解
import numpy as np
def lu(xs):
n = len(xs)
zs = xs.astype(np.float_)
for i in range(n):
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i+1:] -= temp * zs[i, i+1:]
zs[j, i] = temp
return zs
def lu_solver(xs, ys):
n = len(xs)
zs = ys.astype(np.float_)
# 前進代入
for i in range(1, n):
zs[i] -= xs[i, :i] @ zs[:i]
# 後退代入
for i in range(n - 1, -1, -1):
zs[i] -= xs[i, i+1:] @ zs[i+1:]
zs[i] /= xs[i, i]
return zs
# 簡単なテスト
a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
b1 = np.array([100, 272])
b2 = np.array([10, 38, 14])
b3 = np.array([-5, -7, -31, -35])
b4 = np.array([1, 0, 8, 0, 1])
for a, b in [(a1, b1), (a2, b2), (a3, b3), (a4, b4)]:
xs = lu(a)
print(xs)
print(lu_solver(xs, b))
関数 lu() は引数の係数行列 xs を LU 分解します。アルゴリズムはガウスの消去法と同じですが、係数を 0 にするための倍率 temp を zs[j, i] にセットするところが異なります。
関数 lu_solver() は LU 分解した行列 xs と右辺の定数を格納したベクトル ys を受け取り、その解を求めます。前進代入は後退代入とは逆に前から順番に値を求めていきます。行列 L の対角成分 L[i, i] は 1 なので、割り算する必要はありません。後退代入はガウスの消去法のプログラムと同じです。
それでは実行してみましょう。
[[ 1. 1.]
[ 2. 2.]]
[ 64. 36.]
[[ 1. 1. 1.]
[ 2. 2. 4.]
[ 2. -1. 6.]]
[ 3. 5. 2.]
[[ 1. 1. 1. 1.]
[-1. 2. 0. 2.]
[ 8. -2. -6. -3.]
[-8. 6. -1. -6.]]
[ 0. -9. 1. 3.]
[[ 1.00000000e+00 -1.00000000e+00 1.00000000e+00 -1.00000000e+00 1.00000000e+00]
[ 1.20000000e+01 6.00000000e+00 -1.00000000e+01 1.20000000e+01 -1.20000000e+01]
[ 1.00000000e+00 3.33333333e-01 3.33333333e+00 -2.00000000e+00 4.00000000e+00]
[ 1.20000000e+01 3.00000000e+00 6.00000000e+00 -1.20000000e+01 -3.55271368e-15]
[ 4.00000000e+00 1.16666667e+00 2.90000000e+00 2.66666667e-01 -1.60000000e+00]]
[ 0.3125 0. -1.875 3.5 6.0625]
正常に動作していますね。
●ピボット選択を使った LU 分解
LU 分解にピボット選択を適用することも簡単です。次のリストを見てください。
リスト : LU 分解 (ピボット選択バージョン)
import numpy as np
def select_pivot(xs, idx, i):
k = np.abs(xs[i:,i]).argmax() + i
if k != i:
temp = xs[i].copy()
xs[i] = xs[k]
xs[k] = temp
idx[i], idx[k] = idx[k], idx[i]
def lu_pv(xs):
n = len(xs)
zs = xs.astype(np.float_)
idx = list(range(n))
for i in range(n):
select_pivot(zs, idx, i)
for j in range(i + 1, n):
temp = zs[j, i] / zs[i, i]
zs[j, i+1:] -= temp * zs[i, i+1:]
zs[j, i] = temp
return zs, idx
def lu_solver(xs, ys):
n = len(xs)
zs = ys.astype(np.float_)
# 前進代入
for i in range(1, n):
zs[i] -= xs[i, :i] @ zs[:i]
# 後退代入
for i in range(n - 1, -1, -1):
zs[i] -= xs[i, i+1:] @ zs[i+1:]
zs[i] /= xs[i, i]
return zs
def lu_solver_pv(xs, idx, ys):
return lu_solver(xs, ys[idx])
select_pivot() の引数 idx に交換した行の情報をセットします。係数行列の大きさを n とすると、idx の初期値は 0 から n - 1 までの整数を格納したリストになります。i 行と k 行を交換した場合は、idx[i] と idx[k] の値も交換します。関数 lu_pv() はピボット選択しながら LU 分解を行い、L と U を格納した行列 zs とリスト idx を返します。関数 lu_solver_pv() は、lu_pv() が返した行列 xs とインデックスリスト idx を受け取り、引数 ys のベクトルの要素を idx に合わせて並べ替えてから lu_solver() を呼び出すだけです。
それでは実際に試してみましょう。以下に簡単なテストと実行結果を示します。
リスト : 簡単なテスト
a1 = np.array([[1, 1], [2, 4]])
a2 = np.array([[1, 1, 1], [2, 4, 6], [2, 0, 4]])
a3 = np.array([[1, 1, 1, 1], [-1, 1, -1, 1], [8, 4, 2, 1], [-8, 4, -2, 1]])
a4 = np.array([[1, -1, 1, -1, 1], [12, -6, 2, 0, 0], [1, 1, 1, 1, 1], [12, 6, 2, 0, 0], [4, 3, 2, 1, 0]])
a5 = np.array([[0, 2, 4], [1, 1, 1], [4, 2, 6]])
b1 = np.array([100, 272])
b2 = np.array([10, 38, 14])
b3 = np.array([-5, -7, -31, -35])
b4 = np.array([1, 0, 8, 0, 1])
b5 = np.array([14, 10, 38])
for a, b in [(a1, b1), (a2, b2), (a3, b3), (a4, b4), (a5, b5)]:
xs, idx = lu_pv(a)
print(xs)
print(idx)
print(lu_solver_pv(xs, idx, b))
[[ 2. 4. ]
[ 0.5 -1. ]]
[1, 0]
[ 64. 36.]
[[ 2. 4. 6. ]
[ 1. -4. -2. ]
[ 0.5 0.25 -1.5 ]]
[1, 2, 0]
[ 3. 5. 2.]
[[ 8. 4. 2. 1. ]
[-1. 8. 0. 2. ]
[ 0.125 0.0625 0.75 0.75 ]
[-0.125 0.1875 -1. 1.5 ]]
[2, 3, 0, 1]
[ 0. -9. 1. 3.]
[[ 12. -6. 2. 0. 0. ]
[ 1. 12. 0. 0. 0. ]
[ 0.33333333 0.41666667 1.33333333 1. 0. ]
[ 0.08333333 -0.04166667 0.625 -1.625 1. ]
[ 0.08333333 0.125 0.625 -0.23076923 1.23076923]]
[1, 3, 4, 0, 2]
[ 0.3125 0. -1.875 3.5 6.0625]
[[ 4. 2. 6. ]
[ 0. 2. 4. ]
[ 0.25 0.25 -1.5 ]]
[2, 0, 1]
[ 5. 3. 2.]
正常に動作していますね。
●参考文献・URL
- 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
- ガウスの消去法 - Wikipedia
