📄 polynomialmatrices.nb
字号:
a[\([k]\)] \[Rule]
a[\([rmin]\)]}]; \[IndentingNewLine]For[i = k + 1,
i \[LessEqual] p, \(i++\),
If[a[\([i, kcl]\)] =!= 0,
lpc = PolynomialQuotient[a[\([i, kcl]\)],
a[\([k, kcl]\)], x]; \[IndentingNewLine]a =
ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc],
i]]]]; (*Endwhile*) \[IndentingNewLine]If[
a[\([k, kcl]\)] =!= 0,
For[i = 1, i < k, \(i++\),
If[a[\([i, kcl]\)] =!= 0,
While[Exponent[a[\([i, kcl]\)], x] \[GreaterEqual]
Exponent[a[\([k, kcl]\)], x],
lpc = PolynomialQuotient[a[\([i, kcl]\)],
a[\([k, kcl]\)], x]; \[IndentingNewLine]a =
ReplacePart[a, Together[a[\([i]\)] - a[\([k]\)]*lpc],
i];] (*endwhile*) ]]]; \[IndentingNewLine]If[
a[\([k, kcl]\)] =!= 0,
lpc = Coefficient[a[\([k, kcl]\)], x,
Exponent[a[\([k, kcl]\)], x]]; \[IndentingNewLine]If[
lpc =!= 1,
a = ReplacePart[a, Expand[a[\([k]\)]/lpc],
k]]]; \[IndentingNewLine]\(kcl++\)]; (*Endfor*) For[
k = 1, k < p, \(k++\),
If[a[\([k]\)] === Table[0, {m}],
i = p; \[IndentingNewLine]While[
a[\([i]\)] === Table[0, {m}] && i > k,
i = i - 1]; \[IndentingNewLine]If[i > k,
For[j = k, j \[LessEqual] i, \(j++\),
a = ReplacePart[a, a[\([j + 1]\)],
j]]];] (*endif*) ]; (*endfor*) If[ch,
a = a /. {x^\((n : _)\) \[Rule] x^\(-n\),
x \[Rule] x^\(-1\)}]; \[IndentingNewLine]Return[a]] /;
MatrixQ[at] ||
Message[General::mtrx, "\<HermiteForm\>"];\)\)], "Input",
InitializationCell->True],
Cell[BoxData[
\(\(ExtendedMcMillanForm[at_, xt_: {}] :=
Module[{a = at, x = xt, s, u, v, lcm, coef},
If[x === {} && Length[Variables[a]] > 1,
Message[General::badarg]; \[IndentingNewLine]Return[$Failed],
If[x === {},
x = Variabile[
a]]]; \[IndentingNewLine] (*Compute\ the\ monic\ lcm\ of\ \
the\ denominators*) lcm =
Apply[PolynomialLCM,
Flatten[Denominator[Together[a]]]]; \[IndentingNewLine]coef =
Coefficient[lcm, x, Exponent[lcm, x]]; \[IndentingNewLine]If[
coef =!= 1,
lcm = lcm/
coef]; \[IndentingNewLine] (*Compute\ the\ Smith\ form\ s\ \
of\ the\ polynomial\ matrix\ a*lcm*) {s, u, v} =
ExtendedSmithForm[Apart[a*lcm],
x]; \[IndentingNewLine] (*The\ McMillan\ form\ is\ s/
lcm*) {Together[s/lcm], u, v}] /;
MatrixQ[at] ||
Message[General::mtrx, "\<ExtendedMcMillanForm\>"];\)\)], "Input",\
InitializationCell->True],
Cell["Same as above without unimodular matrices", "SmallText"],
Cell[BoxData[
\(\(McMillanForm[at_, xt_: {}] :=
Module[{a = at, x = xt, s, lcm, coef},
If[x === {} && Length[Variables[a]] > 1,
Message[General::badarg]; \[IndentingNewLine]Return[$Failed],
If[x === {}, x = Variabile[a]]]; \[IndentingNewLine]lcm =
Apply[PolynomialLCM,
Flatten[Denominator[Together[a]]]]; \[IndentingNewLine]coef =
Coefficient[lcm, x, Exponent[lcm, x]]; \[IndentingNewLine]If[
coef =!= 1, lcm = lcm/coef]; \[IndentingNewLine]s =
SmithForm[Apart[a*lcm], x]; \[IndentingNewLine]Together[
s/lcm]] /;
MatrixQ[at] ||
Message[General::mtrx, "\<McMillanForm\>"];\)\)], "Input",
InitializationCell->True],
Cell["\<\
The following auxiliary function eliminates the dependent rows of \
the polynomial matrix at\
\>", "SmallText"],
Cell[BoxData[
\(\(EliminateDep[at_, xt_: {}] :=
Module[{a = \(b = at\), x = xt, u, p, m, rmin, k, kcl, ch = False},
If[x === {} && Length[Variables[a]] > 1,
Message[General::badarg]; \[IndentingNewLine]Return[$Failed],
If[x === {}, x = Variabile[a]]]; \[IndentingNewLine]If[
Not[ExtPolQ[a, x]], Message[General::pol];
Return[$Failed]]; \[IndentingNewLine]If[
Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[a]]],
a = a /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \
\[IndentingNewLine]ch = True]; \[IndentingNewLine]{p, m} =
Dimensions[a]; \[IndentingNewLine]u =
IdentityMatrix[p]; \[IndentingNewLine]For[k = \(kcl = 1\),
k \[LessEqual] p && kcl \[LessEqual] m, \(k++\),
While[a[\([Range[k, p], kcl]\)] === Table[0, {p - k + 1}] &&
kcl < m, \(kcl++\)]; \[IndentingNewLine]While[
If[k \[Equal] p, False,
a[\([Range[k + 1, p], kcl]\)] =!= Table[0, {p - k}]],
rmin = FindMin[\(Transpose[a]\)[\([kcl]\)], k, x,
p]; \[IndentingNewLine]If[rmin \[NotEqual] k,
a = a /. {a[\([rmin]\)] \[Rule] a[\([k]\)],
a[\([k]\)] \[Rule]
a[\([rmin]\)]}; \[IndentingNewLine]u =
u /. {u[\([rmin]\)] \[Rule] u[\([k]\)],
u[\([k]\)] \[Rule]
u[\([rmin]\)]}]; \[IndentingNewLine]For[i = k + 1,
i \[LessEqual] p, \(i++\),
If[a[\([i, kcl]\)] =!= 0,
lpc = PolynomialQuotient[a[\([i, kcl]\)],
a[\([k, kcl]\)], x]; \[IndentingNewLine]a =
ReplacePart[a, Expand[a[\([i]\)] - a[\([k]\)]*lpc, x],
i];]]]; (*Endwhile*) \(kcl++\)]; (*endfor*) k =
Flatten[Position[Inverse[u] . a,
n_ /; n === Table[0, {m}]]]; \[IndentingNewLine]If[ch,
b = b /. {x^\((n : _)\) \[Rule] x^\(-n\),
x \[Rule]
x^\(-1\)}]; \[IndentingNewLine]Return[\(IdentityMatrix[
p]\)[\([Complement[Range[1, p], k]]\)] . b]] /;
MatrixQ[at] ||
Message[General::mtrx, "\<EliminateDep\>"];\)\)], "Input",
InitializationCell->True],
Cell["\<\
The following auxiliary function computes a left GCD and a right \
lcm of matrices at and bt\
\>", "SmallText"],
Cell[BoxData[
\(\(AuxLGCDlcm[at_, bt_, xt_] :=
Module[{a = at, b = bt, x = xt, l, m, q, v, f, i, j, done, krow,
lambda, esp, coef, deg, ch = False},
If[Not[ExtPolQ[{a, b}, x]], Message[General::pol];
Return[$Failed]]; \[IndentingNewLine]If[
Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a,
b} = {a,
b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \
\[IndentingNewLine]ch =
True]; \[IndentingNewLine] (*Initialize\ the\ variables\
*) \[IndentingNewLine]{l, m} =
Dimensions[
a]; \[IndentingNewLine]q = \(Dimensions[
b]\)[\([2]\)]; \[IndentingNewLine]v =
IdentityMatrix[m + q]; \[IndentingNewLine]f =
Transpose[
Flatten[{Transpose[a], Transpose[b]},
1]]; \[IndentingNewLine] (*Main\ loop\ on\ k,
columns\ of\ f*) \[IndentingNewLine]For[k = \(krow = 1\),
krow \[NotEqual] l + 1 && k \[NotEqual] m + q, \(k++\),
done = False; \[IndentingNewLine]While[
f[\([krow, Range[k, m + q]]\)] === Table[0, {m + q - k + 1}] &&
krow < l, \(krow++\)]; \[IndentingNewLine] (*Eliminate\ all\ \
the\ elements\ of\ row\ krow, starting\ from\ column\ k + 1*)
While[f[\([krow, Range[k + 1, m + q]]\)] =!=
Table[0, {m + q - k}] && Not[done],
j = FindMin[f[\([krow]\)], k, x,
m + q]; \[IndentingNewLine]If[j \[NotEqual] k,
f = Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], \
#[\([j]\)] \[Rule] #[\([k]\)]} &\)[Transpose[f]]]; \[IndentingNewLine]v =
Transpose[\(# /. {#[\([k]\)] \[Rule] #[\([j]\)], \
#[\([j]\)] \[Rule] #[\([k]\)]} &\)[
Transpose[
v]]]]; \[IndentingNewLine] (*If\ on\ row\ krow\ all\
\ the\ elements\ following\ a[\([krow, k]\)]\ are\ zero\ we\ are\ done,
otherwise\ lower\ the\ degree\ of\ these\ non -
zero\ elements*) If[
Onenonzero[f[\([krow, Range[k, m + q]]\)], x], done = True,
For[i = k + 1, i \[LessEqual] m + q, \(i++\),
If[f[\([krow, i]\)] =!= 0,
lambda =
PolynomialQuotient[f[\([krow, i]\)],
f[\([krow, k]\)], x]; \[IndentingNewLine]f =
Transpose[\(ReplacePart[#,
Expand[#[\([i]\)] - #[\([k]\)] \((lambda)\),
x], i] &\)[
Transpose[f]]]; \[IndentingNewLine]v =
Transpose[\(ReplacePart[#,
Expand[#[\([i]\)] - #[\([k]\)] \((lambda)\),
x], i] &\)[
Transpose[
v]]]]] (*endfor*) ] (*endif*) ] (*endwhile*) \ \
\(krow++\)]; (*endfor*) (*Count\ the\ number\ of\ last\ zero\ columns\ in\ \
f*) j = q + m; \[IndentingNewLine]While[
j \[GreaterEqual] 1 && \(Transpose[f]\)[\([j]\)] ===
Table[0, {l}],
j = j - 1]; \[IndentingNewLine] (*Return\ the\ left\ GCD\ and\ if\
\ Flatten[lambda] =!= {}\ return\ also\ a\ right\ lcm*) If[
ch, {a, v,
f} = {a, v, f} /. {x^\((n : _)\) \[Rule] x^\(-n\),
x \[Rule] x^\(-1\)}]; \[IndentingNewLine]lambda =
v[\([Range[1, m], Range[j + 1, q + m]]\)]; \[IndentingNewLine]{If[
l > \((m + q)\),
Transpose[
Flatten[{Transpose[f[\([Range[1, l], Range[1, m + q]]\)]],
Table[0, {l - m - q}, {l}]}, 1]],
f[\([Range[1, l], Range[1, l]]\)]],
If[Flatten[lambda] =!= {},
If[l < m + q,
Transpose[EliminateDep[Transpose[Expand[a . lambda, x]], x]],
Expand[a . lambda, x], {}], {}]}];\)\)], "Input",
InitializationCell->True],
Cell[BoxData[
\(\(LGCD[at_, bt_, xt_: {}] :=
Module[{a = at, b = bt, x = xt},
If[x === {} && Length[Variables[{a, b}]] > 1,
Message[General::badarg]; \[IndentingNewLine]Return[$Failed],
If[x === {},
x = Variabile[{a,
b}]]]; \[IndentingNewLine]\(If[# =!= $Failed, \
#[\([1]\)], #] &\)[AuxLGCDlcm[a, b, x]]] /;
MatrixQ[at] &&
MatrixQ[bt] && \(Dimensions[at]\)[\([1]\)] \[Equal] \(Dimensions[
bt]\)[\([1]\)];\)\)], "Input",
InitializationCell->True],
Cell[BoxData[
\(\(Rlcm[at_, bt_, xt_: {}] :=
Module[{a = at, b = bt, x = xt},
If[x === {} && Length[Variables[{a, b}]] > 1,
Message[General::badarg]; \[IndentingNewLine]Return[$Failed],
If[x === {},
x = Variabile[{a,
b}]]]; \[IndentingNewLine]\(If[# =!= $Failed, \
#[\([2]\)], #] &\)[AuxLGCDlcm[a, b, x]]] /;
MatrixQ[at] &&
MatrixQ[bt] && \(Dimensions[at]\)[\([1]\)] \[Equal] \(Dimensions[
bt]\)[\([1]\)];\)\)], "Input",
InitializationCell->True],
Cell["\<\
The following auxiliary function computes a right GCD and a left \
lcm of matrices at and bt\
\>", "SmallText"],
Cell[BoxData[
\(\(AuxRGCDlcm[at_, bt_, xt_] :=
Module[{a = at, b = bt, x = xt, l, m, p, v, f, i, j, done, kcl,
lambda, esp, coef, deg, ch = False},
If[Not[ExtPolQ[{a, b}, x]], Message[General::pol];
Return[$Failed]]; \[IndentingNewLine]If[
Not[Apply[And, \(PolynomialQ[#, x] &\) /@ Flatten[{a, b}]]], {a,
b} = {a,
b} /. {\((x^\((n : _)\))\) \[Rule] \((x^\(-n\))\)}; \
\[IndentingNewLine]ch =
True]; \[IndentingNewLine] (*Initialize\ the\ variables*) {l,
m} = Dimensions[
a]; \[IndentingNewLine]p = \(Dimensions[
b]\)[\([1]\)]; \[IndentingNewLine]v =
IdentityMatrix[l + p]; \[IndentingNewLine]f =
Flatten[{a, b}, 1]; \[IndentingNewLine] (*Main\ loop\ on\ k,
rows\ of\ f*) For[k = \(kcl = 1\),
kcl \[NotEqual] m + 1 && k \[NotEqual] l + p, \(k++\),
done = False; \[IndentingNewLine]While[
f[\([Range[k, l + p], kcl]\)] === Table[0, {l + p - k + 1}] &&
kcl < m, \(kcl++\)]; \[IndentingNewLine] (*Eliminate\ all\ \
the\ elements\ of\ column\ kcl, starting\ from\ row\ k + 1*)
While[f[\([Range[k + 1, l + p], kcl]\)] =!=
Table[0, {l + p - k}] && Not[done],
j = FindMin[\(Transpose[f]\)[\([kcl]\)], k, x,
l + p]; \[IndentingNewLine]If[j \[NotEqual] k,
f = f /. {f[\([k]\)] \[Rule] f[\([j]\)],
f[\([j]\)] \[Rule]
f[\([k]\)]}; \[IndentingNewLine]v =
v /. {v[\([k]\)] \[Rule] v[\([j]\)],
v[\([j]\)] \[Rule]
v[\([k]\)]}]; \[IndentingNewLine] (*If\ on\ column\ \
kcl\ all\ the\ elements\ below\ a[\([k, kcl]\)]\ are\ zero\ we\ are\ done,
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -