C..8 中尾・山本 第3章

chapter3ex2.cpp

/*
 * chapter3ex2.cpp
 *
 * g++ -I/usr/local/include/profil -o chapter3ex2 chapter3ex2.cpp -lProfil -lBias -llr
 *
 *  -u" - πu = (π - 1)sinπx を精度保証付きで求める。
 *     (中尾先生の本を少し改変)
 *
 * 近似解もこのプログラム内で求めている。
 * 
 */

#include <iostream>
#include <iomanip>
#include "IntervalMatrix.h"
#include "Constants.h"
#include "Functions.h"
#include "Utilities.h"
#include "LSS.h"

using namespace std;

#define PI Constant::Pi
#define DELTA 10e-10
#define LOOP 300

int main(void)
{

  INT N, done = 0;
  INT i, j, k;
  REAL h;
  INTERVAL Ih;

  cout << setprecision(16);
  
  cout << "Input N" << endl;
  cin >> N;
  
  h = 1.0 / (N + 1);
  Ih = 1.0 / Hull(N + 1);

  INTERVAL I, IPi;
  IPi = 2.0 * ArcSin(Hull(1));

  INTERVAL_MATRIX L(N, N), M(N, N);
  MATRIX D(N, N), Dinv(N, N), AbsDInv(N, N);
  INTERVAL_VECTOR f(N), u(N);
  VECTOR fm(N), x(N+2);

  for (i = 1; i <= N+2; i++)
    x(i) = (i - 1) * h;

  VECTOR Exact(N);
  for (i = 1; i <= N; i++)
    Exact(i) = 1.0 / PI * Sin(PI * h * i);

  Clear(L);  Clear(M);  Clear(D);  Clear(fm);

  D(1, 1) = 2 / h; D(1, 2) = - 1 / h;
  for (i = 2; i < N; i++){
    D(i, i) = 2 / h;
    D(i, i-1) = D(i, i+1) = - 1 / h;
  }
  D(N, N-1) = - 1 / h;
  D(N, N) = 2 / h;

  L(1, 1) = 2 * Ih / 3; L(1, 2) = Ih / 6;
  for (i = 2; i < N; i++){
    L(i, i) = 2 * Ih / 3;
    L(i, i-1) = L(i, i+1) = Ih / 6;
  }
  L(N, N-1) = Ih / 6;
  L(N, N) = 2 * Ih / 3;

  /* MATRIX Inverse (const MATRIX &); */
  Dinv = Inverse(D);
  AbsDInv = Abs(Dinv);

  M = D - PI * L;

  for (i = 1; i <= N; i++)
    fm(i) = (PI - 1) / (h * PI * PI) * 
      (2.0 * Sin(PI * x(i+1)) - (Sin(PI * x(i+2)) + Sin(PI * x(i))));

  f = Hull(fm);

  u = ILSS(M, f, done); 

  INT count;
  INTERVAL Iphi, sin_norm, Ialp, F_ort, sup_norm;
  INTERVAL_VECTOR Phi(N), C(N), Y(N), W(N);
  INTERVAL_VECTOR U(N);
  REAL SUP_norm;

  REAL width, WIDTH;

  Iphi = Sqrt(Ih * 2.0 / 3.0);
  Initialize(Phi, Iphi);
  /* ||sin πx||=1/√2 */
  sin_norm = 1 / Sqrt(Hull(2));

  Clear(W); 
  Ialp = 0;

  for (count = 1; count <= LOOP; count++){
    C = IPi * L * W + Ih * SymHull(Abs(Ialp)) * Phi;     
    Y = AbsDInv * C;

    U = u + W;
    sup_norm = Sqrt(Abs(U * (L * U)));
    SUP_norm = Abs(sup_norm);

    F_ort = Ih / IPi * (IPi * SUP_norm + Ialp + 
		       (IPi - 1.0) * sin_norm);

    if (Y < W && Sup(F_ort) <= Inf(Ialp)) {
      cout << "Iteration Stop --- 反復回数" << count << endl << endl;
      for (i = 1; i <= N; i++)
	cout << i << "  " << U(i) << "   " << Exact(i) << endl;

      WIDTH = 0.0;
      for (i = 1; i <= N; i++) {
	if ((Hull(Exact(i)) <= U(i)) == FALSE)
	  cout << i << " : 厳密解が含まれていません!!" << endl;

	width = Diam(U(i));
	if (WIDTH < width)
	  WIDTH = width;
      }	
      cout << "幅の最大 = " << WIDTH << endl;
      cout << "球の部分の半径 ; alpha = " << Abs(Ialp) << endl;
      exit(0);
    }
    W = (1.0 + DELTA) * Y;
    Ialp = (1.0 + DELTA) * F_ort;
  }
}