Home › Category Archives › 計算

(18) 行列計算で中1レベル連立方程式を解く

2元連立1次方程式をOctaveを使って解いてみる。

まずは話を簡単にするために、以下のような簡単な2式の連立方程式を解いてみる。
y = 2x + 4
y = -2x + 16

Octaveでグラフ表示すると以下のような直線が現れる。

octave:21> x=[-10:10]
x =
  -10   -9   -8   -7   -6   -5   -4   -3   -2   -1    0    1    2    3    4    5    6    7    8    9   10

octave:22> y1 = 2 * x + 4
y1 =
  -16  -14  -12  -10   -8   -6   -4   -2    0    2    4    6    8   10   12   14   16   18   20   22   24

octave:23> y2 = -2 * x + 16
y2 =
   36   34   32   30   28   26   24   22   20   18   16   14   12   10    8    6    4    2    0   -2   -4

octave:24> plot(x,y1,x,y2)
octave:25>

まずは与えられた2式を行列で表せるように各項を並び替える。
y = 2x + 4
y = -2x + 16
 ↓
2x – y = -4
-2x – y = -16

これを以下のような行列式で表し、両辺に逆行列をかけて解を求める。

これをOctaveで書くと以下のようになる。

octave:24> A=[2 -1 ; -2 -1]
A =
   2  -1
  -2  -1

octave:25> B=[-4 ; -16]
B =
   -4
  -16

octave:26> Ainv=inv(A)
Ainv =
   0.25000  -0.25000
  -0.50000  -0.50000

octave:27> xy=Ainv * B
xy =
    3
   10

octave:28>

x=3, y=10 が求まった。

※注意
この方法はいつでも使えるわけではない。
行列Aが逆行列を持たない場合、すなわち両方の式の傾きが同じ場合、2本の直線は交点を持たず、連立方程式として成り立たない。
これをOctaveで確認してみる。

octave:34> A=[1 2;4 8]
A =
   1   2
   4   8

octave:35> Ainv=inv(A)
warning: inverse: matrix singular to machine precision, rcond = 0
Ainv =
   Inf   Inf
   Inf   Inf

octave:36>

逆行列が求まらなかった。
2元1次方程式であればイメージがわきやすいが、3元,4元,…,多元連立方程式になると直感的に逆行列の存在有無を判断できない。

次に以下のような4元1次連立方程式を解いてみる。
a + 2b + 3c + 4d = 8
5a + 4b + 6c + 7d = 3
6a + 3b + 4c + 9d = 9
5a + 7b + c + 5d = 5

octave:37> A=[1 2 3 4;5 4 6 7;6 3 4 9;5 7 1 5]
A =
   1   2   3   4
   5   4   6   7
   6   3   4   9
   5   7   1   5

octave:38> B=[8;3;9;5]
B =
   8
   3
   9
   5

octave:39> i=inv(A)
i =
  -5.5674e-01   2.3404e-01   6.7376e-02  -3.5461e-03
   1.6667e-01  -2.5967e-17  -1.6667e-01   1.6667e-01
  -3.1915e-02   3.1915e-01  -1.8085e-01  -9.5745e-02
   3.2979e-01  -2.9787e-01   2.0213e-01  -1.0638e-02

octave:40> A*i
ans =
   1.00000   0.00000   0.00000  -0.00000
   0.00000   1.00000  -0.00000  -0.00000
   0.00000   0.00000   1.00000  -0.00000
   0.00000   0.00000   0.00000   1.00000

octave:41> i*B
ans =
  -3.16312
   0.66667
  -1.40426
   3.51064

octave:42>

a=-3.16312, b=0.66667, c=-1.40426, d=3.51064 が求まった。

★確認
逆行列と元の行列の積は1(=単位行列)となる。上では「octave:40」の A*i でこれを確認している。

(15) ベクトルの内積

Octaveではベクトルの内積も簡単に求められる。
ベクトルaとベクトルbを以下のように定義する。

octave:20> a=[1 2 3]
a =
   1   2   3
octave:21> b=[2 3 4]
b =
   2   3   4

内積は「a・b」だからと「a*b」と書いてしまうとダメ。

octave:22> a*b
error: operator *: nonconformant arguments (op1 is 1x3, op2 is 1x3)

1行3列の行列の積として認識されてしまうのでエラーが出てしまう。

そこで、行列の積が成り立つように形を変えてやる。
すなわち、aが1行3列であれば、bを3行1列に転置してやれば行列の積が成り立つ。
Octaveでは行列の転置を「’」で記述する。

octave:22> b'
ans =
   2
   3
   4
octave:23> a*b'
ans =  20

a=(1,2,3) と b=(2,3,4) の内積は、1×2+2×3+3×4=20 である。

当然だが、以下のように要素の積の総和を求めても同じ。

octave:24> sum(a.*b)
ans =  20

(3) Octaveで行列計算 (足し算)

Octaveを使うと行列計算が簡単にできる。
個人的にはExcelよりも便利だと思う。

まずは足し算

octave:2> a=[1 2 3;4 5 6;7 8 9]
a =
   1   2   3
   4   5   6
   7   8   9

octave:3> b=[5 6 7;8 9 10;11 12 13]
b =
    5    6    7
    8    9   10
   11   12   13

octave:4> a+b
ans =
    6    8   10
   12   14   16
   18   20   22

octave:5>

ベクトルデータとスカラーデータの足し算もできる。

octave:5> a=[1 2 3;4 5 6;7 8 9]
a =
   1   2   3
   4   5   6
   7   8   9

octave:6> a+1
ans =
    2    3    4
    5    6    7
    8    9   10

octave:7>

これは以下のような書き方と同じ結果となる。

octave:7> a=[1 2 3;4 5 6;7 8 9]
a =
   1   2   3
   4   5   6
   7   8   9

octave:8> b=ones(3,3)
b =
   1   1   1
   1   1   1
   1   1   1

octave:9> a+b
ans =
    2    3    4
    5    6    7
    8    9   10

octave:10>

以下のような書き方でも同じ結果となる。(warningが出るけれども…)

octave:12> a=[1 2 3;4 5 6;7 8 9]
a =
   1   2   3
   4   5   6
   7   8   9

octave:13> b=[1 1 1]
b =
   1   1   1

octave:14> a+b
warning: operator +: automatic broadcasting operation applied
ans =
    2    3    4
    5    6    7
    8    9   10

octave:15> b=[1;1;1]
b =
   1
   1
   1

octave:16> a+b
warning: operator +: automatic broadcasting operation applied
ans =
    2    3    4
    5    6    7
    8    9   10

octave:17>