C..6 二分法による方程式の解法

二分法で方程式の解を求める計算をしてみよう。

test5.cpp

// test5.cpp
// g++ -I/usr/local/include/profil -o test5 test5.cpp -lProfil -lBias -llr

#include <iostream>
#include <iomanip>
#include <Interval.h>
#include <Functions.h>
using namespace std;

INTERVAL f(INTERVAL x)
{
  return x - Cos(x);
}

int DifferentSign(INTERVAL x, INTERVAL y)
{
  return ((Inf(x) > 0 && Sup(y) < 0) || (Sup(x) < 0 && Inf(y) > 0));
}

int DefiniteSign(INTERVAL x)
{
  return (Inf(x) > 0 || Sup(x) < 0);
}

void Check(REAL x)
{
  REAL a, b;
  INTERVAL fa, fb;
  if (DefiniteSign(f(Hull(x)))) {
    cout << "あれ、定符号ですけど";
    return;
  }
  a = x; b = x; fa = f(a); fb = f(b);
  while (!DefiniteSign(fa)) {
    a = Pred(a); fa = f(Hull(a));
    cout << "f(" << a << ")=" << fa << endl;
  }
  cout << "定符号になった。" << endl;
  while (!DefiniteSign(fb)) {
    b = Succ(b); fb = f(Hull(b));
    cout << "f(" << b << ")=" << fb << endl;
  }
  cout << "定符号になった。" << endl;
  INTERVAL I(a, b);
  cout << I << "に解がある" << endl;
  cout << "幅=" << Diam(I) << endl;
}

int main()
{
  REAL a, b, c;
  INTERVAL Ia, Ib, Ic, fa, fb, fc;

  cout << setprecision(18);
  a = 0; b = 1;
  fa = f(Hull(a)); fb = f(Hull(b));
  if (0 <= fa) {
    cout << "f(" << a << ")=" << fa << "は 0 を含む" << endl;
    Check(a);
  }
  if (0 <= fb) {
    cout << "f(" << b << ")=" << fb << "は 0 を含む" << endl;
    Check(b);
  }
  if (DifferentSign(fa, fb)) {
    while (1) {
      cout << "[" << a << "," << b << "] にある, 幅="
           << b-a << endl;
      c = (a + b) / 2; fc = f(c);
      if (0 <= fc) {
        cout << "f(" << c << ")=" << fc << "は 0 を含む" << endl;
        Check(c);
        break;
      }
      if (DifferentSign(fc, fb)) {
        /* [c,b] */
        a = c; fa = fc;
      }
      else if (DifferentSign(fa, fc)) {
        /* [a,c] */
        b = c; fb = fc;
      }
      else {
        cout << "error!" << endl;
        break;
      }
    }
  }
  else {
    cout << "f(" << a << ") と f(" << b << ") は異なる符号とは判定できない"
         << endl;
  }
  return 0;
}

test5 の結果

[katsurada-no-MacBook-Air-4:~/work/tmp] mk% ./test5
[0,1] にある, 幅=1
[0.5,1] にある, 幅=0.5
[0.5,0.75] にある, 幅=0.25
[0.625,0.75] にある, 幅=0.125
[0.6875,0.75] にある, 幅=0.0625
[0.71875,0.75] にある, 幅=0.03125
[0.734375,0.75] にある, 幅=0.015625
[0.734375,0.7421875] にある, 幅=0.0078125
[0.73828125,0.7421875] にある, 幅=0.00390625
[0.73828125,0.740234375] にある, 幅=0.001953125
[0.73828125,0.7392578125] にある, 幅=0.0009765625
[0.73876953125,0.7392578125] にある, 幅=0.00048828125
[0.739013671875,0.7392578125] にある, 幅=0.000244140625
[0.739013671875,0.7391357421875] にある, 幅=0.0001220703125
[0.73907470703125,0.7391357421875] にある, 幅=6.103515625e-05
[0.73907470703125,0.739105224609375] にある, 幅=3.0517578125e-05
[0.73907470703125,0.7390899658203125] にある, 幅=1.52587890625e-05
[0.73908233642578125,0.7390899658203125] にある, 幅=7.62939453125e-06
[0.73908233642578125,0.739086151123046875] にある, 幅=3.814697265625e-06
[0.739084243774414063,0.739086151123046875] にある, 幅=1.9073486328125e-06
[0.739084243774414063,0.739085197448730469] にある, 幅=9.5367431640625e-07
[0.739084720611572266,0.739085197448730469] にある, 幅=4.76837158203125e-07
[0.739084959030151368,0.739085197448730469] にある, 幅=2.384185791015625e-07
[0.739085078239440918,0.739085197448730469] にある, 幅=1.1920928955078125e-07
[0.739085078239440918,0.739085137844085694] にある, 幅=5.9604644775390625e-08
[0.739085108041763306,0.739085137844085694] にある, 幅=2.98023223876953125e-08
[0.7390851229429245,0.739085137844085694] にある, 幅=1.49011611938476563e-08
[0.739085130393505097,0.739085137844085694] にある, 幅=7.45058059692382813e-09
[0.739085130393505097,0.739085134118795395] にある, 幅=3.72529029846191407e-09
[0.739085132256150246,0.739085134118795395] にある, 幅=1.86264514923095704e-09
[0.739085133187472821,0.739085134118795395] にある, 幅=9.31322574615478516e-10
[0.739085133187472821,0.739085133653134108] にある, 幅=4.65661287307739258e-10
[0.739085133187472821,0.739085133420303464] にある, 幅=2.32830643653869629e-10
[0.739085133187472821,0.739085133303888143] にある, 幅=1.16415321826934815e-10
[0.739085133187472821,0.739085133245680482] にある, 幅=5.82076609134674073e-11
[0.739085133187472821,0.739085133216576651] にある, 幅=2.91038304567337037e-11
[0.739085133202024736,0.739085133216576651] にある, 幅=1.45519152283668519e-11
[0.739085133209300694,0.739085133216576651] にある, 幅=7.27595761418342591e-12
[0.739085133212938672,0.739085133216576651] にある, 幅=3.63797880709171296e-12
[0.739085133214757662,0.739085133216576651] にある, 幅=1.81898940354585648e-12
[0.739085133214757662,0.739085133215667157] にある, 幅=9.09494701772928238e-13
[0.739085133214757662,0.739085133215212409] にある, 幅=4.54747350886464119e-13
[0.739085133214985036,0.739085133215212409] にある, 幅=2.2737367544323206e-13
[0.739085133215098722,0.739085133215212409] にある, 幅=1.1368683772161603e-13
[0.739085133215155566,0.739085133215212409] にある, 幅=5.68434188608080149e-14
[0.739085133215155566,0.739085133215183987] にある, 幅=2.84217094304040075e-14
[0.739085133215155566,0.739085133215169777] にある, 幅=1.42108547152020038e-14
[0.739085133215155566,0.739085133215162671] にある, 幅=7.10542735760100186e-15
[0.739085133215159118,0.739085133215162671] にある, 幅=3.55271367880050093e-15
f(0.739085133215160895)=[-0.0000000000000008,0.0000000000000012]は 0 を含む
f(0.739085133215160784)=[-0.0000000000000009,0.0000000000000010]
f(0.739085133215160673)=[-0.0000000000000010,0.0000000000000006]
f(0.739085133215160562)=[-0.0000000000000015,0.0000000000000005]
f(0.739085133215160451)=[-0.0000000000000016,0.0000000000000004]
f(0.73908513321516034)=[-0.0000000000000017,0.0000000000000003]
f(0.739085133215160229)=-0.00000000000000[02,18]
定符号になった。
f(0.739085133215161006)=[-0.0000000000000007,0.0000000000000013]
f(0.739085133215161117)=[-0.0000000000000003,0.0000000000000014]
f(0.739085133215161228)=[-0.0000000000000002,0.0000000000000018]
f(0.739085133215161339)=0.00000000000000[00,19]
f(0.73908513321516145)=0.00000000000000[01,20]
定符号になった。
0.73908513321516[02,15]に解がある
幅=1.2212453270876722e-15
[katsurada-no-MacBook-Air-4:~/work/tmp] mk%