連立方程式を解くときには,ガウスの掃き出し法(Gaussian elimination)が利用される. これはalgorithmとしてはlower, upper triangle matrixに 分ける操作からLU分解(LU decomposition)と呼ばれる.
import numpy as np
import scipy.linalg
np.set_printoptions(precision=3, suppress=True)
a=np.array([[1,1],[2,4]])
print(a)
P, L, U = scipy.linalg.lu(a)
print(P)
print(L)
print(U)
[[1 1] [2 4]] [[ 0. 1.] [ 1. 0.]] [[ 1. 0. ] [ 0.5 1. ]] [[ 2. 4.] [ 0. -1.]]
a=np.array([[1,1],[2,4]])
b = np.array([5,18])
ab = np.column_stack((a,b))
P, L, U = scipy.linalg.lu(ab)
print(P)
print(L)
print(U)
[[ 0. 1.] [ 1. 0.]] [[ 1. 0. ] [ 0.5 1. ]] [[ 2. 4. 18.] [ 0. -1. -4.]]
例えば,次のような連立方程式を解くと $$ x - y -z = 1 \\ x - y +z = -1 \\ x + y -z = -1 $$
a=np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
b = np.array([1,-1,-1])
ab = np.column_stack((a,b))
P, L, U = scipy.linalg.lu(ab)
#print(P)
#print(L)
print(U)
[[ 1. -1. -1. 1.] [ 0. 2. 0. -2.] [ 0. 0. 2. -2.]]
a=np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
b = np.array([1,-1,-1])
scipy.linalg.inv(a).dot(b)
array([-1., -1., -1.])
となり, $$ x=-1, y=-1, z=-1 $$ が得られる.これの意味はつぎのように理解できる.
まず,先ほどの連立方程式を$z$について形式的に解く.すると, $$ z = x - y - 1 \\ z = -1-x+y \\ z = 1+x+y $$ となる.これをplot3dしてみると次の通り,表示される.
%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
def f(x,y):
return x - y - 1
def g(x,y):
return -1-x+y
def h(x,y):
return 1+x+y
x = np.arange(-2, 2, 0.25)
y = np.arange(-2, 2, 0.25)
X, Y = np.meshgrid(x, y)
Z1 = f(X,Y)
Z2 = g(X,Y)
Z3 = h(X,Y)
fig = plt.figure()
plot3d = Axes3D(fig)
plot3d.plot_surface(X,Y,Z1, alpha = 1)
plot3d.plot_surface(X,Y,Z2, color='r', alpha = 1)
plot3d.plot_surface(X,Y,Z3, color='y', alpha = 1)
plt.show()
# x=0.439423 , y=2.21846 , z=-7.77475
pythonのplot3dではrenderingがちゃんとできてないので,見にくいが, ちゃんとrenderingしてくれる画像でみると,
であることがわかると思う.
3次元のちゃんとしたrenderingにはmayaviを使えというお達しだが,
pip install mayavi
でやってもvtkが入ってないとかで怒られる.これは諦める. そうなると,3dplotが貧弱だが,まあ,仕方ないか.
Mapleの優位さを宣伝する部分だね.
sagaを使えということか...そうなるとまた,教研でinstallしなあかんのが増えるな.