[CAS-Lab] 行列の固有値(座標変換と図形変形) Mathematica編

本稿は『数学アラカルト　行列の固有値とリーマン予想』で使用したMathematicaのソースコードの公開(個人での利用を許可します)と簡単な解説です。Mathematicaの環境が無い場合はこちらを参考に環境を整備してください。

サンプルプログラム実行前の準備

1) サンプルプログラム使用パッケージ（描画関数や座標変換関数等）のセーブ

以下ボタンで［.nb］ファイルをダウンロードしMathematicaで開いて実行してみてください。実行後いくつかのパッケージ（xxx.pac) がWindowsの場合デフォルトで”ドキュメント”　ディレクトリにセーブされるはずです。

(*Mathematica*)

(*---------------基本パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------------描画パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------------描画パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：

Save[(*cnd*)"udPlotPac1.pac",
 {
 セーブする関数、変数名のリスト  
 }];

サンプルプログラムの実行

以下のボタンでサンプルプログラムをダウンロードしMathematicaで開いて実行できます。あるいは、自身で新規の .nb ファイルを開き、以下のCodeをセル単位でコピペしながらステップバイステップで実行してみてください。

最初にグローバル変数等の初期化と上記のパッケージ “udPlotPac1.pac” と “udElemntPac1.pac” を読み込ます。これにより以降のプログラムで使用可能となります。

(* Mathematica *)
ClearAll["Global`*"]

Get["udPlotPac1.pac"];
Get["udElemntPac1.pac"];

以下各章節にリンクを貼りました。読み比べながら実行してみてください。

［1. ベクトル空間］の描画に対応

描画出力は代表の画像を1枚のみ貼り付けてあります。実行するとその他も描画されます。

［1.1 基底とベクトルの成分の関係］

(*グラフのフォーム設定*)
Gh = {{-5, 5}, {-5, 5}};
Gt = {{-5, 5}, {-5, 5}};
gX = 3 {{-1, 1}, {-1, 1}};
gU = 2 {{-1, 1}, {-1, 1}};
scText = {{1, ""}, {1, ""}};
(*{{Post Offset},{Pre Offset}}*) 
scOfst = {{0.2, -0.15}, {-0.2, 0.15}};

r = 2 {1, 1}
b = {{1, 0.25}, {0, 1}};
(*b={{1,1/4},{0.2,0.8}};*)
Print["B=", b // mf]
Print["\!\(\*SuperscriptBox[\(B\), \(-1\)]\)=", i[b] // mf]

s = i[b].r;
cdm = dm[s];
% // mf
vb = vec[{cGrnH, Thick}, t[b], {"b1", "b2"}, 1.4 lp[-22]];
vr = vec[{Red, Thick}, {r}, {"R"}, 1.7 lp[-22]];
br2 = (t[b].r)[[2]] t[b][[2]]
lbr2 = Graphics[{cBluH, Dashed, Line[{{0, 0}, br2, r}]}];
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluH, 1];
Show[grXsys, vb, vr, lbr2]
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, b, {2 {-1, 1}, 3 {-1, 1}}, scText, scOfst, 
   cGrnH, 1];
Show[grXsys, grUsys, vb, vr]
vbu = vec[{cGrnH, Thick}, t[i[b].b], {"b1", "b2"}, 1.4 lp[-22]];
vru = vec[{Red, 
    Thick}, {i[b].r}, {"S=\!\(\*SuperscriptBox[\(B\), \(-1\)]\)R"}, 
   1.7 lp[-22]];
grXsys = tfCrdSys2Dg[LTF, i[b], gX, scText, scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, e, {2 {-1, 1}, 3 {-1, 1}}, scText, scOfst, 
   cGrnH, 1];
Show[grXsys, grUsys, vbu, vru]


vr = vec[{Red, Thick}, {r}, {"R"}, 1.7 lp[-22]];
rb1 = t[b.cdm][[1]]; (*cdm=b^-1S成分の注入*)
rb2 = t[b.cdm][[2]];
(*rb=rb1+rb2;*)
rb = b.s;
rl = Graphics[{cBluH, Dashed, Line[{{rb1, rb}, {rb2, rb}}]}];
vrb = vec[{cBluH, Dashed}, {rb1, rb2, rb}, {"rb1", "rb2", "rb1+rb2"}];
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluH, 1];
Show[grXsys, vb, vr]
Show[grXsys, vb, vr, vrb, rl]
Show[grXsys, grUsys, vb, vr]

［1.2 座標変換と固有値の関係］に対応

回転変換の例

(*位置ベクトル p を計算する作業１*)
f[t_] := 2 N[{Cos[t], -Sin[t]}]
pltVf[cc_, txt_, w_: dm[{1, 1}]] := 
 Module[{p12, p, vf, vp, vp12, robj},
  p12 = w.dm[f[1.75 Pi]];
  vp12 = vec[{cc, Thick, Dashed}, t[p12], {txt[[1]], txt[[2]]}, 
    lp[-22]];
  p = t[p12][[1]] + t[p12][[2]];
  
  vp = vec[{cc, Thick}, {p}, {txt[[3]]}, lp[-22]];
  vf = Graphics[{cc, Thick, 
     Line[Table[w.f[t], {t, 0, 1.75 Pi, 0.01}]]}];
  robj = {vp12, vp, vf};
  robj
  ]
B = e;
L = dm[{1, 0.5}];

W = B; Print["W=", 
 W // mf]; Print["\!\(\*SuperscriptBox[\(W\), \(-1\)]\)=", i[W] // mf]

(*基準系を表示*)

grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluH, 1];
vW = vec[{cBluH, Thick}, t[B], {"e1", "e2"}];
pVfe = pltVf[Red, {" Px1", "Px2", "Px"}, B];
Show[grXsys, vW, 
 PlotRange -> {{Gh[[1, 1]], Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.7, 
    Gh[[2, 2]]}}]
Show[grXsys, vW, pVfe, 
 PlotRange -> {{Gh[[1, 1]], Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.7, 
    Gh[[2, 2]]}}]

(*基準系に回転ベクトルBを表示*)
B = rt[1/4 Pi].e;
(*B={{1,1/Sqrt[2]},{0,1/Sqrt[2]}};*)
W = B; Print["W=", 
 W // mf]; Print["\!\(\*SuperscriptBox[\(W\), \(-1\)]\)=", i[W] // mf]

vW = vec[{cGrnH, Thick}, t[B], {"   b1", "b2"}, 0.7 lp[-45]];
Show[grXsys, vW, pVfe, 
 PlotRange -> {{Gh[[1, 1]], Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.7, 
    Gh[[2, 2]]}}]

(*基準系に斜交ベクトルBと斜交座標と変換オブジェクトを表示*)
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, B, {2 {-1, 1}, 3 {-1, 1}}, scText, scOfst, 
   cGrnH, 1];
Show[grXsys, vW, grUsys, pVfe]


(*回転座標系に移動しベクトルBとと変換オブジェクトを表示*)
vW = vec[{cGrnH, Thick}, 
   t[i[B].B], {"\!\(\*SuperscriptBox[\(B\), \(-1\)]\)b1", 
    "       \!\(\*SuperscriptBox[\(B\), \(-1\)]\)b2"}, 0.7 lp[-45]];
grXsys = tfCrdSys2Dg[LTF, i[B], gX, scText, scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, i[B].B, {2 {-1, 1}, 3 {-1, 1}}, scText, 
   scOfst, cGrnH, 1];
Show[grXsys, grUsys, vW, 
 pltVf[Red, {"\!\(\*SuperscriptBox[\(B\), \(-1\)]\)Px1", 
   "\!\(\*SuperscriptBox[\(B\), \(-1\)]\)Px2", 
   "             \!\(\*SuperscriptBox[\(B\), \(-1\)]\)Px"}, i[B]]]

(*回転座標系に移動しLを乗じベクトルBとと変換オブジェクトを表示*)
grXsys = tfCrdSys2Dg[LTF, i[B], {3 {-1, 1}, 3 {-1, 1}}, scText, 
   scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, i[B].B, {2 {-1, 1}, 3 {-1, 1}}, scText, 
   scOfst, cGrnH, 1];
p = pltVf[
   Red, {"\!\(\*SuperscriptBox[\(LB\), \(-1\)]\)Px1           ", 
    "\!\(\*SuperscriptBox[\(LB\), \(-1\)]\)Px2         ", 
    "           \!\(\*SuperscriptBox[\(LB\), \(-1\)]\)Px"}, L.i[B]];
Show[grXsys, vW, grUsys, p[[1]], p[[2]], p[[3]]]



(*基準系に斜交ベクトルBと斜交座標と変換オブジェクトを表示*)
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluH, 1];
grUsys = tfCrdSys2Dg[LTF, B, {2 {-1, 1}, 3 {-1, 1}}, scText, scOfst, 
   cGrnL, 1];
p = pltVf[
   Red, {"", "", 
    "            \!\(\*SuperscriptBox[\(BLB\), \(-1\)]\)Px"}, 
   B.L.i[B]];
Show[grXsys, grUsys, p[[2]], p[[3]]]

斜交座標変換の例

B = e;
L = dm[{0.5, 2/Sqrt[2]}];
W = B; Print["W=", 
 W // mf]; Print["\!\(\*SuperscriptBox[\(W\), \(-1\)]\)=", i[W] // mf]

(*基準系を表示*)

grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluH, 1];
vW = vec[{cBluH, Thick}, t[B], {"e1", "e2"}];
pVfe = pltVf[Red, {" Px1", "Px2", "Px"}, B];
Show[grXsys, vW, 
 PlotRange -> {{Gh[[1, 1]], Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.7, 
    Gh[[2, 2]]}}]
Show[grXsys, vW, pVfe, 
 PlotRange -> {{Gh[[1, 1]], Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.7, 
    Gh[[2, 2]]}}]

(*基準系に斜交ベクトルBを表示*)
B = {{1, 1}, {0, 1}};
(*B={{1,1/Sqrt[2]},{0,1/Sqrt[2]}};*)
W = B; Print["W=", 
 W // mf]; Print["\!\(\*SuperscriptBox[\(W\), \(-1\)]\)=", i[W] // mf]
pVfe = pltVf[RGBColor[1, 0, 0, 0.3], {" Px1", "Px2", "Px"}, e];
vW = vec[{cGrnH, Thick}, t[B], {"b1", "b2             "}, 0.7 lp[-45]];
Show[grXsys, vW, pVfe, 
 PlotRange -> {{Gh[[1, 1]], Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.7, 
    Gh[[2, 2]]}}]

(*基準系に斜交ベクトルBと斜交座標と変換オブジェクトを表示*)
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, B, {2 {-1, 1}, 3 {-1, 1}}, scText, scOfst, 
   cGrnH, 1];
Show[grXsys, vW, grUsys, pVfe]


(*斜交座標系に移動しベクトルBとと変換オブジェクトを表示*)
vW = vec[{cGrnH, Thick}, 
   t[i[B].B], {"\!\(\*SuperscriptBox[\(B\), \(-1\)]\)b1", 
    "\!\(\*SuperscriptBox[\(B\), \(-1\)]\)b2"}, 0.7 lp[-45]];
grXsys = tfCrdSys2Dg[LTF, i[B], gX, scText, scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, i[B].B, {2 {-1, 1}, 3 {-1, 1}}, scText, 
   scOfst, cGrnH, 1];
Show[grXsys, grUsys, vW, 
 pltVf[Red, {"\!\(\*SuperscriptBox[\(B\), \(-1\)]\)Px1", 
   "\!\(\*SuperscriptBox[\(B\), \(-1\)]\)Px2", 
   "             \!\(\*SuperscriptBox[\(B\), \(-1\)]\)Px"}, i[B]]]

(*斜交座標系に移動しLを乗じベクトルBとと変換オブジェクトを表示*)
grXsys = tfCrdSys2Dg[LTF, i[B], {3 {-1, 1}, 3 {-1, 1}}, scText, 
   scOfst, cBluL, 1];
grUsys = tfCrdSys2Dg[LTF, i[B].B, {2 {-1, 1}, 3 {-1, 1}}, scText, 
   scOfst, cGrnH, 1];
p = pltVf[
   Red, {"\!\(\*SuperscriptBox[\(LB\), \(-1\)]\)Px1           ", 
    "\!\(\*SuperscriptBox[\(LB\), \(-1\)]\)Px2         ", 
    "           \!\(\*SuperscriptBox[\(LB\), \(-1\)]\)Px"}, L.i[B]];
Show[grXsys, vW, grUsys, p[[1]], p[[2]], p[[3]]]



(*基準系に斜交ベクトルBと斜交座標と変換オブジェクトを表示*)
grXsys = tfCrdSys2Dg[LTF, e, gX, scText, scOfst, cBluH, 1];
grUsys = tfCrdSys2Dg[LTF, B, {2 {-1, 1}, 3 {-1, 1}}, scText, scOfst, 
   cGrnL, 1];
p = pltVf[
   Red, {"", "", 
    "            \!\(\*SuperscriptBox[\(BLB\), \(-1\)]\)Px"}, 
   B.L.i[B]];
Show[grXsys, grUsys, p[[2]], p[[3]]]

［1.3 時間-周波数領域変換］に対応

アダマール変換

n = 8

B = HadamardMatrix[n];
A = i[B];
(*Print["B=",B//mf,",  B^-1=",i[B]//Simplify//mf];*)
vp = Table[
   DiscretePlot[t[B][[iy, ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2,
     AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-0.5, 0, 0.5}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 0.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, n}];
GraphicsColumn[vp, Frame -> All]

(*参照信号Rの生成*)
R = Table[Cos[2 Pi 1/n t], {t, 0, n - 1}];
(*R=Table[Sin[2Pi 1/n t],{t,0,n-1}];*)

vp = Table[
   DiscretePlot[R[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 1.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]


S = A.R;
vp = Table[
   DiscretePlot[S[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 2 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]

Rd = B.S;
vp = Table[
   DiscretePlot[Rd[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 1.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]

離散コサイン変換(DCT)

(*離散コサイン変換（DCT）*)
n = 8
B = N[FourierDCTMatrix[n]];
A = i[B];
(*Print["B=",B//mf,",  B^-1=",i[B]//Simplify//mf];*)
vp = Table[
   DiscretePlot[t[B][[iy, ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2,
     AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, -0.5, 0, 0.5, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 0.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, n}];
GraphicsColumn[vp, Frame -> All]

(*GraphicsGrid[Table[plotCpSpc[B[[ix,iy]],1/n^0.5],{ix,1,n},{iy,1,n}],\
Frame\[Rule]All]
*)
R = Table[Cos[2 Pi 1/n t], {t, 0, n - 1}];
vp = Table[
   DiscretePlot[R[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 1.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]

S = A.R;
vp = Table[
   DiscretePlot[S[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 3 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]

Rd = B.S;
vp = Table[
   DiscretePlot[Rd[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 1.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]

(*lv=Graphics[{cBluL,rplot[1/3S]}];
grXsys=tfCrdSys2Dg[LTF,e,{{0,8},{-1,1}},{{1,""},{3,""}}
,scOfst,cBluH,1];
Show[{grXsys,lv}]*)

離散フーリエ変換（DFT)

B = FourierMatrix[n];
Print[1/Sqrt[n], Sqrt[n] B // mf]
A = i[B];
GraphicsGrid[
 Table[plotCpSpc[B[[ix, iy]], 1/n^0.5], {ix, 1, n}, {iy, 1, n}], 
 Frame -> All]

GraphicsGrid[
 Table[plotCpSpc[A[[ix, iy]], 1/n^0.5], {ix, 1, n}, {iy, 1, n}], 
 Frame -> All]

(*参照信号Rの生成 Cos*)
R = Table[Cos[2 Pi 1/n t], {t, 0, n - 1}];
GraphicsRow[Table[plotCpSpc[R[[ix]], 1], {ix, 1, n}], Frame -> All]
S = A.R;
GraphicsRow[Table[plotCpSpc[S[[ix]], 1], {ix, 1, n}], Frame -> All]
Rd = B.S;
GraphicsRow[Table[plotCpSpc[Rd[[ix]], 1], {ix, 1, n}], Frame -> All]
(*参照信号Rの生成 Sin*)
R = Table[Sin[2 Pi 1/n t], {t, 0, n - 1}];
vp = Table[
   DiscretePlot[R[[ix + 1]], {ix, 0, n - 1}, ExtentSize -> 0.2, 
    AspectRatio -> 0.125, AxesOrigin -> {0, 0}, 
    Ticks -> {Range[0, n - 1], {-1, 0, 1}}, 
    PlotRange -> {{-1, n - 1 + 0.5}, 1.5 {-1, 1}}, 
    PlotStyle -> {cBluL, PointSize[0.02]}], {iy, 1}];
GraphicsColumn[vp, Frame -> All]
GraphicsRow[Table[plotCpSpc[R[[ix]], 1], {ix, 1, n}], Frame -> All]
S = A.R;
GraphicsRow[Table[plotCpSpc[S[[ix]], 1], {ix, 1, n}], Frame -> All]
Rd = B.S;
GraphicsRow[Table[plotCpSpc[Rd[[ix]], 1], {ix, 1, n}], Frame -> All]

以上

（時間が出来たら、もう少し補足・説明を追加します）