ある行から下が実行されません。

[unittk501st]

85行目以下が実行されません。
試しに　setbuf(stdout, NULL);　を使ってみたのですがダメでした…
どこを直せば良いでしょうか？ご教示下さい。

コード:

#include <stdio.h>
#include <math.h>

double func1(double x){
  double y;
  y=x+1/2-sqrt(1-x*x);
  return y;
}

double func2(double x){
  double y;
  y=x+1/2+sqrt(1-x*x);
  return y;
}

int i, j;/*グローバル変数*/

double NewtonRaphson(double (*func_p)(double), double init, double eps, int imax) {/*ニュートン法*/
  double x, x_n, f, df, h;
  x=init;
  i=0;
  h=1e-4;
  while (i<imax) {
    f=func_p(x);
    df=(func_p(x+h)-func_p(x))/h;
    x_n=x - f/df;/*ニュートン反復により次の点x_n を計算する部分*/
    if (fabs(x_n - x)<eps){/* 収束条件を満たしたら解を返す*/
      return x_n;
    }
    i++;/*カウンター*/
    x=x_n;
  }
  return 0.0/0.0; /* NaN を返す*/
}

double Bisection(double (*func_p)(double), double a, double b, double eps){/*二分法*/
  double c;
  c = (a+b)/2;
  j=1;
  while(fabs(a-b)>=2*eps){/*収束前のときループ*/
    c = (a+b)/2;
    j++;/*カウンター*/
    if(func_p(c)>0){/*ｂとf(c)が同符号のとき*/
      c = b;
    }
    if(func_p(c)<0){/*aとf(c)が同符号のとき*/
      c = a;
    }
    if(func_p(c)==0.0){/*f(c)=0のとき*/
      return c;
    }
  }
  return c;
}



int main(void){
  double x;
setbuf(stdout, NULL);
  printf("Newton-Raphson method:\n");/*ニュートン法*/

  printf("Upper part: ");/*関数１*/
  x = NewtonRaphson(func1, 0, 1e-6, 20);
  if (!isnan(x)){
	  printf("(%d times):( %+.6f, %+.6f )\n", i+1, x, func1(x));
  }else{
	  printf("not converged\n");
  }
  x=0;
  printf("Lower part: ");/*関数２*/
  x = NewtonRaphson(func2, 0, 1e-6, 20);
  if (!isnan(x)){
	  printf("(%d times):( %+.6f, %+.6f )\n", i+1, x, func2(x));
  }else{
	  printf("not converged\n");
  }
  x=0;

  printf("Bisection method:\n"); /*二分法*/
  printf("Upper part: "); /*関数１*/

/*ここから下が実行されません*/

  x = Bisection(func1, -1, 1, 1e-6);
  printf("(%d times):( %+.6f, %+.6f )\n", j, x, func1(x));
  x=0.0;
  printf("Lower part: ");/*関数２*/
  x = Bisection(func2, -1, 1, 1e-6);
  printf("(%d times):( %+.6f, %+.6f )\n", j, x, func2(x));
  x=0.0;


return 0;
}

[h2so5]

「実行されない」とはどういう状況ですか？
途中でプログラムが終了するのか、それともブレークポイントを入れても一時停止しないのか、単に出力が無いのかどれでしょう。

[unittk501st]

返信ありがとうございます。

「実行されない」とは、「プログラムが終了せず85行目以下が表示されない」という意味です。

[h2so5]

それなら、Bisection関数の中で無限ループしている可能性が高いですね。

[box]

コード:

  y=x+1/2-sqrt(1-x*x);
  y=x+1/2+sqrt(1-x*x);

この2行、1/2は整数どうしの割り算ですから結果はゼロです。つまり、

コード:

  y=x-sqrt(1-x*x);
  y=x+sqrt(1-x*x);

こう書いたのと同じです。
そうなることをお望みでしょうか？

[box]

85行目以下が実行されません。

コード:

    if(func_p(c)==0.0){/*f(c)=0のとき*/
      return c;
    }

ここの判定において、「厳密に」0.0と比較しなければなりませんか？

[みけCAT]

二分探索の終了条件をEPSを用いて判断すると終了しないことがあるため、
「値が変化しなくなったら終了」という条件をおすすめします。
コード例（平方根を求める）

コード:

double mysqrt(double target) {
	double left,right,mid,now;
	left=0;right=(target<=1?1:target);
	while(1) {
		mid=(left+right)/2;
		now=mid*mid;
		if(now<target) {
			if(left==mid)break;
			left=mid;
		} else {
			if(right==mid)break;
			right=mid;
		}
	}
	return mid;
}

この場合

コード:

double Bisection(double (*func_p)(double), double a, double b, double eps){/*二分法*/
  double c;
  c = (a+b)/2;
  j=1;
  while(fabs(a-b)>=2*eps){/*収束前のときループ*/
    c = (a+b)/2;
    j++;/*カウンター*/
    if(func_p(c)>0){/*ｂとf(c)が同符号のとき*/
      if(c==b)return c;
      c = b;
    }
    if(func_p(c)<0){/*aとf(c)が同符号のとき*/
      if(c==a)return c;
      c = a;
    }
    if(func_p(c)==0.0){/*f(c)=0のとき*/
      return c;
    }
  }
  return c;
}

とするということです。

重要
Bisection関数内では、条件分岐のあとcではなくaやbを更新しないといけません。

[box]

85行目以下が実行されません。

コード:

    if(func_p(c)>0){/*ｂとf(c)が同符号のとき*/
      c = b;
    }

ここで、c の値を上書きする可能性がありますね。その上書きした c の値を

コード:

    if(func_p(c)<0){/*aとf(c)が同符号のとき*/
      c = a;
    }

ここで本当に使っていいのかな？という疑問があります。
本当は

コード:

    else if(func_p(c)<0){/*aとf(c)が同符号のとき*/
      c = a;
    }

だったりするのではないでしょうか。

[box]

どうやら、私の回答はいいかげんな内容みたいですので、
他の回答者さんの助言に従ってください。

[みけCAT]

大事なことなので二度言います。

コード:

if(func_p(c)>0){/*ｂとf(c)が同符号のとき*/
  if(c==b)return c;
  c = b;
}
if(func_p(c)<0){/*aとf(c)が同符号のとき*/
  if(c==a)return c;
  c = a;
}

ここでcの値を更新しても

コード:

c = (a+b)/2;

でもとのcの値に戻り、無限ループが発生します。
本当は、cではなくaやbの値を更新するべきではありませんか？