#!/usr/bin/env python
# coding: utf-8
#
# 線形代数(Linear Algebra I)
#
#
# cc by Shigeto R. Nishitani, 2018-03-15
#
#
# * file: /Users/bob/python/doing_math_with_python/linear_algebra/LA-I.ipynb
# # Table of Contents
#
# # LU 分解
#
# 連立方程式を解くときには,**ガウスの掃き出し法(Gaussian elimination)**が利用される.
# これはalgorithmとしてはlower, upper triangle matrixに
# 分ける操作から**LU分解(LU decomposition)**と呼ばれる.
#
# * **階段行列(echelon form)**は左側に連続して0が並んでいて,行とともにその数が増えていく行列
# * r+1以下の行の成分が全て0であるとき,**階数(rank)**がrであるという.
# * **係数行列**$A$,**拡大係数行列**$[Ab]$
# * 連立一次方程式,
# * 解の存在定理$rank[Ab]=rank A$
# * $n$次正方行列が**正則(regular)**,**次数**と階数が一致する行列$rank A = n$.連立方程式がただ1組の解を持つ条件.
# * **同次連立一次方程式** $Ax = 0$
# * **自明解**$x=0$,**非自明解**$x \neq 0$,
# * 連立方程式$Ax=b$の**一般解**$x$について,$x=x_0+\alpha x_1+\beta x_2$としたときの$x_0$が**特殊解**,$x_1, x_2$は**同伴な同次連立方程式**$Ax=0$の**基本解**.$x$はそれらの**線型結合(linear combination)**で表されている.
# # 連立方程式,掃き出しとplot
#
# ![LA-I](./figs/LinearAlgebra-I.png)
# ## 拡大係数行列の作成と掃き出し(LU分解)
# In[1]:
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)
# In[2]:
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)
# ## 3次元の連立方程式の意味とplot
# 例えば,次のような連立方程式を解くと
# $$
# x - y -z = 1 \\
# x - y +z = -1 \\
# x + y -z = -1
# $$
# In[3]:
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)
# In[4]:
a=np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
b = np.array([1,-1,-1])
scipy.linalg.inv(a).dot(b)
# となり,
# $$
# x=-1, y=-1, z=-1
# $$
# が得られる.これの意味はつぎのように理解できる.
#
# まず,先ほどの連立方程式を$z$について形式的に解く.すると,
# $$
# z = x - y - 1 \\
# z = -1-x+y \\
# z = 1+x+y
# $$
# となる.これをplot3dしてみると次の通り,表示される.
# In[2]:
get_ipython().run_line_magic('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してくれる画像でみると,
# 1. それぞれの方程式が表す面が表示され,
# 1. 面同士がintersectするところが直線をなし.
# 1. それらが交わった点が$x=-1, y=-1, z=-1$
# であることがわかると思う.
# 3次元のちゃんとしたrenderingにはmayaviを使えというお達しだが,
# > pip install mayavi
#
# でやってもvtkが入ってないとかで怒られる.これは諦める.
# そうなると,3dplotが貧弱だが,まあ,仕方ないか.
#
# Mapleの優位さを宣伝する部分だね.
#
# ![three_planes_plot](./figs/three_planes_plot.png)
#
# sagaを使えということか...そうなるとまた,教研でinstallしなあかんのが増えるな.
# In[ ]: