二次曲面の標準形を求めるMapleプログラムの例

二次曲面の標準形を求める Maple プログラムの例。

小寺平治『明解演習線形代数』(共立出版) の 182 ページにある「2次曲面の主軸変換の手順」をプログラム化しました。

# 多変数多項式の係数の抽出
myCoeffTable := proc()
  local (
    f,     # 与えられた多項式
    coef,  # 係数
    var,   # 変数
    dim,   # 変数の数
    deg,   # 多項式の次数
    i, j, k, B, T, nv, nvseq, c
  );

  f := args[1];
  var := args[2];
  dim := nops(var);
  deg := degree(f);

  B := [seq([seq(i, i=0..deg)], i=1..dim)];
  T := combinat[cartprod](B);  # [0,..,0]→…→[deg,..,deg]
  while not T[finished] do
    nv := ListTools[Reverse](T[nextvalue]());
    nvseq := seq(nv[i], i=1..nops(nv));
    c[0] := f;
    for j from 1 to dim do
      c[j] := coeff(c[j-1], var[j], nv[j]);
    end do;
    coef[nvseq] := c[dim];
  end do:
  coef;
end proc:

# 列ベクトルの成分ごとに式を代入する
myVecSubs := proc()
  local f, v, w, vdim, wdim, i;

  f := args[2];
  w := lhs(args[1][1]);
  wdim := LinearAlgebra[Dimension](w);
  v := rhs(args[1][1]);
  vdim := LinearAlgebra[Dimension](v);
  if wdim != vdim then
    error "両辺の列ベクトルの次元が一致しません。";
  end if;
  subs({seq(w[i]=v[i], i=1..wdim)}, f)
end proc:

# 二次曲面の標準形を求める
normalizeQuadSurface := proc(arg)
  local (
    f, A, b, c, s, coefxyz,
    T, VAL, VEC, eigenvecs,
    R, cs, sn
  );
  f := args[1];
  coefxyz := myCoeffTable(f, [x,y,z]);

  A := Matrix([
    [coefxyz[2,0,0], coefxyz[1,1,0]/2, coefxyz[1,0,1]/2],
    [coefxyz[1,1,0]/2, coefxyz[0,2,0], coefxyz[0,1,1]/2],
    [coefxyz[1,0,1]/2, coefxyz[0,1,1]/2, coefxyz[0,0,2]]
  ]);
  b := <coefxyz[1,0,0]/2, coefxyz[0,1,0]/2, coefxyz[0,0,1]/2>;
  c := coefxyz[0,0,0];

  if LinearAlgebra[Rank](A) = 3 then
    # 平行移動
    s := LinearAlgebra[LinearSolve](A, -b);  # As + b = 0
    f := simplify(myVecSubs({<x,y,z> = <x,y,z> + s}, f));
  elif LinearAlgebra[Rank](A) = 2 then
    # nothing to do
  elif LinearAlgebra[Rank](A) = 1 then
    # nothing to do
  else
    print("Error");
  end if;

  # 回転
  VAL, VEC := LinearAlgebra[Eigenvectors](A);
  eigenvecs := LinearAlgebra[Column](VEC,
    [1 .. LinearAlgebra[ColumnDimension](VEC)]);
  T := Matrix(LinearAlgebra[GramSchmidt]([eigenvecs], normalized));
  f := simplify(myVecSubs({<x,y,z> = T.<x,y,z>}, f));

  if LinearAlgebra[Rank](A) = 2 then
    if degree(f, y) = 2 and degree(f, z) = 2 then
      f := subs({x = z, z = x}, f);
    elif degree(f, x) = 2 and degree(f, z) = 2 then
      f := subs({y = z, z = y}, f);
    end if;

    # 平行移動
    coefxyz := myCoeffTable(f, [x,y,z]);
    f := simplify(subs({x = x - coefxyz[1,0,0]/(2*coefxyz[2,0,0])}, f));
    f := simplify(subs({y = y - coefxyz[0,1,0]/(2*coefxyz[0,2,0])}, f));

    # 平行移動
    coefxyz := myCoeffTable(f, [x,y,z]);
    if coefxyz[0,0,1] != 0 then
      f := simplify(subs({z = z - coefxyz[0,0,0]/coefxyz[0,0,1]}, f));
    end if;
  elif LinearAlgebra[Rank](A) = 1 then
    if degree(f, y) = 2 then
      f := subs({x = y, y = x}, f);
    elif degree(f, z) = 2 then
      f := subs({x = z, z = x}, f);
    end if;

    # 平行移動
    coefxyz := myCoeffTable(f, [x,y,z]);
    f := simplify(subs({x = x - coefxyz[1,0,0]/(2*coefxyz[2,0,0])}, f));  

    # 平行移動
    coefxyz := myCoeffTable(f, [x,y,z]);
    if coefxyz[0,0,1] != 0 then
      f := simplify(subs({z = z - coefxyz[0,0,0]/coefxyz[0,0,1]}, f));
    end if;

    # 回転
    coefxyz := myCoeffTable(f, [x,y,z]);
    if coefxyz[0,0,1] != 0 then
      cs := coefxyz[0,1,0]/sqrt(coefxyz[0,1,0]^2+coefxyz[0,0,1]^2);
      sn := coefxyz[0,0,1]/sqrt(coefxyz[0,1,0]^2+coefxyz[0,0,1]^2);
      R := Matrix([[1, 0, 0], [0, cs, -sn], [0, sn, cs]]);
      f := simplify(myVecSubs({<x,y,z> = R.<x,y,z>}, f));
    end if;
  else
       # nothing to do
  end if;
  sort(f, [x, y, z]);
end proc:

f := 2*x^2 - 3*y^2 + z^2 - 8*y*z - 12*z*x + 4*x*y + 20*x + 12*y + 8*z - 36;
normalizeQuadSurface(f);  #=> 2葉双曲面

f := x^2 + 2*y^2 + z^2 - 2*y*z - 2*x*y - 2*x - 6*y + 2*z + 6;
normalizeQuadSurface(f);  #=> 楕円的放物面

f := x^2 + 4*y^2 + 4*z^2 + 8*y*z + 4*z*x + 4*x*y + 6*x - 10*y - 8*z + 1;
normalizeQuadSurface(f);  #=> 放物柱面

f := x^2 + 4*y^2 + 4*z^2 + 8*y*z + 4*z*x + 4*x*y - 2*x - 4*y - 4*z + 1;
normalizeQuadSurface(f);  #=> 2重1平面

f := y*z + z*x + x*y - x - 2*y - 3*z + 2;
normalizeQuadSurface(f);  #=> 錘面

※ プロシージャを実行するたびに、標準形における係数 (A の固有値) が入れ替わる場合がある。これは、LinearAlgebra[Eigenvectors] の出力結果が一定でないことによる。

参考文献

  • 小寺平治:明解演習線形代数, 共立出版, 1982.

【theme : 数学
【genre : 学問・文化・芸術

プロフィール

よしいず

Author:よしいず
MATHEMATICS.PDFというウェブサイトを運営しています。

管理の都合上、トラックバックとコメントはオフにしてあります。ブログ経験者なら分かっていただけると思いますが、スパム(アダルトやその他の宣伝)ばかりなのが現実です。

リンクは自由です。当サイトの記事に対する間違いの指摘・意見・感想などを述べた記事からのリンクは歓迎です。ただし、ブログ記事アップ直後はミスが多く、頻繁に修正します。場合によっては削除する可能性もあります。その際、何も断りもなく修正・削除しますがご了承ください。内容を参考にする場合には投稿後一週間ほど様子を見てからにしてください(笑)。

記事の間違いを指摘するときは、その具体的箇所、理由(仕様に反するなど)・根拠(参考にした文献など)、代替案(同じ結果を得るための正しいやり方)も教えてください。そうしないと、(指摘される側および第三者はその時点では無知の状態なので、)どこが間違いなのか分かりませんし、本当に間違っているのかどうかが判断・検証できません。実際、間違いだと指摘されたことが結局は正しかったというケースもありますので。

このブログのタイトル一覧

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示