f(x) = 0の解をもとめる

２分法を使った、f(x)=0 となるこの x の値を求めるのに参考になるソースプログラムを探してます。(>_<)
直接回答でも結構ですので、お分かりになる方は、
よろしく御願いいたします。

No.5 ベストアンサー 回答者： tatsu99 回答日時： 2003/09/16 22:58

以下のようにしてください。

------------------------------------------
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
//２分法による解の算出
double func (double x){
//f(x)= x*x -1
//printf("func:x=%lf\n",x);
return x*x -1;
}

int bisection_method(double (*Fx)(double x),double *kai,double start1,double end1,int maxTime,double band)
{
double low = start1;
double high = end1;
double mid;
double fx,fx1,fx2;
int Count = 0;
*kai = -1;
// maxTime は正の値であることが前提
// Band は正の実数であることが前提
//下限でFx(x)を求める
fx1 = (*Fx)(low);
if (fabs(fx1) <= band){
*kai = low;
return 0;
}
//上限でFx(x)を求める
fx2 = (*Fx)(high);
if (fabs(fx2) <= band){
*kai = high;
return 0;
}
if (fx1 * fx2 > 0.0){
//解なし
//なぜならFx(low)とFx(high)が同じ符号だから。つまりFx(?)=0になるポイントがないから
return -1;
}
//見出しを印字
printf("loop x1 x3 x2 f(x3)\n");
//以前に計算したのは、計算回数には含めないものとする
for (Count =0; Count < maxTime; Count++){
mid = (low + high)/2.0;
fx = (*Fx)(mid);
printf("%d %lf %lf %lf %lf\n",Count,low,mid,high,fx);
if (fabs(fx) <= band){
*kai = mid;
return Count;
}
if (fx1 * fx > 0.0){
//fxの符号がfx1と同じ場合
low = mid;
}else{
//fxの符号がfx1と異なる場合
high = mid;
}
}
//実行回数オーバー
return -100;
}
int main (int argc, char **argv)
{
double x;
// int loop = bisection_method(&func, &x, -10, 10, 100, 0.001);
int loop = bisection_method(&func, &x, -0, 1000, 100, 0.001);
if (loop >= 0){
printf("func(%f)=%f loop=%d\n", x, func(x), loop);
}else{
printf("解なし\n");
}
return 0;
}
---------------------------------------
funcの解は-1と+1ですので、この２点が含まれないように
開始と終わりを指定してください
例の-10,10は解なしとなります。
0,1000は解ありとなります。

この回答へのお礼

誠にありがとうございましたっ！　まさかここまでの形でお答えいただけるとは思いもしませんでした。すごーいっの一言です。自分なら１日はたっぷりかかるだろうプログラムもすらすら書ける人がうらやましいです。動作も確認済みです。（今回は甘えすぎなのは反省しています。。）
おそらくは貴重なお時間をお使いいただいて、わざわざ作成、実行確認と最後までご回答くださったことに感謝いたします。回答者さまのお名前は忘れません。いつかポイント以外の形でお返しできればと思う所存です。

お礼日時：2003/09/17 20:49

No.4 回答者： tatsu99 回答日時： 2003/09/16 18:30

No2の続きです。

>>int loop = bisection_method(&func, &x, -10, 10, 100, 0.001);
とありますが、bisection_methodの内部のソースを公開して下さい。またこの関数の外部仕様（この関数の機能及び各々のパラメータの意味（入力か出力かも含めて））を公開して下さい。
>>？？？(なにかの処理)
とは何でしょうか。
>>printf("func(%f)=%f loop=%d\n", x, func(x), loop);
は何の為に行っているのでしょうか。
>>せめて最後の行くらい出れば・・・・
この意味は何でしょうか。
func(0.999756)=-0.000488 loop=13
を印字した後、return 0；のステートメント迄
実行されていないということでしょうか。
>2 0.000000 2.500000 5.000000 5.250000
>3 0.000000 1.250000 2.500000 0.562500
>4 0.000000 0.625000 1.250000 -0.609375
>5 0.625000 0.937500 1.250000 -0.121094
これは、どこで出力していますか。
また、この各々の値の意味はなんですか。
上記の件が不明ですので、この資料のみでは、
判断が出来ません。

この回答への補足

bisection_method関数については、これそのものを作るのが目的です。ですからソースがありません。
main関数については（なにかの処理）というのはなくして
以下のように決め打ちで訂正します。


int main() {

　double x;
　int loop = bisection_method(&func, &x, -10, 10, 100, 0.001);
　printf("func(%f)=%f loop=%d\n", x, func(x), loop);
　return 0;

}

ちなみにfunc関数というのはこのように定義されて
います。

double func (double x){

　return x*x -1;

}



-------------------------------------------
次に、
printf("func(%f)=%f loop=%d\n", x, func(x), loop); については、main で一番下の

func(0.999756)=0-0.000488 loop=13

の部分を出力してます。求めるのが x3 とf(x3) 、それとループ回数というわけです。

bisection_method()関数の中で
2 0.000000 2.500000 5.000000 5.250000
3 0.000000 1.250000 2.500000 0.562500
4 0.000000 0.625000 1.250000 -0.609375
5 0.625000 0.937500 1.250000 -0.121094

の一覧を出力し、最終的にmain に帰ってきて、さきほどの一行これを func(0.999756)=0-0.000488 loop=13
出すという仕様です。
なんで一覧を出力しているのかというと、徐々に間を
詰めて様子を確認したいからです。

bisection_method の引数の説明ですが、

　　　　　　　　解を求める
　　　　　　　　func関数のポインタ　　　　細大繰返し回数
　　　　　　　　　　　　↓　　　　　　　　　　　↓
bisection_method (&func, &x, -10, 10, 100, 0.001);
　　　　　　　　　　　　　　　　　↑　↑　↑　　　　　↑
　　　　　　　　　　　　　　　　　解　xの初期値　　　収束判定値


収束判定値が |f(x)| < error となったとき、解が得られたとし、解が求められなかった場合は負、解がまとまった
時は繰り返し回数を返す。
ということです。
私が知っていることはこれだけで・・・ようは研修の問題なんですが、ここにきてつまってしまって、質問させていただいた次第です。

補足日時：2003/09/16 20:44

No.3 回答者： episteme 回答日時： 2003/09/16 14:21

Re:#1

> f(x1)<0, f(x2)>0 のとき、x3=(x1+x2)/2 として f(x3) を計算して
> f(x3)<0 なら x1=x3, f(x3)>0 なら x2=x3 として x1とx2の差を徐々に
> 縮めていくアルゴリズムのようです。。

それがわかっているなら、そのままコードに落とすだけではないかと。

double low = /*初期値*/
double high = /*初期値*/
while ( 適当な条件 ) {
double mid = (low+high)/2.0;
/*
f(low), f(mid), f(high) を求め、
符号の反転している側(low,mid)または(mid,high)
を新たな(low,high)とする。
*/
}

この回答へのお礼

ありがとうございます。
Cをはじめたのが、昨日今日って状態なので、お手数かけてすみません。
教えていただいたとおりちょっとやってみます。

お礼日時：2003/09/16 14:55

No.2 回答者： tatsu99 回答日時： 2003/09/15 15:12

前提として以下の条件が必要です。

開始点：ａと終端：ｂの間に解が１つだけ存在すること。
f(a)～f(b)まで連続してf(x)が定義されていること
以下は、開始点：ａと終端ｂ：が与えられた時に、
f(x)=0となるxを求める方法です。
f(x)=A*x*x + B*x + C (A=4 B=6 C=-10）としてあります。
--------------------------------
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
//２分法による解の算出
double Fx(double x)
{
double A = 4.0;
double B = 6.0;
double C = -10.0;
double fx = A*x*x + B*x + C;
return fx;
}
int main (int argc, char **argv)
{
double a = atof(argv[1]);
double b = atof(argv[2]);
double c;
double fx,fx1,fx2;
double low = a;
double high = b;
int ctr = 0;
//始点ａと終端ｂを印字
printf ("a=%lf b= %lf\n",a,b);
fx1 = Fx(a);
fx2 = Fx(b);
if (fx1 == 0.0){
printf("解=%lf\n",a);
return(0);
}
if (fx2 == 0.0){
printf("解=%lf\n",b);
return(0);
}
if (fx1 * fx2 > 0){
printf("解無し\n");
return(0);
}
while(1){
ctr++;
c = (low + high)/2.0;
fx = Fx(c);
if (fabs(fx) < 0.000000001){
printf("実行回数=%d 解=%lf\n",ctr,c);
return(0);
}
if (fx1 * fx > 0.0){
low = c;
}else{
high = c;
}
if (ctr > 10000){
printf("実行回数オーバー\n");
return(0);
}
}
return 0;
}
--------------------------------------

この回答へのお礼

ソースまで見せていただき本当にありがとうございます。
今日、ようやく実機にて確かめることができました。

最終的な結果としては以下のような形が望ましいのですが、
なにかよいアドバイスはないでしょうか。
せめて最後の行くらい出れば・・・・
（ずうずうしいお願いなのでもしお手が空いてればで結構です。）


double func(double x) {
return x * x - 1;
}


int main() {
double x;
int loop = bisection_method(&func, &x, -10, 10, 100, 0.001);


？？？(なにかの処理)

printf("func(%f)=%f loop=%d\n", x, func(x), loop);


？？？(なにかの処理)

return 0;
}


loop x1 x3 x2 f(x3)
0 -10.000000 0.000000 10.000000 -1.000000
1 0.000000 5.000000 10.000000 24.000000
2 0.000000 2.500000 5.000000 5.250000
3 0.000000 1.250000 2.500000 0.562500
4 0.000000 0.625000 1.250000 -0.609375
5 0.625000 0.937500 1.250000 -0.121094
6 0.937500 1.093750 1.250000 0.196289
7 0.937500 1.015625 1.093750 0.031494
8 0.937500 0.976563 1.015625 -0.046326
9 0.976563 0.996094 1.015625 -0.007797
10 0.996094 1.005859 1.015625 0.011753
11 0.996094 1.000977 1.005859 0.001954
12 0.996094 0.998535 1.000977 -0.002928
13 0.998535 0.999756 1.000977 -0.000488
func(0.999756)=-0.000488 loop=13

お礼日時：2003/09/16 14:52

No.1 回答者： episteme 回答日時： 2003/09/15 15:06

二分法のアルゴリズムはおわかりですか？

この回答への補足

f(x1)<0, f(x2)>0 のとき、x3=(x1+x2)/2 として f(x3) を計算して
f(x3)<0 なら x1=x3, f(x3)>0 なら x2=x3 として x1とx2の差を徐々に
縮めていくアルゴリズムのようです。。
x1 とか f とかで考えると難しいから、それぞれ変数と考えるのだと・・・

補足日時：2003/09/16 14:01