AL-LIN3.C |
/* * AL-LIN3.cpp * g++ -I/usr/local/include/profil -o AL-LIN3 AL-LIN3.cpp -lProfil -lBias -llr */ /* Ax=b の A,近似逆行列R,右辺b から精度保証つきの解の検証 (中尾・山本 p.14)。 教科書中の点行列、ベクトルの B,Z,x を区間で考えた (IB,IZ,xt)。 */ #include <iostream> #include <iomanip> /* setprecision */ #include "Constants.h" /* Machine::Epsilon */ #include "IntervalMatrix.h" /* 区間行列 */ using namespace std; //#define PRINT #define ZPRINT #define LOOP 300 #define DELTA Machine::Epsilon int main(void) { int i, j, N, k; INTERVAL delta(1.0 - DELTA, 1.0 + DELTA); cout << setprecision(16) ; cout << "Machine ε= " << Machine::Epsilon << endl; cout << " δ= " << DELTA << endl; cout << "何次元ですか ?" << endl; cin >> N; VECTOR x(N), b(N), Z(N); MATRIX A(N,N), R(N,N), I(N,N); INTERVAL_VECTOR xt(N), IZ(N), Y(N), X(N); INTERVAL_MATRIX IR(N,N), IA(N,N), IB(N,N); cout << "行列A" << endl; cin >> A; cout << A << endl; /* IA = Hull(A) */ IA = A; cout << "近似逆行列R" << endl; cin >> R; cout << R << endl; /* IR = Hull(R); */ IR = R; cout << "右辺b" << endl; cin >> b; cout << b << endl; /* 単位行列 */ Clear(I); for (i = 1; i <= N; i++) I(i,i) = 1.0; IB = I - IR * IA; xt = IR * b; IZ = IR * (b - IA * xt); X = Z; #ifdef PRINT cout << "区間行列 B" << endl << IB << endl; cout << "区間ベクトル ~x" << endl << xt << endl; cout << "区間行列 Z" << endl << IZ << endl; cout << "区間ベクトル X" << endl << X << endl; cout << "δ" << endl << delta << endl; #endif for (k = 1; k <= LOOP; k++) { /* X のデルタ拡大 */ Y = delta * X; X = IZ + IB * Y; #ifdef ZPRINT cout << k << "回目" << endl ; cout << "X=" << X << endl; cout << "Y= " << Y << endl; #endif /* (Inf(Y) < Inf(X) && Sup(X) < Sup(Y)) としていた */ if (X < Y) { cout << endl << "反復終了" << endl; cout << xt << endl << "+" << endl << X << endl << "=" << endl << xt+X << endl << "の中に解がある。(反復 " << k << "回)" << endl << "Xの幅は" << endl << Diam(X) << endl; exit(0); } } cout << endl << "反復回数オーバー" << endl << xt << endl << "+" << endl << X << endl << "=" << endl << xt+X << endl << "の中には解はありません。" << endl; return 0; } |