Faddeev.java

/*
 * $Id: Faddeev.java,v 1.23 2008/03/23 14:29:08 koga Exp $
 *
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.control;

import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.scalar.DoublePolynomial;


/**
 * ファデーブのアルゴリズムを用いて特性多項式を求めるクラスです。
 * 
 * <p>Characteristic polynomial by Faddeev's algorithm
 * 
 * @author koga
 * @version $Revision: 1.23 $
 * @see org.mklab.tool.matrix.Makepoly
 * @see org.mklab.tool.control.Resolvent
 */
public class Faddeev {

  /**
   * メインメソッド
   * 
   * @param args コマンドライン引数
   */
  public static void main(String[] args) {
    DoubleMatrix A = new DoubleMatrix(new double[][] { {0, 1, 0}, {0, 0, 1}, {-2, -3, -4}});
    FaddeevSolution tmp = faddeev(A);
    DoubleMatrix NN = tmp.getNN();
    DoublePolynomial dd = tmp.getDd();
    NN.print("NN"); //$NON-NLS-1$
    dd.print("dd"); //$NON-NLS-1$
  }

  /**
   * 行列 <code>A</code> のリゾルベント行列
   * 
   * <pre><code>
   * 
   * NN(s) NN_(n-1) s&circ;(n-1) + ... + NN1 s + NN0 -1 ----- = ------------------------------------------- = (s*I - A) dd(s) s&circ;n + dd_(n-1) s&circ;(n-1) + ... + dd1 s + dd0
   * 
   * </code></pre>
   * 
   * の分子行列多項式の係数行列と分母多項式の係数を返します。
   * 
   * <code>NN</code> は、分子多項式行列の係数行列を
   * 
   * <pre><code>
   * 
   * NN = [NN_0, NN_1, ..., NN_(n-1)]
   * 
   * </code></pre>
   * 
   * のように含む。
   * 
   * @param A A行列
   * @return {NN, dd} (伝達関数行列の分子多項式行列, 伝達関数の分母多項式) characteristic polynomial
   */
  public static FaddeevSolution faddeev(DoubleMatrix A) {
    int n = A.getRowSize();
    DoubleMatrix AA = A.createClone();
    DoubleMatrix NN = AA.createZero(AA.getRowSize(), n * AA.getColumnSize());
    DoublePolynomial dd = new DoublePolynomial(AA.createZero(1, n + 1));

    dd.setCoefficient(n, 1);

    NN.setSubMatrix(1, n, AA, A.createUnit(n));

    dd.setCoefficient(n - 1, AA.trace().unaryMinus());

    for (int i = n - 2; i >= 0; i--) {
      NN.setSubMatrix(1, i + 1, AA, AA.multiply(NN.getSubMatrix(1, i + 2, AA)).add(A.createUnit(n).multiply(dd.getCoefficient(i + 1))));
      dd.setCoefficient(i, AA.multiply(NN.getSubMatrix(1, i + 1, AA)).trace().unaryMinus().divide(n - i));
    }

    return new FaddeevSolution(NN,dd);
  }
}