はさみうち法のプログラム（C言語）について

質問

調べたりしながら私なりに書いてみたのですが上手くいきません。
「ans=1.048404」とならなければいけないのに
「nan」で返ってきてしまうのです。
どこをどの様に変えたら良いのでしょうか・・・？
※課題の締め切りが近いので,早急にご回答宜しくお願い致します！

以下、はさみうちのプログラム↓
#include<stdio.h>
#include<math.h>

/*初期値x0,*/
/*f(x)を定義*/
double f(double x){
return sqrt(x*x+10.0) - 10.0*exp(-x*x);
}

/*はさみうち法の関数*/
double hasamiuchi(double x0,double x1,double eps,double imax){
double x2; /*x軸との交点*/
int i; /*繰り返しの回数*/

if(f(x0)*f(x1)>0){/*x0とx1が同符号の場合x1をx0とする*/
x0=x1;
}
for(i=0;i<imax;i++){
x2=(x1*f(x0)-x0*f(x1))/(f(x0)-f(x1)); /*x0とx1で出来る直線とx軸の交点x2を求める*/
x1=x2; /*x2の値を次のループにおけるx1とする*/
}
if(fabs(x1-x0)<eps) return x2; /*収束したら解x2を返す*/
return 0.0/0.0; /*収束しなかったらnanを返す*/
}

int main(void){
printf("ans = %f\n", hasamiuchi(2.0,1.0,1e-6,20));/*値を代入してはさみうちの実行*/
return 0;
}

回答No.5

anmochi · ベストアンサー

問題の箇所を質問文のプログラムでズバリ言うと
if(fabs(x1-x0)<eps) return x2; /*収束したら解x2を返す*/
この行だ。

はさみうち法は２つ方法があって、それぞれ何が何に収束していくかが違う。
両方が移動する方法と片方が移動する方法で、前者は２点の距離が0に、後者は移動する方の点x'におけるf(x')が0に収束する。で、あなたのプログラムは計算ロジックは後者なのに収束判定に前者のものを使っている。

～～～～～～

余談ではあるけど、移動するのは片方だけでいいんだけど、どっちを移動させるかは自動でより適切な方が選ばれるようにした方が良いだろう。あなたのプログラムは、x1が移動する方だと20回では1e-6以内に収まらない。x0の方だと5回目で0.000000に収束する。

お礼

なるほど！！
簡潔でとても分かりやすかったです！
ありがとうございました！

回答No.4

f272

http://www-it.sci.waseda.ac.jp/teachers/w405201/OLD/CPR2-2010/cprogram10.pdf
これのリスト26 falseposition.cでは
if ( fx*f0 > 0 ) x0 = xp;
と書いてあるよね。
また収束の判定についても
if ( fabs(x-xp) < DBL_EPSILON || fabs(fx) < DBL_EPSILON ) break;
となってるよね。

回答No.3

hashioogi

求める点はどこだかわからないけど決まっている。
その点を二人の人がいて一人は上から下へ、もう一人の人が下から上へ追い込んで行って、二人の距離が十分小さくなったときに挟み込めたと判断するのが挟み込みだと思うが、プログラムを見ると片方の人（x1）しか移動していない。これだとx0が移動していないからいつになっても挟み込めない。と思う。

お礼

なるほど！
確かにそうですね！！
ご指摘,ありがとうございます

回答No.2

επιστημη（@episteme）

初期値 : x0, x1 (x0 < x1) の範囲で f(x) == 0 となるxを求める

[1] xm = (x0 + x1) /2
[2.0] x１-x０ が十分小さければ xm が解である → おわり
[2.1] 解が x0 と　xm の間にあるならば x1 = xm
[2.2] 解が xm と x1 の間にあるならば x0 = xm
[2.3] どちらでもなければ失敗 → おわり
[3] [1]に戻る

補足

ご回答ありがとうございます
ですが,[1]の式を使って方程式の解を求めるのは,二分法の時だけでなく,
はさみうちもなのでしょうか・・・？

回答No.1

Tacosan

「はさみうち法」の原理は理解できていますか?

補足

http://www-it.sci.waseda.ac.jp/teachers/w405201/OLD/CPR2-2010/cprogram10.pdf
http://shiotani.cn/shiotanicnb4a/20111212cpptwo.html

以上の2つのページを基に
　(1)x0を基準点として既定のx1とx0を通る直線を考える
　(2)(1)の直線とx軸とが交わる点を次のx1とする
　(3)|x1-x0|<εとなるまで(1)(2)を繰り返す
　　|x1-x0|<εとなる時,x1が解
と考えるのかなと思っていました。
解釈が違っているのかもしれません・・・
どの様に考えれば良いのか教えていただけませんか？