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.]])
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]]
array([[ 6.], [-1.], [ 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.]])
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]])
np.dot(L,U)
array([[ 2., -2., -1., 11.], [ 1., 4., 3., 11.], [ 1., -2., 1., 11.]])
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.]])
np.dot(L,A)
array([[ 1., 4., 3.], [ 1., -2., 1.], [ 2., -2., -1.]])
pprint(L)
array([[ 1. , 0. , 0. ], [ 1. , 1. , 0. ], [ 2. , 1.66666667, 1. ]])
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.]])
[-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]
# 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]
# 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]