【数式処理 MuPAD 入門】
平田 浩一 / 2000年01月31日 
hirata@edserv.ed.ehime-u.ac.jp 
愛媛大学教育学部数学教室 


行列とベクトル (3)



ベクトル

今回はベクトルの計算について解説します。

まずは、行列を扱うための準備をしましょう。

● loadlib("linalg"): export(linalg);

● MAT := Dom::Matrix();

               Dom::Matrix(Dom::ExpressionField(id, iszero))

■ ベクトルの計算

ベクトルの計算のために次のような関数が用意されています。これらの関数は行ベクトルにも列ベクトルにも用いることができます。

norm(u, 2)     ベクトル u の大きさ
normalize(u)     ベクトル u と同じ向きで長さ 1 のベクトル
scalarProduct(u, v)     ベクトル u, v の内積 u・v
angle(u, v)     ベクトル u, v のなす角 (ラジアン)

● u := MAT([[1,2]]);

                                 +-    -+
                                 | 1, 2 |
                                 +-    -+

● v := MAT([[2,3]]);

                                 +-    -+
                                 | 2, 3 |
                                 +-    -+

● norm(u,2);

                                    1/2
                                   5

● normalize(u);

                            +-              -+
                            |   1/2     1/2  |
                            |  5     2 5     |
                            |  ----, ------  |
                            |   5      5     |
                            +-              -+

● scalarProduct(u,v);

                                     8

● a := angle(u,v);

                               /    1/2   1/2 \
                               | 8 5    13    |
                           acos| ------------ |
                               \      65      /

ベクトル u, v のなす角の計算は、 ほとんどの場合上のような式になってしまい、 きちんとした値がでてきません。 このようなときは、近似計算をしましょう。

● float(a);

                               0.1243549945

単位はラジアンです。

ラジアンではなく度で求めたいときは次のようにします。

● float(a*180/PI);

                                7.125016349

以上は、2次元ベクトルだけでなく、一般のn次元ベクトルに対しても計算できます。

3次元ベクトルに対しては、外積を計算できます。

crossProduct(u, v)     ベクトル u, v の外積 u×v

● u := MAT([[1,0,1]]);

                                +-       -+
                                | 1, 0, 1 |
                                +-       -+

● v := MAT([[0,1,2]]);

                                +-       -+
                                | 0, 1, 2 |
                                +-       -+

● crossProduct(u,v);

                               +-         -+
                               | -1, -2, 1 |
                               +-         -+
1次変換

次は1次変換です。 曲線を一次変換してえられるグラフと方程式を取り上げてみましょう。

■ パラメタ表示 x = f(t), y = g(t)

例として、パラメタ表示された曲線 x = cos3t, y = sin3t  (0 ≦ t ≦ 2π) を取り上げます。

まずはグラフを書いてみましょう。 曲線はパラメタ t を含むベクトルとして入力します。

● c := MAT([cos(t)^3, sin(t)^3]);

                               +-         -+
                               |        3  |
                               |  cos(t)   |
                               |           |
                               |        3  |
                               |  sin(t)   |
                               +-         -+

● plot2d(Axes=Box, [Mode=Curve, [c[1], c[2]], t=[0,2*PI], Grid=[50]]);

関数 plot2d 内の c[1] と c[2] は、 ベクトル c の第1成分と第2成分を表しています。

グラフ1

行列 A := MAT([[1, 2], [3, 1]]) による1次変換でグラフを移動してみましょう。

● A := MAT([[1, 2], [3, 1]]);

                                +-      -+
                                |  1, 2  |
                                |        |
                                |  3, 1  |
                                +-      -+

● c1 := A * c;

                         +-                     -+
                         |        3           3  |
                         |  cos(t)  + 2 sin(t)   |
                         |                       |
                         |          3         3  |
                         |  3 cos(t)  + sin(t)   |
                         +-                     -+

● plot2d(Axes=Box, [Mode=Curve, [c1[1], c1[2]], t=[0,2*PI], Grid=[50]]);

グラフ2

■ 方程式 f(x, y) = 0

次に、曲線が方程式 f(x, y) = 0 で表されている場合の1次変換を考えてみましょう。

方程式は f(x, y) = (x2+y2)3 - (x2-y2)2 = 0 とし、これを π/4 ( 45°) 回転してみましょう。

変換する前のグラフを求めておきましょう。

● f := func((x^2+y^2)^3 - (x^2-y^2)^2, x, y);

                 func((x^2 + y^2)^3 - (x^2 - y^2)^2, x, y)

● plotlib::implicitplot(f, -1..1, -1..1, 7):

描画されるまで数分かかります。

グラフ3

1次変換は a = π/4 とおくとき、

     +-  -+     +-              -+  +- -+
     | x1 |     | cos a   -sin a |  | x |
     |    |  =  |                |  |   |
     | y1 |     | sin a    cos a |  | y |
     +-  -+     +-              -+  +- -+

となります。これを linearSolve で解きましょう。

● a := PI/4;

                                    PI
                                    --
                                    4

● A := MAT([[cos(a),-sin(a)],[sin(a),cos(a)]]);

                            +-              -+
                            |   1/2     1/2  |
                            |  2       2     |
                            |  ----, - ----  |
                            |   2       2    |
                            |                |
                            |   1/2    1/2   |
                            |  2      2      |
                            |  ----,  ----   |
                            |   2      2     |
                            +-              -+

● V1 := MAT([x1, y1]);

                                 +-    -+
                                 |  x1  |
                                 |      |
                                 |  y1  |
                                 +-    -+

● v := linearSolve(A,V1);

                          +-                   -+
                          |    1/2 / x1   y1 \  |
                          |   2    | -- + -- |  |
                          |        \ 2    2  /  |
                          |                     |
                          |   1/2               |
                          |  2    (- x1 + y1 )  |
                          |  -----------------  |
                          |          2          |
                          +-                   -+

これを関数 f に代入し、展開します。

● g := f(v[1],v[2]);

    /          2                  \3   /                           2 \2
    | (y1 - x1)      / x1   y1 \2 |    |   / x1   y1 \2   (y1 - x1)  |
    | ---------- + 2 | -- + -- |  |  - | 2 | -- + -- |  - ---------- |
    \     2          \ 2    2  /  /    \   \ 2    2  /        2      /

● expand(g);

                 6     6       2   2       2   4       4   2
               x1  + y1  - 4 x1  y1  + 3 x1  y1  + 3 x1  y1

これで、f(x, y) = 0 を π/4 回転したときの方程式 g(x, y) = 0 の左辺がもとまりました。

グラフを求めてみましょう。関数 g の式は手で直接入力します。

● g := func(x1^6 + y1^6 - 4*x1^2*y1^2 + 3*x1^2*y1^4 + 3*x1^4*y1^2, x1, y1);

  func(((x1^6 + y1^6) - 4*x1^2*y1^2) + 3*x1^2*y1^4 + 3*x1^4*y1^2, x1, y1)

● plotlib::implicitplot(g, -1..1, -1..1, 7);

グラフ4
固有値と固有ベクトル

次に、固有値と固有ベクトルの計算です。 固有値と固有ベクトルについては既知と仮定します。 詳しいことは線形代数の教科書を参考ください。

■ 固有値

行列 A の固有値 t は、方程式 det(t E - A) = 0 の解です。 ここで、E は単位行列です。

行列 A := MAT([[2, 3], [1, 4]]) の固有値を求めてみましょう。

● A := MAT([[2,3],[1,4]]);

                                +-      -+
                                |  2, 3  |
                                |        |
                                |  1, 4  |
                                +-      -+

● E2 := MAT(2,2,1,Diagonal);

                                +-      -+
                                |  1, 0  |
                                |        |
                                |  0, 1  |
                                +-      -+

● cp := det(t * E2 - A);

                            (t - 2) (t - 4) - 3

● solve(cp = 0, t);

                                  {1, 5}

この計算にでてくる、行列 t E - A を 特性行列 といいます。 またその行列式 det(t E - A) を 特性多項式 といいます。

これらを簡単に計算する関数が用意されています。

charMatrix(A, t)     行列 A の特性行列
charPolynomial(A, t)     行列 A の特性多項式
eigenValues(A)     行列 A の固有値

● charMatrix(A, t);

                            +-              -+
                            |  t - 2,   -3   |
                            |                |
                            |    -1,  t - 4  |
                            +-              -+

● charPolynomial(A, t);

                                       2
                              - 6 t + t  + 5

● eigenValues(A);

                                  {1, 5}

■ 固有ベクトル

固有ベクトルを求める関数が用意されています。

eigenVectors(A)     行列 A の固有値、重複度、固有ベクトル

● A := MAT([[2,3],[1,4]]);

                                +-      -+
                                |  2, 3  |
                                |        |
                                |  1, 4  |
                                +-      -+

● eigenVectors(A);

        -- --       -- +-    -+ -- --  --       -- +-   -+ -- -- --
        |  |        |  |  -3  |  |  |  |        |  |  1  |  |  |  |
        |  |  1, 1, |  |      |  |  |, |  5, 1, |  |     |  |  |  |
        |  |        |  |   1  |  |  |  |        |  |  1  |  |  |  |
        -- --       -- +-    -+ -- --  --       -- +-   -+ -- -- --

関数 eigenVectors が計算して求めた結果は、
  (1) 固有値 1、重複度 1、固有ベクトル [-3, 1]
  (2) 固有値 5、重複度 1、固有ベクトル [1, 1]
です。

もう1つ計算してみましょう。

● B := MAT([[3, 2], [0, 3]]);

                                +-      -+
                                |  3, 2  |
                                |        |
                                |  0, 3  |
                                +-      -+

● eigenVectors(B);

                      -- --       -- +-   -+ -- -- --
                      |  |        |  |  1  |  |  |  |
                      |  |  3, 2, |  |     |  |  |  |
                      |  |        |  |  0  |  |  |  |
                      -- --       -- +-   -+ -- -- --

計算結果は、
  (1) 固有値 3、重複度 2、固有ベクトル [1, 0]
です。

■ ジョルダン標準形

ジョルダン標準形を求める関数は jordanForm です。

jordanForm(A)     行列 A のジョルダン標準形

● A := MAT([[2,3],[1,4]]);

                                +-      -+
                                |  2, 3  |
                                |        |
                                |  1, 4  |
                                +-      -+

● jordanForm(A);

                                +-      -+
                                |  1, 0  |
                                |        |
                                |  0, 5  |
                                +-      -+

【演習問題】

(1) ベクトル u = [3, 5], v = [1, -2] について、次を計算せよ。

  1.  u, v の大きさ。
  2.  内積 u・v の値。
  3.  u, v のなす角のラジアンでの近似値。

(2) ベクトル u = [3, 5, 1], v = [-1, 1, -2] について、次を計算せよ。

  1.  u, v の大きさ。
  2.  内積 u・v の値。
  3.  外積 u×v。

(3) 曲線 x = 3 cos t + cos 3t, y = 3 sin t - sin 3t (0 ≦ t ≦ 2π) を 原点のまわりに 60°回転させたときの曲線のグラフを求めよ。

(4) 2次関数 y = x2 + x を原点のまわりに 30°回転させた 曲線を C とする。 曲線 C の方程式を求め、グラフに表せ。

(5) 次の行列の固有値、固有ベクトル、ジョルダン標準形を求めよ。

           +-    -+
           | 1  2 |
      1.   |      |
           | 4  3 |
           +-    -+


           +-            -+
           |   6  -1   5  |
           |              |
      2.   |  -3   2  -3  |
           |              |
           |  -7   1  -6  |
           +-            -+

【数式処理 MuPAD 入門】