第10回：行列・線形代数

▶ ベクトルの線形結合からなる格子点

複数のベクトルの線形結合とは， それらのベクトルのスカラー倍を加え合わせたものを，それらのベクトルの線形結合という．

二つのベクトル $a_1=\begin{bmatrix} 1 \\\\ 0 \end{bmatrix}, a_2=\begin{bmatrix} 0 \\\\ 1 \end{bmatrix}$ の「整数」係数の線形結合による格子点を描く．

さらに， $b_1=\begin{bmatrix} \dfrac{1}{2} \\\\ \dfrac{1}{2} \end{bmatrix}, b_2=\begin{bmatrix} \dfrac{1}{2} \\\\ -\dfrac{1}{2} \end{bmatrix}$ の「整数」係数の線形結合からなる格子点を重ねる．

using PyPlot
plt.axes().set_aspect("equal")
a1 = [1, 0];
a2 = [0, 1];
for m = -3:3, n = -3:3
   r = m * a1 + n * a2
   plt.plot(r[1], r[2], "bo")
end

b1 = [1 / 2, 1 / 2];
b2 = [1 / 2, -1 / 2];
for m = -3:3, n = -3:3
   r = m * b1 + n * b2
   plt.plot(r[1], r[2], "r.")
end

plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

どちらも正方格子（cubic lattice）であるが，座標系の取り方が異なる． 基底 $b_1, b_2$ で張られる格子点は， 基底 $a_1, a_2$ で張られる格子点の中央の点も含んでいることが観察できる．

今度は， $c_1=\begin{bmatrix} 1 \\ 0 \end{bmatrix}, c_2=\begin{bmatrix} -\dfrac{1}{2} \\ \dfrac{\sqrt{3}}{2} \end{bmatrix}$ で張られる格子点を描いてみる． これは，六方格子（hexagonal lattice）と呼ばれる．

using PyPlot
plt.axes().set_aspect("equal")

c1 = [1, 0];
c2 = [-1 / 2, sqrt(3) / 2];
for m = -3:3, n = -3:3
   r = m * c1 + n * c2
   plt.plot(r[1], r[2], "g.")
end

plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

■ 平面ベクトルの内積

関数 dot(a,b) は，ベクトル a と b との「内積（inner product）」を返す．

中置き演算子 ⋅ を用いて a⋅b と書くこともできる． 「 ⋅ 」 は，バックスラッシュ \ に続けて cdots と入力してから，TABキーを押すことによって入力できる． かな漢字変換システムで入力できる「・」（中黒=なかぐろ）とは，別の文字である．

▶ 平面ベクトルどうしのなす角を求める

ベクトル $a$ と $b$の内積は，$a$ と $b$ のなす角$\theta$ を用いて，以下のように定義される．

\[a\cdot b = \left\vert{a}\right\vert \left\vert{b}\right\vert \cos\theta\]

これから，$\theta$ を求めるには，次の式を用いればよい．

\[\cos\theta = \dfrac{a\cdot b}{ \left\vert{a}\right\vert \left\vert{b}\right\vert }\]

また，内積の定義から，自分自身の内積は，ノルムの二乗であることも分かる．

\[a\cdot a = \left\vert{a}\right\vert^2\]

▶ 例：ベクトルどうしのなす角度を求める

上で出てきたベクトルのうち，a1, a2, c1, c2 のノルムは $1$ である．

julia> using LinearAlgebra
julia> #
       a1 = [1, 0];
julia> a2 = [0, 1];
julia> b1 = [1 / 2, 1 / 2];
julia> b2 = [1 / 2, -1 / 2];
julia> c1 = [1, 0];
julia> c2 = [-1 / 2, sqrt(3) / 2];
julia> norm(a1)1.0
julia> norm(a2)1.0
julia> norm(c1)1.0
julia> norm(c2)0.9999999999999999

b1 と b2 のノルムは $\dfrac{1}{\sqrt{2}}$ である． 自分自身の内積の値と比較しよう

julia> b1 ⋅ b10.5
julia> norm(b1)0.7071067811865476
julia> norm(b1)^20.5000000000000001
julia> b1 ⋅ b10.5
julia> norm(b2)0.7071067811865476
julia> norm(b2)^20.5000000000000001

内積から算出した $\cos\theta$から角度 $\theta$ を得るには，関数 acos() を用いる． 関数 acos(r) は $\cos \theta = r$ となる $\theta$ をラジアンで返す． 関数 acosd(r) は，$\theta$ を角度で返す．

これらのベクトルのなす角度を算出しよう． a1 と a2，および，b1 と b2 は直交している． a1 と b1 は，45° をなしている． c1 と c2 は，120° をなしている，ことが計算できた．

julia> acosd(a1 ⋅ a2)90.0
julia> acosd(b1 ⋅ b2 / norm(b1) / norm(b2))90.0
julia> acosd(a1 ⋅ b1 / norm(a1) / norm(b1))45.00000000000001
julia> acosd(c1 ⋅ c2)120.00000000000001

■ タプル

タプル（tuple）は，複数の値をカンマ , で区切って並べ，括弧 () で囲んだものである． ベクトルと似たように使えるが，要素を更新することはできない．

julia> # 要素 1つのタプル
       (1,)(1,)
julia> # 要素 2つのタプル
       (1, 2)(1, 2)
julia> # 要素 3つのタプル
       a = (1, 2, 3)(1, 2, 3)
julia> # タプルの長さ
       length(a)3
julia> # タプルの要素
       a[2]2
julia> # 更新はできない
       a[2] = 3 # => MethodErrorERROR: MethodError: no method matching setindex!(::Tuple{Int64, Int64, Int64}, ::Int64, ::Int64)

関数には，複数の値を返すものがある．このとき，タプルが用いられる．

例えば，divrem(x,d) は，div(x,d) と rem(x,d) の二つの値を返す．

julia> divrem(5, 3)(1, 2)

タプルを右辺において，複数の変数に同時に代入できる．

julia> x, y = (1, 2, 3)(1, 2, 3)
julia> x1
julia> y2

■ 行列

要素を ; で区切って列挙したものを， 大かっこ [] で囲むと，行列を作ることができる．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

3行 4列の行列

■ ベクトルも，■ 行列も，「配列（array）」として表されている． ベクトルと同じ関数が用いられる．

julia> # 寸法 => タプル
       size(a)(3, 4)
julia> # 第1軸 = 列の寸法
       size(a, 1)3
julia> # 第2軸 = 行の寸法
       size(a, 2)4
julia> # 全要素数
       length(a)12

■ 行列の転置

transpose(a) は，行列 a を転置する．すなわち，行と列を入れ替える．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34
julia> size(a)(3, 4)
julia> b = transpose(a)4×3 transpose(::Matrix{Int64}) with eltype Int64:
 11  21  31
 12  22  32
 13  23  33
 14  24  34
julia> size(b)(4, 3)

▶ 行列のスカラー倍・スカラーの和差

以下，しばらく，2x2 の正方行列を例に説明する．

julia> a = [11 12; 21 22]2×2 Matrix{Int64}:
 11  12
 21  22
julia> a * 22×2 Matrix{Int64}:
 22  24
 42  44
julia> a .+ 22×2 Matrix{Int64}:
 13  14
 23  24
julia> a .- 22×2 Matrix{Int64}:
  9  10
 19  20

▶ 行列に列ベクトルを加減

julia> a = [11 12; 21 22]2×2 Matrix{Int64}:
 11  12
 21  22
julia> v = [1, -1]2-element Vector{Int64}:
  1
 -1
julia> a .+ v2×2 Matrix{Int64}:
 12  13
 20  21
julia> a .- v2×2 Matrix{Int64}:
 10  11
 22  23

▶ 行列どうしの加減

julia> a = [11 12; 21 22]2×2 Matrix{Int64}:
 11  12
 21  22
julia> b = a * 22×2 Matrix{Int64}:
 22  24
 42  44
julia> a + b2×2 Matrix{Int64}:
 33  36
 63  66
julia> a - b2×2 Matrix{Int64}:
 -11  -12
 -21  -22

■ 添字を用いた行列の要素の読み書き

行列の添字は， 第1軸（列）と第2軸（行）の番号を，カンマ , で区切って並べ，大かっこ [] で囲んだものである．

ベクトルと同じように，添字で示された要素の読み出し， 添字で示された要素の書き換えができる．

julia> # 添字による要素の読み出し
       a[2, 2]22
julia> # 行列の要素の更新
       a[1, 2] = 3030
julia> a2×2 Matrix{Int64}:
 11  30
 21  22

■ 部分行列

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34
julia> # 列を取り出す
       a[:, 2]3-element Vector{Int64}:
 12
 22
 32
julia> # 行を取り出す
       a[2, :]4-element Vector{Int64}:
 21
 22
 23
 24
julia> # 部分行列
       a[1:2, 1:2]2×2 Matrix{Int64}:
 11  12
 21  22
julia> a[2:3, 2:3]2×2 Matrix{Int64}:
 22  23
 32  33

▶ 行列に入れた点座標で図形を描画する

using PyPlot
plt.axes().set_aspect("equal")

xy = [1 2 2 1; 1 1 3 1]
@show xy
plt.plot(xy[1, :], xy[2, :])
xy = xy .+ [1 / 2, 1 / 2]
plt.plot(xy[1, :], xy[2, :])
plt.xlim(0, 4)
plt.ylim(0, 4)
xy = [1 2 2 1; 1 1 3 1]

■ 行列とベクトルの積

julia> a = [11 12; 21 22]2×2 Matrix{Int64}:
 11  12
 21  22
julia> v = [1, -1]2-element Vector{Int64}:
  1
 -1
julia> a * v2-element Vector{Int64}:
 -1
 -1

▶ 回転行列とベクトルの積

以下の形の行列を，平面空間に対する回転行列という．

\[R(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos \theta \end{bmatrix}\]

回転行列とベクトルの積は， そのベクトルを，原点の周りに 反時計方向に角 $\theta$ だけ回転する写像に対応する．

\[x^{\prime} = R(\theta) x\]

julia> # 回転行列
       r15 = [cosd(15) -sind(15); sind(15) cosd(15)]2×2 Matrix{Float64}:
 0.965926  -0.258819
 0.258819   0.965926
julia> xy = [1, 0]2-element Vector{Int64}:
 1
 0
julia> xy = r15 * xy2-element Vector{Float64}:
 0.9659258262890683
 0.25881904510252074
julia> xy = r15 * xy2-element Vector{Float64}:
 0.8660254037844387
 0.49999999999999994

これらを描こう．軌跡は円を描いた．

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1, 0]

for i = 0:20
   global xy
   plt.plot(xy[1], xy[2], "o")
   xy = r15 * xy
end

plt.xlim(-1.2, 1.2)
plt.ylim(-1.2, 1.2)
plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

原点以外の点 $c$ の周りで回転する場合は，回転の中心をずらして，

\[\begin{aligned} (x^{\prime}-c) & = R(\theta) (x-c) \\ x^{\prime} & = c + R(\theta) (x-c) \end{aligned}\]

とすればよい．

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1, 0]
c = [1 / 2, 1 / 2]

for i = 0:20
   global xy
   plt.plot(xy[1], xy[2], "o")
   xy = c + r15 * (xy - c)
end

plt.axvline(c[1], color = "k", lw = 0.5)
plt.axhline(c[2], color = "k", lw = 0.5)

plt.xlim(-1, 2)
plt.ylim(-1, 2)

■ 行列と行列の積

julia> a = [1 2; 3 4]2×2 Matrix{Int64}:
 1  2
 3  4
julia> b = [5 6; 7 8]2×2 Matrix{Int64}:
 5  6
 7  8
julia> a * b2×2 Matrix{Int64}:
 19  22
 43  50

▶ 座標を行列に格納した図形を回転する

using PyPlot
plt.axes().set_aspect("equal")
r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1 2 2 1; 1 1 3 1]

for i = 0:20
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = r15 * xy
end

plt.xlim(-4, 4)
plt.ylim(-4, 4)
plt.axhline(0, color = "k", lw = 0.5)
plt.axvline(0, color = "k", lw = 0.5)

回転中心をずらしてみる

plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
xy = [1 2 2 1; 1 1 3 1]
c = [1 / 2, 1 / 2]

for i = 0:20
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = c .+ r15 * (xy .- c)
end

plt.axvline(c[1], color = "k", lw = 0.5)
plt.axhline(c[2], color = "k", lw = 0.5)

plt.xlim(-4, 4)
plt.ylim(-4, 4)

■ いろいろな行列の生成

■ 要素が 0 の行列を作る

関数 zeros は，要素が零 $0$ の行列を作るのに使える．

  • 関数 zeros(n,m) は，要素の型が浮動小数点で，寸法 (n,m) の行列を作る．
  • 関数 zeros(T, n,m) は，要素の型が T で，寸法 (n,m) の行列を作る．
julia> zeros(3, 4) # 要素は浮動小数点3×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
julia> zeros(Float64, 3, 4) # 上と同じ3×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
julia> zeros(Int64, 3, 4) # 要素は整数3×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0

行列 a と同じ寸法を持つ 0 ベクトルを作るには，

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34
julia> zeros(Int64, size(a))3×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0
julia> zeros(Float64, size(a))3×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

関数 zero.() を以下のように用いれば 行列 a の要素の型と同じ要素の型を持ち，a と寸法が等しい 0 ベクトルを作れる．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34
julia> zero.(a)3×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0

■ 要素が 1 の行列を作る

関数 ones は，要素が零 $1$ の行列を作るのに使える．

  • 関数 ones(n,m) は，要素の型が浮動小数点で，寸法 (n,m) の行列を作る．
  • 関数 ones(T, n,m) は，要素の型が T で，寸法 (n,m) の行列を作る．
julia> ones(3, 4) # 要素は浮動小数点3×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
julia> ones(Float64, 3, 4) # 上と同じ3×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
julia> ones(Int64, 3, 4) # 要素は整数3×4 Matrix{Int64}:
 1  1  1  1
 1  1  1  1
 1  1  1  1

行列 a と同じ寸法を持つ 1 行列を作るには，

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34
julia> ones(Int64, size(a))3×4 Matrix{Int64}:
 1  1  1  1
 1  1  1  1
 1  1  1  1
julia> ones(Float64, size(a))3×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

関数 one.() を以下のように用いれば 行列 a の要素の型と同じ要素の型を持ち，a と寸法が等しい 1 行列を作れる．

julia> a = [11 12 13 14; 21 22 23 24; 31 32 33 34]3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34
julia> one.(a)3×4 Matrix{Int64}:
 1  1  1  1
 1  1  1  1
 1  1  1  1

■ 対角要素を指定して，正方行列をつくる

julia> using LinearAlgebra
julia> Diagonal([1, 2, 3])3×3 LinearAlgebra.Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

■ 疑似乱数を要素とする行列を作る

julia> rand(3, 3)3×3 Matrix{Float64}:
 0.488932  0.790778  0.932983
 0.950551  0.815109  0.746246
 0.415108  0.60719   0.451398

■ ベクトルをまとめて行列を作る

julia> a = [1, 2]2-element Vector{Int64}:
 1
 2
julia> b = [3, 4]2-element Vector{Int64}:
 3
 4
julia> c = [5, 6]2-element Vector{Int64}:
 5
 6
julia> [a b c]2×3 Matrix{Int64}:
 1  3  5
 2  4  6
julia> m = hcat(a, b, c)2×3 Matrix{Int64}:
 1  3  5
 2  4  6
julia> size(m)(2, 3)

さらに転置をとろう．

julia> mt = transpose(m)3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  2
 3  4
 5  6
julia> size(mt)(3, 2)

● 内包表現を用いて行列を作る

julia> [[t, t .^ 2] for t in [0, 2, 4]]3-element Vector{Vector{Int64}}:
 [0, 0]
 [2, 4]
 [4, 16]
julia> hcat([[t, t .^ 2] for t in [0, 2, 4]]...)2×3 Matrix{Int64}:
 0  2   4
 0  4  16

▶ 楕円を描く・回転する

using PyPlot
plt.axes().set_aspect("equal")

ts = 0:pi/18:2pi
xs = 2 * cos.(ts)
ys = sin.(ts)
xy = transpose(hcat(xs, ys))
plt.plot(xy[1, :], xy[2, :])

plt.xlim(-2.2, 2.2)
plt.ylim(-2.2, 2.2)

Note

楕円の座標を含む行列 xy を作るのに，以下のように● 内包表記を使ってもよい．

julia> xy=hcat([ [2*cos.(t), sin(t)] for t=0:pi/18:2pi]...)2×37 Matrix{Float64}:
 2.0  1.96962   1.87939  1.73205  …   1.87939   1.96962    2.0
 0.0  0.173648  0.34202  0.5         -0.34202  -0.173648  -2.44929e-16

回転する

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
ts = 0:pi/18:2pi
xy = transpose(hcat(2 * cos.(ts), sin.(ts)))

for i = 0:5
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = r15 * xy
end

■ 行列式と逆行列

julia> a = [1 2; 2 3]2×2 Matrix{Int64}:
 1  2
 2  3
julia> det(a) # 行列式-1.0
julia> inv(a) # 逆行列2×2 Matrix{Float64}:
 -3.0   2.0
  2.0  -1.0
julia> a^(-1) # 逆行列2×2 Matrix{Float64}:
 -3.0   2.0
  2.0  -1.0

■ 行列の商（2x2 行列）

行列 $A$とベクトル $x$ に対して，$x=Ay$ を満たすようなベクトル $y$を，$x$を$A$で割った商という．商を求める演算が用意されている．この演算は，$A$の逆行列を計算するよりも，計算コストの上で有利である．

julia> a = [1 2; 2 3]2×2 Matrix{Int64}:
 1  2
 2  3
julia> v = [1, 1]2-element Vector{Int64}:
 1
 1
julia> a \ v2-element Vector{Float64}:
 -1.0
  1.0
julia> w = [3, 5]2-element Vector{Int64}:
 3
 5
julia> a \ w2-element Vector{Float64}:
 1.0
 1.0

行列 $A$ が正則でない（逆行列が存在しない）場合には，例外が発生する．

julia> c = [1 2; 2 4]2×2 Matrix{Int64}:
 1  2
 2  4
julia> det(c)0.0
julia> c \ v # => SingularExceptionERROR: LinearAlgebra.SingularException(2)

同様に，行列 $X$ に対して，$X=AY$ を満たすような行列 $Y$ を，$X$ を $A$ で割った商という．

julia> b = [1 3; 1 5]2×2 Matrix{Int64}:
 1  3
 1  5
julia> a \ b2×2 Matrix{Float64}:
 -1.0  1.0
  1.0  1.0

▶ 楕円を，逆に回転する

行列の商を用いて，楕円を逆回転してみよう．

using PyPlot
plt.axes().set_aspect("equal")

r15 = [cosd(15) -sind(15); sind(15) cosd(15)]
ts = 0:pi/18:2pi
xy = transpose(hcat(2 * cos.(ts), sin.(ts)))

for i = 0:5
   global xy
   plt.plot(xy[1, :], xy[2, :])
   xy = r15 \ xy
end

行列式が 0 の行列は，正則ではない

julia> a = [1 2; 2 4]2×2 Matrix{Int64}:
 1  2
 2  4
julia> det(a)0.0
julia> v = [1, 1]2-element Vector{Int64}:
 1
 1
julia> # 例外を発生する
       a \ vERROR: LinearAlgebra.SingularException(2)
julia> # 例外を発生する
       inv(a)ERROR: LinearAlgebra.SingularException(2)

▶ 空間ベクトル：なす角を求める

空間ベクトルとは，寸法 3 のベクトルである． 内積が $0$ なら，それらのベクトルは直交である．

例: 以下の３つのベクトルが，互いに直行することを示せ．

julia> using LinearAlgebra
julia> #
       a = [1 / 2, 1 / 2 + sqrt(2) / 4, 1 / 2 - sqrt(2) / 4]3-element Vector{Float64}:
 0.5
 0.8535533905932737
 0.1464466094067262
julia> b = [-1 / 2, 1 / 2 - sqrt(2) / 4, 1 / 2 + sqrt(2) / 4]3-element Vector{Float64}:
 -0.5
  0.1464466094067262
  0.8535533905932737
julia> c = [1 / sqrt(2), -1 / 2, 1 / 2]3-element Vector{Float64}:
  0.7071067811865475
 -0.5
  0.5
julia> norm(a)1.0
julia> norm(b)1.0
julia> norm(c)0.9999999999999999
julia> dot(a, b)-5.2974517141865445e-17
julia> a ⋅ b-5.2974517141865445e-17
julia> b ⋅ c0.0
julia> c ⋅ a-2.7755575615628914e-17

二つのベクトルのなす角を求めよ．

julia> a = [-3, 1, 2]3-element Vector{Int64}:
 -3
  1
  2
julia> b = [2, -3, 1]3-element Vector{Int64}:
  2
 -3
  1
julia> ab = a ⋅ b-7
julia> na = norm(a)3.7416573867739413
julia> nb = norm(b)3.7416573867739413
julia> r = ab / na / nb-0.5
julia> # ラジアン単位
       acos(r)2.0943951023931957
julia> # 角度単位
       acosd(r)120.00000000000001
julia> # 一行で書ける
       acosd(a ⋅ b / norm(a) / norm(b))120.00000000000001

■ 空間ベクトルの外積

関数 cross(a,b) は，空間ベクトル a と b との外積または「ベクトル積（outer product）」を返す．

中置き演算子 × を用いて a×b と書く こともできる． × は，バックスラッシュ \ に times と入力してから， TABキーを押すことによって入力できる． かな漢字変換システムで入力できる「✕」 とは，別の文字である．

外積 $a\times b$ の向きは，$a$ から $b$ へ回転したとき，右ねじが進む方向である．

外積 $a\times b$ の大きさは，$a$ と $b$ のなす角 $\theta$ を用いて，以下のように定義される． これは，ベクトル a と b がなす平行四辺形の面積である．

\[\left\vert a\times b\right\vert = \left\vert{a}\right\vert \left\vert{b}\right\vert \sin\theta\]

▶ 空間座標の基本単位ベクトル

julia> a = [1, 0, 0]3-element Vector{Int64}:
 1
 0
 0
julia> b = [0, 1, 0]3-element Vector{Int64}:
 0
 1
 0
julia> c = [0, 0, 1]3-element Vector{Int64}:
 0
 0
 1
julia> cross(a, b)3-element Vector{Int64}:
 0
 0
 1
julia> # a×b = c
       a × b3-element Vector{Int64}:
 0
 0
 1
julia> # b×c = a
       b × c3-element Vector{Int64}:
 1
 0
 0
julia> # c×a = b
       c × a3-element Vector{Int64}:
 0
 1
 0

別の正規直交系の例

julia> a = [1 / 2, 1 / 2 + sqrt(2) / 4, 1 / 2 - sqrt(2) / 4]3-element Vector{Float64}:
 0.5
 0.8535533905932737
 0.1464466094067262
julia> b = [-1 / 2, 1 / 2 - sqrt(2) / 4, 1 / 2 + sqrt(2) / 4]3-element Vector{Float64}:
 -0.5
  0.1464466094067262
  0.8535533905932737
julia> c = [1 / sqrt(2), -1 / 2, 1 / 2]3-element Vector{Float64}:
  0.7071067811865475
 -0.5
  0.5
julia> # a×b = c
       a × b3-element Vector{Float64}:
  0.7071067811865475
 -0.5
  0.5
julia> # b×c = a
       b × c3-element Vector{Float64}:
 0.5
 0.8535533905932737
 0.14644660940672627
julia> # c×a = b
       c × a3-element Vector{Float64}:
 -0.5
  0.14644660940672627
  0.8535533905932737

ベクトル３重積

３つの空間ベクトルに対して，一般に，以下が成り立つ．

\[{a}\times ({b} \times {c}) = ({a}\cdot{c}){b} - ({a}\cdot {b}){c}\]

例： 具体的なベクトルの例で，上式が成り立つことを示せ．

julia> a = [-3, 1, 2]3-element Vector{Int64}:
 -3
  1
  2
julia> b = [2, -3, 1]3-element Vector{Int64}:
  2
 -3
  1
julia> c = [1, 2, -3]3-element Vector{Int64}:
  1
  2
 -3
julia> # 左辺
       lhs = a × (b × c)3-element Vector{Int64}:
  -7
  35
 -28
julia> # 右辺
       rhs = (a ⋅ c) * b - (a ⋅ b) * c3-element Vector{Int64}:
  -7
  35
 -28

▶ 行列の商（3x3 行列）

行列 A と行列（またはベクトル） Yに対して， 商 A\Y は，$AX-Y$の最小二乗ノルムが最小となる行列（またはベクトル） X を返す． 行列 A が正則なら，A の逆行列を 左から Y に乗じた行列（またはベクトル）と一致する．

julia> #
       a = [-3, 1, 2]3-element Vector{Int64}:
 -3
  1
  2
julia> b = [2, -3, 1]3-element Vector{Int64}:
  2
 -3
  1
julia> w = [a b]3×2 Matrix{Int64}:
 -3   2
  1  -3
  2   1
julia> c = [1, 2, -3]3-element Vector{Int64}:
  1
  2
 -3
julia> v = w \ c2-element Vector{Float64}:
 -1.0000000000000004
 -1.0000000000000002
julia> w * v - c3-element Vector{Float64}:
  8.881784197001252e-16
  4.440892098500626e-16
 -8.881784197001252e-16

▶ 行列の階数（ランク）

関数 rank(a) は，行列 a の階数（ランク，rank）を返す．

julia> a = [1 / 2, 1 / 2 + sqrt(2) / 4, 1 / 2 - sqrt(2) / 4]3-element Vector{Float64}:
 0.5
 0.8535533905932737
 0.1464466094067262
julia> b = [-1 / 2, 1 / 2 - sqrt(2) / 4, 1 / 2 + sqrt(2) / 4]3-element Vector{Float64}:
 -0.5
  0.1464466094067262
  0.8535533905932737
julia> c = [1 / sqrt(2), -1 / 2, 1 / 2]3-element Vector{Float64}:
  0.7071067811865475
 -0.5
  0.5
julia> v = [a b c]3×3 Matrix{Float64}:
 0.5       -0.5        0.707107
 0.853553   0.146447  -0.5
 0.146447   0.853553   0.5
julia> rank(v)3
julia> #
       a = [-3, 1, 2]3-element Vector{Int64}:
 -3
  1
  2
julia> b = [2, -3, 1]3-element Vector{Int64}:
  2
 -3
  1
julia> c = [1, 2, -3]3-element Vector{Int64}:
  1
  2
 -3
julia> v = [a b c]3×3 Matrix{Int64}:
 -3   2   1
  1  -3   2
  2   1  -3
julia> rank(v)2

▼ 行列の固有値・固有ベクトル

\[Ax = \lambda x\]

julia> a = [4 1; 2 3]2×2 Matrix{Int64}:
 4  1
 2  3
julia> #
       using LinearAlgebra
julia> # 固有値
       lam = eigvals(a)2-element Vector{Float64}:
 2.0
 5.0
julia> # 固有ベクトル
       evs = eigvecs(a)2×2 Matrix{Float64}:
 -0.447214  0.707107
  0.894427  0.707107
julia> evs[:, 1]2-element Vector{Float64}:
 -0.4472135954999579
  0.8944271909999159
julia> evs[:, 2]2-element Vector{Float64}:
 0.7071067811865475
 0.7071067811865475
julia> # 確認する a x - lambda x == 0 になるべき
       a * evs[:, 1] - lam[1] * evs[:, 1]2-element Vector{Float64}:
 0.0
 2.220446049250313e-16
julia> a * evs[:, 2] - lam[2] * evs[:, 2]2-element Vector{Float64}:
 0.0
 0.0

★ 今回のまとめ

  • ベクトルの内積
  • 行列の生成
  • 行列に対する関数
  • 行列とベクトルの演算
  • 行列と行列の演算
  • 部分行列
  • 2次元の回転行列