線形システムの伝達関数行列

   Func void main()
   {
     Matrix A, B, C, G;    // 変数の宣言
     Matrix Faddeev();
05:
     read A, B, C;         // 行列の読み込み
     G = Faddeev(A, B, C); // Faddeev()を呼ぶ
     print G;              // G を表示
   }
10:
   Func Matrix Faddeev(A, B, C)
     Matrix A, B, C;
   {
     Matrix Gamma, N;
15:  Polynomial s, ch;
     Integer i, n;
   
     s = $;            // 多項式の変数
     n = Rows(A);      // システムの次数
20:
     Gamma = Z(n,A);
     ch = Polynomial(Z(1,n+1)); 
                       // 特性多項式
     ch(n) = 1.0;      // 最高次の係数は 1
25:
     Gamma(n-1,A) = I(n);
     ch(n-1) = - trace(A);
   
     for (i = n-2; i >= 0; i--) {
30:      Gamma(i,A) = A*Gamma(i+1,A) + ch(i+1)*I(n);
         ch(i) = - trace(A*Gamma(i,A))/(n - i);
     }
   
     N = Z(Rows(C), Cols(B));   
35:  for (i = 0; i < n; i++) {
         N = N + C*Gamma(i,A)*B * s^i;
     }
   
     return N/ch;      // 伝達関数を返す
40:}