スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

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

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

小寺平治『明解演習線形代数』(共立出版) の 181 ページにある「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:

# 二次曲線の標準形を求める
normalizeQuadCurve := proc()
  local (
    f, A, b, c, s, coefxy,
    T, VAL, VEC, eigenvecs
  );

  f := args[1];
  coefxy := myCoeffTable(f, [x,y]);
  A := Matrix([
    [coefxy[2,0], coefxy[1,1]/2],
    [coefxy[1,1]/2, coefxy[0,2]]
  ]);
  b := <coefxy[1,0]/2, coefxy[0,1]/2>;
  c := coefxy[0,0];

  if LinearAlgebra[Rank](A) = 2 then
    # 平行移動
    s := LinearAlgebra[LinearSolve](A, -b);  # As + b = 0
    f := simplify(myVecSubs({<x,y> = <x,y> + s}, f));
  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> = T.<x,y>}, f));

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

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

    # 平行移動
    coefxy := myCoeffTable(f, [x,y]);
    if coefxy[0,1] <> 0 then
      f := simplify(subs({y = y - coefxy[0,0]/coefxy[0,1]}, f));
    end if;
  else
    # nothing to do
  end if;
  sort(f, [x,y]);
end proc:

f := x^2 - 4*x*y - 2*y^2 + 8*x + 20*y - 32;
normalizeQuadCurve(f);  #=> 双曲線

f := 5*x^2 - 6*x*y+5*y^2-4*x-4*y-4;
normalizeQuadCurve(f);  #=> 楕円

f := 9*x^2 + 24*x*y + 16*y^2 - 26*x + 7*y - 9;
normalizeQuadCurve(f);  #=> 放物線

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

参考文献

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

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

プロフィール

よしいず

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

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

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

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

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

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。