3x3行列のLU分解と解

In [1]:
from pprint import pprint
import scipy.linalg as linalg   # SciPy Linear Algebra Library
import numpy as np

# A = np.array([[1,1,-2,-4],[1,-2,1,5],[2,-2,-1,2]])
A = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])
b = np.array([[11.0],[11],[11]])
pprint(A)
pprint(b)
array([[ 1.,  4.,  3.],
       [ 1., -2.,  1.],
       [ 2., -2., -1.]])
array([[ 11.],
       [ 11.],
       [ 11.]])
In [2]:
inv_A = linalg.inv(A)
print(inv_A)
np.dot(inv_A,b)
[[ 0.18181818 -0.09090909  0.45454545]
 [ 0.13636364 -0.31818182  0.09090909]
 [ 0.09090909  0.45454545 -0.27272727]]
Out[2]:
array([[ 6.],
       [-1.],
       [ 3.]])
In [3]:
e_A = np.hstack((A, b))
pprint(e_A)
pprint(np.column_stack((A,b)))
array([[  1.,   4.,   3.,  11.],
       [  1.,  -2.,   1.,  11.],
       [  2.,  -2.,  -1.,  11.]])
array([[  1.,   4.,   3.,  11.],
       [  1.,  -2.,   1.,  11.],
       [  2.,  -2.,  -1.,  11.]])
In [4]:
P, L, U = linalg.lu(e_A)

print("P:")
pprint(P)

print("L:")
pprint(L)

print("U:")
pprint(U)
P:
array([[ 0.,  1.,  0.],
       [ 0.,  0.,  1.],
       [ 1.,  0.,  0.]])
L:
array([[ 1. ,  0. ,  0. ],
       [ 0.5,  1. ,  0. ],
       [ 0.5, -0.2,  1. ]])
U:
array([[  2. ,  -2. ,  -1. ,  11. ],
       [  0. ,   5. ,   3.5,   5.5],
       [  0. ,   0. ,   2.2,   6.6]])
In [5]:
np.dot(L,U)
Out[5]:
array([[  2.,  -2.,  -1.,  11.],
       [  1.,   4.,   3.,  11.],
       [  1.,  -2.,   1.,  11.]])
In [1]:
from pprint import pprint
import scipy.linalg as linalg   # SciPy Linear Algebra Library
import numpy as np


A = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])
A0 = np.array([[1.0,4,3],[1,-2,1],[2,-2,-1]])
b = np.array([[11.0],[11],[11]])


n = 3
L = np.identity(n)
T = []
for i in range(n):        #i行目
    T.append(np.identity(n))
    for j in range(i+1, n):
        am = A[j,i]/A[i,i]    #i行の要素を使って,i+1行目の先頭を消す係数を求める
        T[i][j,i]=-am         #i番目の消去行列に要素を入れる
        L[j,i]=am             #LTMの要素
        for k in range(n):
            A[j,k] -= am*A[i,k] #もとの行列をUTMにしていく
        b[j] -= b[i]*am   #数値ベクトルも操作
pprint(A)
pprint(b)
array([[ 1.        ,  4.        ,  3.        ],
       [ 0.        , -6.        , -2.        ],
       [ 0.        ,  0.        , -3.66666667]])
array([[ 11.],
       [  0.],
       [-11.]])
In [2]:
np.dot(L,A)
Out[2]:
array([[ 1.,  4.,  3.],
       [ 1., -2.,  1.],
       [ 2., -2., -1.]])
In [3]:
pprint(L)
array([[ 1.        ,  0.        ,  0.        ],
       [ 1.        ,  1.        ,  0.        ],
       [ 2.        ,  1.66666667,  1.        ]])

pivot付きのLU分解

In [5]:
from pprint import pprint
import scipy.linalg as linalg   # SciPy Linear Algebra Library
import numpy as np

A=np.array([[3,2,2,1],[3,2,3,1],[1,-2,-3,1],[5,3,-2,5]])
b=np.array([-6,2,-9,2])
ex_A = np.column_stack((A,b))
P,L,U = linalg.lu(ex_A)
pprint(P)
pprint(L)
pprint(U)

pprint(np.dot(P,A))
pprint(np.dot(L,U))
array([[ 0.,  0.,  0.,  1.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 1.,  0.,  0.,  0.]])
array([[ 1.      ,  0.      ,  0.      ,  0.      ],
       [ 0.2     ,  1.      ,  0.      ,  0.      ],
       [ 0.6     , -0.076923,  1.      ,  0.      ],
       [ 0.6     , -0.076923,  0.75    ,  1.      ]])
array([[ 5.      ,  3.      , -2.      ,  5.      ,  2.      ],
       [ 0.      , -2.6     , -2.6     , -0.      , -9.4     ],
       [ 0.      ,  0.      ,  4.      , -2.      ,  0.076923],
       [ 0.      ,  0.      ,  0.      , -0.5     , -7.980769]])
array([[ 5.,  3., -2.,  5.],
       [ 1., -2., -3.,  1.],
       [ 3.,  2.,  3.,  1.],
       [ 3.,  2.,  2.,  1.]])
array([[ 5.,  3., -2.,  5.,  2.],
       [ 1., -2., -3.,  1., -9.],
       [ 3.,  2.,  3.,  1.,  2.],
       [ 3.,  2.,  2.,  1., -6.]])

Jacobi vs Gauss-Seidel

[-1.2       0.666667  0.777778  0.6     ]
[-1.608889  0.607407  0.562963  0.751111]
[-1.584296  0.764938  0.54749   0.782519]
[-1.618989  0.751429  0.518705  0.676892]
[-1.589405  0.807797  0.506116  0.680422]
In [6]:
# Jacobi iterative method for solving Ax=b by bob
import numpy as np
np.set_printoptions(precision=6, suppress=True)

A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b=np.array([-6,2,-7,3])
n=4
x0=np.zeros(n)
x1=np.zeros(n)

for iter in range(0, 5):
    for i in range(0, n):
        x1[i]=b[i]
        for j in range(0, n):
            x1[i]=x1[i]-A[i][j]*x0[j]
        x1[i]=x1[i]+A[i][i]*x0[i]
        x1[i]=x1[i]/A[i][i]
        x0[i]=x1[i]
    print(x0)
[-1.2       1.066667  0.407407  0.362963]
[-1.567407  0.932346  0.436763  0.528779]
[-1.579578  0.871345  0.46739   0.580064]
[-1.58376   0.845435  0.478382  0.600844]
[-1.584932  0.835236  0.482827  0.608976]
In [15]:
# Jacobi iterative method for solving Ax=b by bob
import numpy as np
np.set_printoptions(precision=6, suppress=True)

A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b=np.array([-6,2,-7,3])
n=4
x0=np.zeros(n)

for iter in range(0, 5):
    for i in range(0, n):
        x1 = b[i]
        for j in range(0, n):
            x1 -= A[i][j]*x0[j]
        x1 += A[i][i]*x0[i]
        x1 /= A[i][i]
        x0[i] = x1
    print(x0)
[-1.2       1.066667  0.407407  0.362963]
[-1.567407  0.932346  0.436763  0.528779]
[-1.579578  0.871345  0.46739   0.580064]
[-1.58376   0.845435  0.478382  0.600844]
[-1.584932  0.835236  0.482827  0.608976]
In [ ]: