連続系リカッティ方程式

• 1行と11行: 関数の始まり。
• 3行と14行: 局所変数の宣言。
• 6行: キーボードから行列を読み込む。
• 7行: 関数CRiccati()を呼ぶ。
• 8行: 解行列Pを表示する。
• 12行: 関数の引数の宣言。
• 16行: ハミルトン行列を作る[6]。
• 19行: 固有ベクトルの計算。
• 22行と23行: 行列vecから安定な固有値に対する固有ベクトルを取り出す。 vec(1:n,n+1:2*n)はvecの1行からn行，n+1列から2n列の部分行列， vec(n+1:2*n,n+1:2*n)はvecのn+1行から2n行，n+1列から2n列の部分 行列である。
• 25行と26行: 関数isreal()は与えられた行列が実行列なら1を， そうでなければ0を返す。
• 27行: 関数Re()は行列の実部を返す。

   Func void main()
   {
     Matrix A, B, Q, R, P;     // 局所変数の宣言
     Matrix CRiccati();        // 関数の宣言
05:
     read A, B, Q, R;          // 行列エディタを起動
     P = CRiccati(A, B, Q, R); // CRiccati()を呼ぶ
     print P;                  // Pを表示
   }
10:
   Func Matrix CRiccati(A, B, Q, R)
     Matrix A, B, Q, R;
   {
     Integer n;
15:  Matrix H, U, V, vec, P;

     n = Rows(A);
     H = [[ A  , -B*R~*B#]      // ハミルトン行列
          [-Q  , -A#     ]];
20:   
     vec = eigvec(H);           // H の固有ベクトル
     U = vec(    1:n, n+1:2*n); // 安定な固有値に対応
     V = vec(n+1:2*n, n+1:2*n); // する固有ベクトルを
                                // 取り出す
25:  if (isreal(A) && isreal(B) &&
         isreal(Q) && isreal(R)) {
         P = Re(V*U~);
     } else {
         P = V*U~;
30:  }
   
     return P;
   }