ComplexNonlinearEquationVerifier.java

  1. /*
  2.  * Created on 2004/12/02
  3.  * Copyright (C) 2004 Koga Laboratory. All rights reserved.
  4.  *
  5.  */
  6. package org.mklab.cga.nonlinear;

  8. import java.util.List;

  10. import org.mklab.cga.derivative.IntervalDerivativeComplexMatrix;
  11. import org.mklab.cga.derivative.IntervalDerivativeComplexNumber;
  12. import org.mklab.cga.derivative.IntervalDerivativeFunction;
  13. import org.mklab.cga.derivative.IntervalDerivativeFunctionEvaluation;
  14. import org.mklab.cga.derivative.IntervalDerivativeRealMatrix;
  15. import org.mklab.cga.derivative.IntervalDerivativeRealNumber;
  16. import org.mklab.cga.interval.matrix.IntervalComplexNumericalMatrix;
  17. import org.mklab.cga.interval.matrix.IntervalRealNumericalMatrix;
  18. import org.mklab.cga.interval.scalar.IntervalComplexNumericalScalar;
  19. import org.mklab.cga.interval.scalar.IntervalRealNumericalScalar;
  20. import org.mklab.cga.util.IntervalUtil;
  21. import org.mklab.nfc.matrix.ComplexNumericalMatrix;
  22. import org.mklab.nfc.matrix.NormType;
  23. import org.mklab.nfc.matrix.RealNumericalMatrix;
  24. import org.mklab.nfc.scalar.ComplexNumericalScalar;
  25. import org.mklab.nfc.scalar.RealNumericalScalar;


  28. /**
  29.  * 非線形方程式の精度保証を行うクラスです。
  30.  * 
  31.  * @author hiroki
  32.  * @version $Revision: 1.24 $.2004/12/02
  33.  * @param <RIDS> 実区間微分スカラーの型
  34.  * @param <RIDM> 実区間微分行列の型
  35.  * @param <CIDS> 複素区間微分スカラーの型
  36.  * @param <CIDM> 複素区間微分行列の型
  37.  * @param <RIS> 実区間スカラーの型
  38.  * @param <RIM> 実区間行列の型
  39.  * @param <CIS> 複素区間スカラーの型
  40.  * @param <CIM> 複素区間行列の型
  41.  * @param <RS> 実スカラーの型
  42.  * @param <RM> 実行列の型
  43.  * @param <CS> 複素スカラーの型
  44.  * @param <CM> 複素行列の型
  45.  */
  46. public class ComplexNonlinearEquationVerifier<RIDS extends IntervalDerivativeRealNumber<RIDS,RIDM,CIDS,CIDM,RIS,RIM,CIS,CIM,RS,RM,CS,CM>, RIDM extends IntervalDerivativeRealMatrix<RIDS,RIDM,CIDS,CIDM,RIS,RIM,CIS,CIM,RS,RM,CS,CM>, CIDS extends IntervalDerivativeComplexNumber<RIDS,RIDM,CIDS,CIDM,RIS,RIM,CIS,CIM,RS,RM,CS,CM>, CIDM extends IntervalDerivativeComplexMatrix<RIDS,RIDM,CIDS,CIDM,RIS,RIM,CIS,CIM,RS,RM,CS,CM>, RIS extends IntervalRealNumericalScalar<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, RIM extends IntervalRealNumericalMatrix<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, CIS extends IntervalComplexNumericalScalar<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, CIM extends IntervalComplexNumericalMatrix<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, RS extends RealNumericalScalar<RS,RM,CS,CM>, RM extends RealNumericalMatrix<RS,RM,CS,CM>, CS extends ComplexNumericalScalar<RS,RM,CS,CM>, CM extends ComplexNumericalMatrix<RS,RM,CS,CM>>{

  48.   /** 関数オブジェクト */
  49.   private IntervalDerivativeFunction<CIDS,CIDM,CIS,CIM,CS,CM> function;

  51. //  /** 初期値 */
  52. //  private M initialValue;
  53. //
  54. //  /** 初期値区間 */
  55. //  private IM intervalInitialValue;

  57.   /** 非線形方程式の解 */
  58.   private CIM solution;

  60.   /**
  61.    * コンストラクタ
  62.    * 
  63.    * @param function 関数オブジェクト
  64.    */
  65.   public ComplexNonlinearEquationVerifier(IntervalDerivativeFunction<CIDS,CIDM,CIS,CIM,CS,CM> function) {
  66.     this.function = function;
  67.   }

  69.   /**
  70.    * 非線型方程式の解を精度保証付きで求めます。
  71.    * 
  72.    * @param xs 初期値
  73.    * @param idm idm
  74.    * @return 非線形方程式の解
  75.    */
  76.   public CIM solve(CM xs, CIDM idm) {
  77.     CIDS ids = idm.getElement(1,1);
  78.     CIS is = ids.getX();
  79.     CIM im = idm.getX();
  80.     
  81.     int n = xs.getRowSize();
  82.     double sqrtN = Math.sqrt(n);

  84.     CM xss = xs.createClone();
  85.     CM xsold = xss;

  87.     CS element = xs.getElement(1);

  89.     CM dxs = xs.createZero();
  90.     CM dxsnew = (xs).absElementWise();

  92.     int k = 0;

  94.     CM y_x;
  95.     CM y_dx = null;

  97.     CIDS[] idxss = ids.createArray(xs.getRowSize());
  98.     IntervalDerivativeFunctionEvaluation<CIDS,CIDM,CIS,CIM,CS,CM> nee = new IntervalDerivativeFunctionEvaluation<>(this.function);

  100.     List<CIM> value1;

  102.     while (dxsnew.compareElementWise(".<", dxs.multiply(5)).anyTrue() //$NON-NLS-1$
  103.         && (dxsnew).norm(NormType.TWO).isGreaterThan((xs).norm(NormType.TWO).multiply(sqrtN * 1e-14)) && k < 100 || k < 3) {
  104.       k = k + 1;
  105.       dxs = dxsnew.createClone();
  106.       xsold = xss.createClone();

  108.       for (int i = 0; i < xss.getRowSize(); i++) {
  109.         idxss[i] = ids.create(is.create((xss).getElement(i + 1)), is.create(element.createUnit()));
  110.       }
  111.       value1 = nee.getValues(idxss);
  112.       // f(x)
  113.       y_x = value1.get(0).getMiddle();
  114.       // ヤコビ行列
  115.       y_dx = value1.get(2).getMiddle();
  116.       xss = xss.subtract(y_dx.leftDivide(y_x));
  117.       dxsnew = xss.subtract(xsold).absElementWise();
  118.     }

  120.     CIM ixss = im.create(xss);

  122.     for (int i = 0; i < xss.getRowSize(); i++) {
  123.       idxss[i] = ids.create(ixss.getElement(i + 1),  is.create(element.createUnit()));
  124.     }

  126.     // f(xss)とdf(xss)とヤコビ行列を求める
  127.     List<CIM> values = nee.getValues(idxss);
  128.     // apporoximate inver of IntervalDerivativeFunctionEvaluation
  129.     // TODO 区間にしてから、逆行列？
  130.     CIM R = im.create(y_dx.inverse());
  131.     // initialization of interval iteration
  132.     // TODO 多倍長だと？
  133.     CIM Z = R.unaryMinus().multiply(values.get(0));
  134.     CIM X = Z.createClone();

  136.     boolean ready = false;
  137.     k = 0;
  138.     CIM Y = null;

  140.     while (ready == false && k < 10) {
  141.       CM xRadius = X.getSupremum().subtract(X.getInfimum()).divide(2);
  142.       
  143.       CIM E = im.create(xRadius.multiply(element.createUnit().divide(element.createUnit().multiply(10)))).multiply(is.createInfSup(element.createUnit().unaryMinus(), element.createUnit()))
  144.             .addElementWise(is.createInfSup(element.getMachineEpsilon().unaryMinus(), element.getMachineEpsilon()));

  146.       k = k + 1;
  147.       // TODO 本当はhull(X+E,0)?
  148.       Y = X.add(E);
  149.       CIM Yold = Y.createClone();

  151.       for (int i = 0; i < xs.getColumnSize(); i++) {
  152.         idxss[i] = ids.create(Y.add(im.create(xss)).getElement(i + 1), is.create(element.createUnit()));
  153.       }

  155.       values = nee.getValues(idxss);
  156.       // f(xs+Y)
  157.       // ? IntervalMatrix Iy_x = values[0];
  158.       // ヤコビ行列
  159.       CIM Iy_dx = values.get(2);
  160.       CIM C = im.create(xs.createUnit(n)).subtract(R.multiply(Iy_dx));

  162.       int i = 0;
  163.       // improved interval iteration
  164.       while (!ready && i < 2) {
  165.         i = i + 1;
  166.         X = Z.add(C.multiply(Y));
  167.         ready = IntervalUtil.in0forRump(X, Y).allTrue();
  168.         Y = IntervalUtil.intersect(X, Yold);
  169.       }
  170.     }

  172.     if (ready) {
  173.       X = Y.add(im.create(xss));
  174.     } else {
  175.       throw new IllegalArgumentException();
  176.       // TODO 例外投げるか、NaNを返すようにする
  177.     }

  179.     return X;
  180.   }

  182.   /**
  183.    * 非線型方程式の解を精度保証付きで求めます。
  184.    * 
  185.    * <p>初期値区間に真の解が含まれる場合、非線形方程式の真の解を含むより小さな区間を返します。
  186.    * 
  187.    * @param xs 初期値区間
  188.    * @param idm idm
  189.    * @return 非線形方程式の解
  190.    */
  191.   public CIM solve(CIM xs, CIDM idm) {
  192.     CIDS ids = idm.getElement(1,1);
  193.     CIS is = ids.getX();
  194.     CIM im = idm.getX();
  195.     
  196.     IntervalDerivativeFunctionEvaluation<CIDS,CIDM,CIS,CIM,CS,CM> nee = new IntervalDerivativeFunctionEvaluation<>(this.function);
  197.     CIM X = xs.createClone();

  199.     CS element = (xs.getMiddle()).getElement(1);
  200.     
  201.     boolean ready = false;
  202.     int k = 0;
  203.     // TODO kの値はどれくらい？
  204.     while (ready == false && k < 10) {
  205.       k = k + 1;

  207.       CIDS[] idxs = ids.createArray(X.getRowSize());
  208.       for (int i = 0; i < X.getRowSize(); i++) {
  209.         idxs[i] = ids.create(is.create(X.getMiddle().getElement(i + 1)), is.create(element.createUnit()));
  210.       }

  212.       // f(xs)とdf(xs)とヤコビ行列を求める
  213.       List<CIM> values = nee.getValues(idxs);
  214.       CM y_dx = values.get(2).getMiddle();
  215.       // apporoximate inver of IntervalDerivativeFunctionEvaluation
  216.       CIM R = im.create(y_dx.inverse());
  217.       // initialization of interval iteration
  218.       CIM Z = R.unaryMinus().multiply(values.get(0));
  219.       CIM Y = X.subtract(im.create(X.getMiddle()));

  221.       CIDS[] idX = ids.createArray(X.getRowSize());
  222.       for (int i = 0; i < X.getRowSize(); i++) {
  223.         idX[i] = ids. create(X.getElement(i + 1),  is.create(element.createUnit()));
  224.       }

  226.       values = nee.getValues(idX);
  227.       // ヤコビ行列 f'(X)
  228.       CIM Iy_dx = values.get(2);
  229.       int n = X.getRowSize();
  230.       CIM C = im.create(xs.getMiddle().createUnit(n)).subtract(R.multiply(Iy_dx));

  232.       CIM nextX = Z.add(C.multiply(Y)).add(im.create(X.getMiddle()));
  233.       ready = IntervalUtil.in0forRump(nextX, X).allTrue();

  235.       // 不動点定理によれば、元の集合を返すべき      
  236.       //      if (ready) {
  237.       //        X = (IntervalMatrix)nextX.clone();
  238.       //      }

  240.     }

  242.     //    if (ready == false) {
  243.     //      throw new IllegalArgumentException();
  244.     //    }

  246.     return X;
  247.   }

  249. //  /**
  250. //   * 
  251. //   */
  252. //  public void solve() {
  253. //    if (this.initialValue != null) {
  254. //      this.solution = this.solve(this.initialValue);
  255. //    } else {
  256. //      this.solution = this.solve(this.intervalInitialValue);
  257. //    }
  258. //  }

  260. //  /**
  261. //   * 初期値を設定します。
  262. //   * 
  263. //   * @param initialValue 初期値
  264. //   */
  265. //  public void setInitialValue(M initialValue) {
  266. //    this.initialValue = initialValue;
  267. //  }
  268. //
  269. //  /**
  270. //   * 初期値区間を設定します。
  271. //   * 
  272. //   * @param intervalInitialValue 初期値区間
  273. //   */
  274. //  public void setIntervalInitialValue(IM intervalInitialValue) {
  275. //    this.intervalInitialValue = intervalInitialValue;
  276. //  }
  277. }
Created with JaCoCo 0.8.3.201901230119