RealNonlinearEquationVerifier.java
/*
* Created on 2004/12/02
* Copyright (C) 2004 Koga Laboratory. All rights reserved.
*
*/
package org.mklab.cga.nonlinear;
import java.util.List;
import org.mklab.cga.derivative.IntervalDerivativeComplexMatrix;
import org.mklab.cga.derivative.IntervalDerivativeComplexNumber;
import org.mklab.cga.derivative.IntervalDerivativeFunction;
import org.mklab.cga.derivative.IntervalDerivativeFunctionEvaluation;
import org.mklab.cga.derivative.IntervalDerivativeRealMatrix;
import org.mklab.cga.derivative.IntervalDerivativeRealNumber;
import org.mklab.cga.interval.matrix.IntervalComplexNumericalMatrix;
import org.mklab.cga.interval.matrix.IntervalRealNumericalMatrix;
import org.mklab.cga.interval.scalar.IntervalComplexNumericalScalar;
import org.mklab.cga.interval.scalar.IntervalRealNumericalScalar;
import org.mklab.cga.util.IntervalUtil;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
/**
* 非線形方程式の精度保証を行うクラスです。
*
* @author hiroki
* @version $Revision: 1.24 $.2004/12/02
* @param <RIDS> 実区間微分スカラーの型
* @param <RIDM> 実区間微分行列の型
* @param <CIDS> 複素区間微分スカラーの型
* @param <CIDM> 複素区間微分行列の型
* @param <RIS> 実区間スカラーの型
* @param <RIM> 実区間行列の型
* @param <CIS> 複素区間スカラーの型
* @param <CIM> 複素区間行列の型
* @param <RS> 実スカラーの型
* @param <RM> 実行列の型
* @param <CS> 複素スカラーの型
* @param <CM> 複素行列の型
*/
public class RealNonlinearEquationVerifier<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>> {
/** 関数オブジェクト */
private IntervalDerivativeFunction<RIDS,RIDM,RIS,RIM,RS,RM> function;
// /** 初期値 */
// private M initialValue;
//
// /** 初期値区間 */
// private IM intervalInitialValue;
/** 非線形方程式の解 */
private RIM solution;
/**
* コンストラクタ
*
* @param function 関数オブジェクト
*/
public RealNonlinearEquationVerifier(IntervalDerivativeFunction<RIDS,RIDM,RIS,RIM,RS,RM> function) {
this.function = function;
}
/**
* 非線型方程式の解を精度保証付きで求めます。
*
* @param xs 初期値
* @param idm idm
* @return 非線形方程式の解
*/
public RIM solve(RM xs, RIDM idm) {
RIDS ids = idm.getElement(1,1);
RIS is = ids.getX();
RIM im = idm.getX();
int n = xs.getRowSize();
double sqrtN = Math.sqrt(n);
RM xss = xs.createClone();
RM xsold = xss;
RS element = xs.getElement(1);
RM dxs = xs.createZero();
RM dxsnew = (xs).absElementWise();
int k = 0;
RM y_x;
RM y_dx = null;
RIDS[] idxss = ids.createArray(xs.getRowSize());
IntervalDerivativeFunctionEvaluation<RIDS,RIDM,RIS,RIM,RS,RM> nee = new IntervalDerivativeFunctionEvaluation<>(this.function);
List<RIM> value1;
while (dxsnew.compareElementWise(".<", dxs.multiply(5)).anyTrue() //$NON-NLS-1$
&& (dxsnew).norm(NormType.TWO).isGreaterThan((xs).norm(NormType.TWO).multiply(sqrtN * 1e-14)) && k < 100 || k < 3) {
k = k + 1;
dxs = dxsnew.createClone();
xsold = xss.createClone();
for (int i = 0; i < xss.getRowSize(); i++) {
idxss[i] = ids.create(is.create((xss).getElement(i + 1)), is.create(element.createUnit()));
}
value1 = nee.getValues(idxss);
// f(x)
y_x = value1.get(0).getMiddle();
// ヤコビ行列
y_dx = value1.get(2).getMiddle();
xss = xss.subtract(y_dx.leftDivide(y_x));
dxsnew = xss.subtract(xsold).absElementWise();
}
RIM ixss = im.create(xss);
for (int i = 0; i < xss.getRowSize(); i++) {
idxss[i] = ids.create(ixss.getElement(i + 1), is.create(element.createUnit()));
}
// f(xss)とdf(xss)とヤコビ行列を求める
List<RIM> values = nee.getValues(idxss);
// apporoximate inver of IntervalDerivativeFunctionEvaluation
// TODO 区間にしてから、逆行列？
RIM R = im.create(y_dx.inverse());
// initialization of interval iteration
// TODO 多倍長だと？
RIM Z = R.unaryMinus().multiply(values.get(0));
RIM X = Z.createClone();
boolean ready = false;
k = 0;
RIM Y = null;
while (ready == false && k < 10) {
RM xRadius = X.getSupremum().subtract(X.getInfimum()).divide(2);
RIM E = im.create(xRadius.multiply(element.createUnit().divide(element.createUnit().multiply(10)))).multiply(is.createInfSup(element.createUnit().unaryMinus(), element.createUnit()))
.addElementWise(is.createInfSup(element.getMachineEpsilon().unaryMinus(), element.getMachineEpsilon()));
k = k + 1;
// TODO 本当はhull(X+E,0)?
Y = X.add(E);
RIM Yold = Y.createClone();
for (int i = 0; i < xs.getColumnSize(); i++) {
idxss[i] = ids.create(Y.add(im.create(xss)).getElement(i + 1), is.create(element.createUnit()));
}
values = nee.getValues(idxss);
// f(xs+Y)
// ? IntervalMatrix Iy_x = values[0];
// ヤコビ行列
RIM Iy_dx = values.get(2);
RIM C = im.create(xs.createUnit(n)).subtract(R.multiply(Iy_dx));
int i = 0;
// improved interval iteration
while (!ready && i < 2) {
i = i + 1;
X = Z.add(C.multiply(Y));
ready = IntervalUtil.in0forRump(X, Y).allTrue();
Y = IntervalUtil.intersect(X, Yold);
}
}
if (ready) {
X = Y.add(im.create(xss));
} else {
throw new IllegalArgumentException();
// TODO 例外投げるか、NaNを返すようにする
}
return X;
}
/**
* 非線型方程式の解を精度保証付きで求めます。
*
* <p>初期値区間に真の解が含まれる場合、非線形方程式の真の解を含むより小さな区間を返します。
*
* @param xs 初期値区間
* @param idm idm
* @return 非線形方程式の解
*/
public RIM solve(RIM xs, RIDM idm) {
RIDS ids = idm.getElement(1,1);
RIS is = ids.getX();
RIM im = idm.getX();
IntervalDerivativeFunctionEvaluation<RIDS,RIDM,RIS,RIM,RS,RM> nee = new IntervalDerivativeFunctionEvaluation<>(this.function);
RIM X = xs.createClone();
RS element = (xs.getMiddle()).getElement(1);
boolean ready = false;
int k = 0;
// TODO kの値はどれくらい？
while (ready == false && k < 10) {
k = k + 1;
RIDS[] idxs = ids.createArray(X.getRowSize());
for (int i = 0; i < X.getRowSize(); i++) {
idxs[i] = ids.create(is.create(X.getMiddle().getElement(i + 1)), is.create(element.createUnit()));
}
// f(xs)とdf(xs)とヤコビ行列を求める
List<RIM> values = nee.getValues(idxs);
RM y_dx = values.get(2).getMiddle();
// apporoximate inver of IntervalDerivativeFunctionEvaluation
RIM R = im.create(y_dx.inverse());
// initialization of interval iteration
RIM Z = R.unaryMinus().multiply(values.get(0));
RIM Y = X.subtract(im.create(X.getMiddle()));
RIDS[] idX = ids.createArray(X.getRowSize());
for (int i = 0; i < X.getRowSize(); i++) {
idX[i] = ids. create(X.getElement(i + 1), is.create(element.createUnit()));
}
values = nee.getValues(idX);
// ヤコビ行列 f'(X)
RIM Iy_dx = values.get(2);
int n = X.getRowSize();
RIM C = im.create(xs.getMiddle().createUnit(n)).subtract(R.multiply(Iy_dx));
RIM nextX = Z.add(C.multiply(Y)).add(im.create(X.getMiddle()));
ready = IntervalUtil.in0forRump(nextX, X).allTrue();
// 不動点定理によれば、元の集合を返すべき
// if (ready) {
// X = (IntervalMatrix)nextX.clone();
// }
}
// if (ready == false) {
// throw new IllegalArgumentException();
// }
return X;
}
// /**
// *
// */
// public void solve() {
// if (this.initialValue != null) {
// this.solution = this.solve(this.initialValue);
// } else {
// this.solution = this.solve(this.intervalInitialValue);
// }
// }
// /**
// * 初期値を設定します。
// *
// * @param initialValue 初期値
// */
// public void setInitialValue(M initialValue) {
// this.initialValue = initialValue;
// }
//
// /**
// * 初期値区間を設定します。
// *
// * @param intervalInitialValue 初期値区間
// */
// public void setIntervalInitialValue(IM intervalInitialValue) {
// this.intervalInitialValue = intervalInitialValue;
// }
}